format: text

text:
# Macaulay 2 and the Geometry of Schemes

## Distinguished open sets

We begin with a simple example involving the definition of an affine scheme; see section I.1.4 in Eisenbud and Harris (1999). This example also indicates some of the subtleties involved in working with arithmetic schemes.

###### Question 1.

Let $`R=[x,y,z]`$ and $`X=\mathrm{Spec}(R)`$; in other words, $`X`$ is affine $`3`$-space over the integers. Let $`f=x`$ and consider the basic open subset $`X_fX`$.

1. If $`e_1=x+y+z`$, $`e_2=xy+xz+yz`$ and $`e_3=xyz`$ are the elementary symmetric functions then the set $`\left\{X_{e_i}\right\}_{1i3}`$ is an open cover of $`X_f`$.
2. If $`p_1=x+y+z`$, $`p_2=x^2+y^2+z^2`$ and $`p_3=x^3+y^3+z^3`$ are the power sum symmetric functions then $`\left\{X_{p_i}\right\}_{1i3}`$ is NOT an open cover of $`X_f`$.

###### Solution.

By Lemma I-6 in Eisenbud and Harris (1999), it suffices to show that $`e_1`$, $`e_2`$ and $`e_3`$ generate the unit ideal in $`R_f`$. This is equivalent to showing that $`x^m`$ belongs to the $`R`$-ideal $`e_1,e_2,e_3`$ for some $`m`$. In particular, the saturation $`\left(e_1,e_2,e_3:x^{\mathrm{}}\right)`$ is the unit ideal if and only if $`\left\{X_{e_i}\right\}_{1i3}`$ is an open cover of $`X_f`$. Macaulay2 allows us to work with homogenous ideals over $``$ and we obtain:

```
    i1 : R = ZZ[x, y, z];

```

```
    i2 : elementaryBasis = ideal(x+y+z, x*y+x*z+y*z, x*y*z);

```

```
    o2 : Ideal of R

```

```
    i3 : saturate(elementaryBasis, x)

```

```
    o3 = ideal 1

```

```
    o3 : Ideal of R

```

Similarly, to prove that $`\left\{X_{p_i}\right\}_{1i3}`$ is not an open cover of $`X_f`$, it is enough to show that $`\left(p_1,p_2,p_3:x^{\mathrm{}}\right)`$ is not the unit ideal. We compute this saturation:

```
    i4 : powerSumBasis = ideal(x+y+z, x^2+y^2+z^2, x^3+y^3+z^3);

```

```
    o4 : Ideal of R

```

```
    i5 : saturate(powerSumBasis, x)

```

```
                                2           2
    o5 = ideal (6, x + y + z, 2y  - y*z + 2z , 3y*z)

```

```
    o5 : Ideal of R

```

However, working over the field $``$, we find that $`\left(p_1,p_2,p_3:x^{\mathrm{}}\right)`$ is the unit ideal.

```
    i6 : S = QQ[x, y, z];

```

```
    i7 : powerSumBasis = ideal(x+y+z, x^2+y^2+z^2, x^3+y^3+z^3);

```

```
    o7 : Ideal of S

```

```
    i8 : saturate(powerSumBasis, x)

```

```
    o8 = ideal 1

```

```
    o8 : Ideal of S

```

## Irreducibility

The study of complex semisimple Lie algebras gives rise to an important family of algebraic varieties called nilpotent orbits. To illustrate one of the properties appearing in section I.2.1 of Eisenbud and Harris (1999), we examine the irreducibility of a particular nilpotent orbit.

###### Question 2.

Let $`X`$ be the set of complex $`3\times 3`$ matrices which are nilpotent. Show that $`X`$ is an irreducible algebraic variety.

###### Solution.

A $`3\times 3`$ matrix $`M`$ is nilpotent if and only if its minimal polynomial divides $`𝖳^k`$, for some $`k`$. Since each irreducible factor of the characteristic polynomial of $`M`$ is also a factor of the minimal polynomial, we conclude that the characteristic polynomial of $`M`$ is $`𝖳^3`$. It follows that the coefficients of the characteristic polynomial (except for the leading coefficient which is $`1`$) of a generic $`3\times 3`$ matrix define the algebraic variety $`X`$.

To show $`X`$ is irreducible over $``$, it is enough to construct a rational parameterization of $`X`$; see Proposition 4.5.6 in Cox, Little, and O’Shea (1996). To achieve this, observe that $`\mathrm{GL}_n()`$ acts on $`X`$ by conjugation. Jordan’s canonical form theorem implies that there are exactly three orbits; one for each of the following matrices:

$$N_0=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],N_1=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\text{ and }N_2=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 0& 0& 0\end{array}\right].$$

Each orbit is defined by a rational parameterization, so it suffices to show that the closure of the orbit containing $`N_2`$ is the entire variety $`X`$. In Macaulay2, this calculation can be done as follows:

```
    i1 : S = QQ[t, y_0..y_8, a..i, MonomialOrder => Eliminate 10];

```

```
    i2 : N2 = (matrix {{0,1,0},{0,0,1},{0,0,0}}) ** S

```

```
    o2 = {0} | 0 1 0 |
         {0} | 0 0 1 |
         {0} | 0 0 0 |

```

```
                 3       3
    o2 : Matrix S  <--- S

```

```
    i3 : G = genericMatrix(S, y_0, 3, 3);

```

```
                 3       3
    o3 : Matrix S  <--- S

```

To determine the entries in $`GN_2G^1`$, we use the classical adjoint to construct the inverse of the generic matrix $`G`$.

```
    i4 : adj = (G,i,j) -> (
              n := degree target G;
              m := degree source G;
              (-1)^(i+j)*det(submatrix(G, {0..(i-1), (i+1)..(n-1)},
                        {0..(j-1), (j+1)..(m-1)}))
              );

```

```
    i5 : classicalAdjoint = (G) -> (
              n := degree target G;
              matrix table(n, n, (i, j) -> adj(G, j, i))
              );

```

```
    i6 : numerators = G*N2*classicalAdjoint(G);

```

```
                 3       3
    o6 : Matrix S  <--- S

```

```
    i7 : D = det(G);

```

```
    i8 : M = genericMatrix(S, a, 3, 3);

```

```
                 3       3
    o8 : Matrix S  <--- S

```

The entries in $`GN_2G^1`$ give a rational parameterization of the orbit generated by $`N_2`$. We give an “implicit representation” of this variety by using elimination theory; see section 3.3 in Cox, Little, and O’Shea (1996).

```
    i9 : elimIdeal = minors(1, (D*id_(S^3))*M-numerators) + ideal(1 - D*t);

```

```
    o9 : Ideal of S

```

```
    i10 : closureOfOrbit = ideal(selectInSubring(1, gens gb elimIdeal));

```

```
    o10 : Ideal of S

```

Finally, we check that the closure of this orbit is equal to $`X`$ scheme-theoretically.

```
    i11 : X = ideal submatrix((coefficients({0},
                       det(M - t*id_(S^3))))_1, {1,2,3})

```

```
    o11 = ideal (a + e + i, b*d - a*e + c*g + f*h - a*i - e*i,
               - c*e*g + b*f*g + c*d*h - a*f*h - b*d*i + a*e*i)

```

```
    o11 : Ideal of S

```

```
    i12 : closureOfOrbit == X

```

```
    o12 = true

```

## Singular Points

Section I.2.2 in Eisenbud and Harris (1999) provides the definition of a singular point of a scheme. In our third question, we study the singular locus of a family of elliptic curves. Section V.3 in Eisenbud and Harris (1999) also contains related material.

###### Question 3.

Consider a general form of degree $`3`$ in $`[x,y,z]`$;

$$F=ax^3+bx^2y+cx^2z+dxy^2+exyz+fxz^2+gy^3+hy^2z+iyz^2+jz^3.$$

Give necessary and sufficient conditions in terms of $`a,\mathrm{},j`$ for the cubic curve $`\mathrm{Proj}\left([x,y,z]/F\right)`$ to have a singular point.

###### Solution.

A time consuming elimination gives the degree $`12`$ polynomial which defines the singular locus of a general form of degree $`3`$. This can be done in Macaulay2 as follows. We have not displayed the output o6, as this discriminant has $`2040`$ terms in the $`10`$ variables $`a,\mathrm{},j`$.

```
    i1 : S = QQ[x, y, z, a..j, MonomialOrder => Eliminate 2];

```

```
    i2 : F = a*x^3+b*x^2*y+c*x^2*z+d*x*y^2+e*x*y*z+
         f*x*z^2+g*y^3+h*y^2*z+i*y*z^2+j*z^3;

```

```
    i3 : partials = submatrix(jacobian matrix{{F}}, {0..2}, {0})

```

```
    o3 = {1} | 3x2a+2xyb+y2d+2xzc+yze+z2f |
         {1} | x2b+2xyd+3y2g+xze+2yzh+z2i |
         {1} | x2c+xye+y2h+2xzf+2yzi+3z2j |

```

```
                 3       1
    o3 : Matrix S  <--- S

```

```
    i4 : singularities = ideal(partials) + ideal(F);

```

```
    o4 : Ideal of S

```

```
    i5 : elimDiscr = ideal selectInSubring(1, gens gb singularities);

```

```
    o5 : Ideal of S

```

```
    i6 : elimDiscr = substitute(elimDiscr, {z => 1});

```

```
    o6 : Ideal of S

```

There is also a simple and more useful determinantal formula for this discriminant. It is a specialization of the formula (2.8) in section 3.2 in Cox, Little, and O’Shea (1998):

```
    i7 : hessian = det submatrix(jacobian ideal partials, {0..2}, {0..2});

```

```
    i8 : A = (coefficients({0,1,2}, submatrix(
                       jacobian matrix{{F}}, {0..2}, {0})))_1;

```

```
                 3       6
    o8 : Matrix S  <--- S

```

```
    i9 : B = (coefficients({0,1,2}, submatrix(
                       jacobian matrix{{hessian}}, {0..2}, {0})))_1;

```

```
                 3       6
    o9 : Matrix S  <--- S

```

```
    i10 : detDiscr = ideal det (A || B);

```

```
    o10 : Ideal of S

```

```
    i11 : detDiscr == elimDiscr

```

```
    o11 = true

```

## Fields of Definition

Schemes over non-algebraically closed fields arise in number theory. Our solution to Exercise II-6 in Eisenbud and Harris (1999) indicates one technique for working over a number field in Macaulay2.

###### Question 4.

An inclusion of fields $`KL`$ induces a map $`𝔸_L^n𝔸_K^n`$. Find the images in $`𝔸_{}^2`$ of the following points of $`𝔸_\overline{}^2`$ under this map.

1. $`x\sqrt{2},y\sqrt{2}`$;
2. $`x\sqrt{2},y\sqrt{3}`$;
3. $`x\zeta ,y\zeta ^1`$ where $`\zeta `$ is a $`5`$-th root of unity ;
4. $`\sqrt{2}x\sqrt{3}y`$;
5. $`\sqrt{2}x\sqrt{3}y1`$.

###### Solution.

The images can be determined by (1) replacing coefficients not belonging to $`K`$ with indeterminates, (2) adding the minimal polynomials of these coefficients to the given ideal in $`𝔸_\overline{}^2`$ and (3) eliminating the new indeterminates. Here are the five examples:

```
    i1 : S = QQ[a, b, x, y, MonomialOrder => Eliminate 2];

```

```
    i2 : Ia = ideal(x-a, y-a, a^2-2);

```

```
    o2 : Ideal of S

```

```
    i3 : ideal selectInSubring(1, gens gb Ia)

```

```
                        2
    o3 = ideal (x - y, y  - 2)

```

```
    o3 : Ideal of S

```

```
    i4 : Ib = ideal(x-a, y-b, a^2-2, b^2-3);

```

```
    o4 : Ideal of S

```

```
    i5 : ideal selectInSubring(1, gens gb Ib)

```

```
                 2       2
    o5 = ideal (y  - 3, x  - 2)

```

```
    o5 : Ideal of S

```

```
    i6 : Ic = ideal(x-a, y-a^4, a^4+a^3+a^2+a+1);

```

```
    o6 : Ideal of S

```

```
    i7 : ideal selectInSubring(1, gens gb Ic)

```

```
                          2    2               3    2
    o7 = ideal (x*y - 1, x  + y  + x + y + 1, y  + y  + x + y + 1)

```

```
    o7 : Ideal of S

```

```
    i8 : Id = ideal(a*x+b*y, a^2-2, b^2-3);

```

```
    o8 : Ideal of S

```

```
    i9 : ideal selectInSubring(1, gens gb Id)

```

```
                2   3  2
    o9 = ideal(x  - -*y )
                    2

```

```
    o9 : Ideal of S

```

```
    i10 : Ie = ideal(a*x+b*y-1, a^2-2, b^2-3);

```

```
    o10 : Ideal of S

```

```
    i11 : ideal selectInSubring(1, gens gb Ie)

```

```
                 4     2 2   9  4    2   3  2   1
    o11 = ideal(x  - 3x y  + -*y  - x  - -*y  + -)
                             4           2      4

```

```
    o11 : Ideal of S

```

## Multiplicity

The multiplicity of a zero-dimensional scheme $`X`$ at a point $`pX`$ is defined to be the length of the local ring $`𝒪_{X,p}`$. Unfortunately, we cannot work directly in the local ring in Macaulay2. What we can do, however, is to compute the multiplicity by computing the degree of the component of $`X`$ supported at $`p`$; see page 66 in Eisenbud and Harris (1999).

###### Question 5.

What is the multiplicity of the origin $`(0,0,0)`$ as a zero of the polynomial equations

$$x^5+y^3+z^3=x^3+y^5+z^3=x^3+y^3+z^5=0\mathrm{?}$$

###### Solution.

If $`I`$ is the ideal generated by $`x^5+y^3+z^3`$, $`x^3+y^5+z^3`$ and $`x^3+y^3+z^5`$ in $`[x,y,z]`$, then the multiplicity of the origin is

$$dim_{}\frac{[x,y,z]_{x,y,z}}{I[x,y,z]_{x,y,z}}.$$

It follows that the multiplicity is the vector space dimension of the ring $`[x,y,z]/\phi ^1(I[x,y,z]_{x,y,z})`$ where $`\phi :[x,y,z][x,y,z]_{x,y,z}`$ is the natural map. Moreover, we can express this using ideal quotients:

$$\phi ^1(I[x,y,z]_{x,y,z})=(I:(I:x,y,z^{\mathrm{}})).$$

Carrying out this calculation in Macaulay2, we obtain:

```
    i1 : S = QQ[x, y, z];

```

```
    i2 : I = ideal(x^5+y^3+z^3, x^3+y^5+z^3, x^3+y^3+z^5);

```

```
    o2 : Ideal of S

```

```
    i3 : multiplicity = degree(I : saturate(I))

```

```
    o3 = 27

```

## Flat Families

Non-reduced schemes arise naturally as the flat limit of a family of reduced schemes. Exercise III-68 in Eisenbud and Harris (1999) illustrates how a family of skew lines in $`^3`$ gives a double line with an embedded point.

###### Question 6.

Let $`L`$ and $`M`$ be the lines in $`_{k[t]}^3`$ given by $`x=y=0`$ and $`xtz=y+t^2w=0`$ respectively. Show that the flat limit as $`t0`$ of the union $`LM`$ is the double line $`x^2=y=0`$ with an embedded point of degree $`1`$ located at the point $`(0:0:0:1)`$.

###### Solution.

We find the flat limit by saturating the intersection ideal:

```
    i1 : PP3 = QQ[x, y, z, w];

```

```
    i2 : S = QQ[t, x, y, z, w];

```

```
    i3 : phi = map(PP3, S, 0 | vars PP3 );

```

```
    o3 : RingMap PP3 <--- S

```

```
    i4 : L = ideal(x, y);

```

```
    o4 : Ideal of S

```

```
    i5 : M = ideal(x-t*z, y-t^2*w);

```

```
    o5 : Ideal of S

```

```
    i6 : X = intersect(L, M);

```

```
    o6 : Ideal of S

```

```
    i7 : Xzero = trim phi substitute(saturate(X, t), {t => 0})

```

```
                      2        2
    o7 = ideal (y*z, y , x*y, x )

```

```
    o7 : Ideal of PP3

```

This is the union of a double line and an embedded point of degree $`1`$.

```
    i8 : use PP3;

```

```
    i9 : intersect(ideal(x^2, y), ideal(x, y^2, z))

```

```
                      2        2
    o9 = ideal (y*z, y , x*y, x )

```

```
    o9 : Ideal of PP3

```

```
    i10 : degree( ideal(x^2, y) / ideal(x, y^2, z))

```

```
    o10 = 1

```

## Bézout’s Theorem

Bézout’s theorem (Theorem III-78 in Eisenbud and Harris, 1999) fails without the Cohen-Macaulay hypothesis. Following Exercise III-81 in Eisenbud and Harris (1999), we illustrate this in Macaulay2.

###### Question 7.

Find irreducible closed subvarieties $`X`$ and $`Y`$ in $`^4`$ with

$`\mathrm{codim}(XY)`$ $`=`$ $`\mathrm{codim}(X)+\mathrm{codim}(Y)\text{ and}`$
$`\mathrm{deg}(XY)`$ $`>`$ $`\mathrm{deg}(X)\mathrm{deg}(Y).`$
###### Solution.

We show that the assertion holds when $`X`$ is the cone over the nonsingular rational quartic curve in $`^3`$ and $`Y`$ is a two-plane passing through the vertex of the cone. The computation is done as follows:

```
    i1 : S = QQ[a, b, c, d, e];

```

```
    i2 : quarticCone = trim minors(2,
              matrix{{a, b^2, b*d, c}, {b, a*c, c^2, d}})

```

```
                            3      2     2    2    3    2
    o2 = ideal (b*c - a*d, c  - b*d , a*c  - b d, b  - a c)

```

```
    o2 : Ideal of S

```

```
    i3 : plane = ideal(a, d);

```

```
    o3 : Ideal of S

```

```
    i4 : codim quarticCone + codim plane == codim (quarticCone + plane)

```

```
    o4 = true

```

```
    i5 : (degree quarticCone) * (degree plane)

```

```
    o5 = 4

```

```
    i6 : degree (quarticCone + plane)

```

```
    o6 = 5

```

## Constructing Blow-ups

The blow-up of a scheme $`X`$ along a subscheme $`Y`$ can be constructed from the Rees algebra associated to the ideal sheaf of $`Y`$ in $`X`$; see Theorem IV-22 in Eisenbud and Harris (1999). Gröbner basis techniques allow one to express the Rees algebra in terms of generators and relations. We demonstrate this by solving Exercise IV-43 in Eisenbud and Harris (1999).

###### Question 8.

Find the blow-up of the affine plane $`𝔸^2=\mathrm{Spec}\left(k[x,y]\right)`$ along the subscheme defined by $`x^3,xy,y^2`$.

###### Solution.

We first provide a general function which given an ideal and a list of variables returns the ideal of relations for the Rees algebra.

```
    i1 : blowUpIdeal = method();

```

```
    i2 : blowUpIdeal(Ideal, List) := (I, L) -> (
             r := numgens I;
             S := ring I;
             kk := coefficientRing S;
             n := numgens S;
             y := symbol y;
             St := kk[t, gens S , y_1..y_r, MonomialOrder => Eliminate 1];
             phi := map(St, S, submatrix(vars St, {1..n}));
             F := phi gens I;
             local J;
             J = ideal apply(1..r, j -> y_j - t*F_(0, (j-1)));
             J = ideal selectInSubring(1, gens gb J);
             if (#L < r) then error "not enough variables";
             R := kk[gens S, L];
             theta := map(R, St, 0 | vars R);
             theta J
         );

```

Applying the function to our specific case yields:

```
    i3 : S = QQ[x, y];

```

```
    i4 : I = ideal(x^3, x*y, y^2);

```

```
    o4 : Ideal of S

```

```
    i5 : blowUpIdeal(I, {A, B, C})

```

```
                              2         2          3     2
    o5 = ideal (y*B - x*C, x*B  - A*C, x B - y*A, x C - y A)

```

```
    o5 : Ideal of QQ [x, y, A, B, C]

```

## A Classic Blow-up

We consider the blow-up of the projective plane $`^2`$ at a point. Many related examples appear in section IV.2.2 of Eisenbud and Harris (1999).

###### Question 9.

Show that the following varieties are isomorphic.

1. the image of the rational map from $`^2`$ to $`^4`$ given by

$$(r:s:t)(r^2:s^2:rs:rt:st);$$
2. the blow-up of the plane $`^2`$ at the point $`(0:0:1)`$;
3. the determinantal variety defined by the $`2\times 2`$ minors of the matrix

$$\left[\begin{array}{ccc}a& c& d\\ b& d& e\end{array}\right]$$

where $`^4=\mathrm{Proj}\left(k[a,b,c,d,e]\right)`$.

This surface is called the cubic scroll in $`^4`$.

###### Solution.

We find the ideal in part (a) by elimination theory.

```
    i1 : PP4 = QQ[a, b, c, d, e];

```

```
    i2 : S = QQ[r, s, t, A..E, MonomialOrder => Eliminate 3 ];

```

```
    i3 : I = ideal(A - r^2, B - s^2, C - r*s, D - r*t, E - s*t);

```

```
    o3 : Ideal of S

```

```
    i4 : phi = map(PP4, S, (matrix{{0, 0, 0}}**PP4) | vars PP4)

```

```
    o4 = map(PP4,S,{0, 0, 0, a, b, c, d, e})

```

```
    o4 : RingMap PP4 <--- S

```

```
    i5 : surfaceA = phi ideal selectInSubring(1, gens gb I)

```

```
                                             2
    o5 = ideal (c*d - a*e, b*d - c*e, a*b - c )

```

```
    o5 : Ideal of PP4

```

We determine the surface in part (b) by constructing the blow-up of $`^2`$ in $`^2\times ^1`$ and then projecting its Segre embedding from $`^5`$ into $`^4`$. Notice that its image under the Segre map lies on a hyperplane in $`^5`$.

```
    i6 : R = QQ[t, x, y, z, u, v, MonomialOrder => Eliminate 1];

```

```
    i7 : blowUpIdeal = ideal selectInSubring(1, gens gb ideal(u-t*x, v-t*y))

```

```
    o7 = ideal(y*u - x*v)

```

```
    o7 : Ideal of R

```

```
    i8 : PP2xPP1 = QQ[x, y, z, u, v];

```

```
    i9 : psi = map(PP2xPP1, R, 0 | vars PP2xPP1);

```

```
    o9 : RingMap PP2xPP1 <--- R

```

```
    i10 : blowUp = PP2xPP1 / psi(blowUpIdeal);

```

```
    i11 : PP5 = QQ[A, B, C, D, E, F];

```

```
    i12 : segre = map(blowUp, PP5, matrix{{x*u, y*u, z*u, x*v, y*v, z*v}});

```

```
    o12 : RingMap blowUp <--- PP5

```

```
    i13 : ker segre

```

```
                                    2
    o13 = ideal (B - D, C*E - D*F, D  - A*E, C*D - A*F)

```

```
    o13 : Ideal of PP5

```

```
    i14 : theta = map( PP4, PP5, matrix{{a, c, d, c, b, e}})

```

```
    o14 = map(PP4,PP5,{a, c, d, c, b, e})

```

```
    o14 : RingMap PP4 <--- PP5

```

```
    i15 : surfaceB = trim theta ker segre

```

```
                                              2
    o15 = ideal (c*d - a*e, b*d - c*e, a*b - c )

```

```
    o15 : Ideal of PP4

```

Finally, we compute the surface in part (c) and apply a permutation of the variables to obtain the desired isomorphisms

```
    i16 : determinantal = minors(2, matrix{{a, c, d},{b, d, e}})

```

```
                                              2
    o16 = ideal (- b*c + a*d, - b*d + a*e, - d  + c*e)

```

```
    o16 : Ideal of PP4

```

```
    i17 : sigma = map(PP4, PP4, matrix{{d, e, a, c, b}});

```

```
    o17 : RingMap PP4 <--- PP4

```

```
    i18 : surfaceC = sigma determinantal

```

```
                                              2
    o18 = ideal (c*d - a*e, b*d - c*e, a*b - c )

```

```
    o18 : Ideal of PP4

```

```
    i19 : surfaceA == surfaceB

```

```
    o19 = true

```

```
    i20 : surfaceB == surfaceC

```

```
    o20 = true

```

## Fano Schemes

Our last example concerns the family of Fano schemes associated to a flat family of quadrics. We solve Exercise IV-69 in Eisenbud and Harris (1999).

###### Question 10.

Consider the one-parameter family of quadrics

$$Q=V(tx^2+ty^2+tz^2+w^2)_{k[t]}^3=\mathrm{Proj}\left(k[t][x,y,z,w]\right).$$

As the fiber $`Q_t`$ tends to the double plane $`Q_0`$, what is the flat limit of the Fano scheme $`F_1(Q_t)`$ of lines lying on these quadric surfaces ?

###### Solution.

We first make the homogeneous coordinate ring of the ambient projective $`3`$-space and the ideal of our family of quadrics.

```
    i1 : PP3overBase = QQ[t, x, y, z, w];

```

```
    i2 : Qt = ideal(t*x^2+t*y^2+t*z^2+w^2);

```

```
    o2 : Ideal of PP3overBase

```

We construct an indeterminate line in $`_{[t]}^3`$ by adding parameters $`u,v`$ and two points $`\left(A:B:C:D\right)`$ and $`\left(E:F:G:H\right)`$. The map $`\varphi `$ sends the variables to the coordinates of the general point on this line.

```
    i3 : S = QQ[t, u, v, A..H];

```

```
    i4 : phi = map(S, PP3overBase, matrix{{t}} |
              u*matrix{{A, B, C, D}}+v*matrix{{E, F, G, H}});

```

```
    o4 : RingMap S <--- PP3overBase

```

The indeterminate line is contained in our family of quadrics if and only if $`\varphi (tx^2+ty^2+tz^2+w^2)`$ vanishes identically in $`u`$ and $`v`$. Thus, we extract the coefficients of $`u`$ and $`v`$.

```
    i5 : imageOfQt = phi Qt;

```

```
    o5 : Ideal of S

```

```
    i6 : coeff = (coefficients({1,2}, gens imageOfQt))_1;

```

```
                 1       3
    o6 : Matrix S  <--- S

```

We no longer need the variables $`u`$ and $`v`$.

```
    i7 : Sprime = QQ[t, A..H];

```

```
    i8 : coeff = substitute(coeff, Sprime);

```

```
                      1            3
    o8 : Matrix Sprime  <--- Sprime

```

```
    i9 : Sbar = Sprime / (ideal coeff);

```

To move to the Grassmannian over $`[t]`$, we introduce a polynomial ring in $`6`$ new variables corresponding to the minors of the matrix $`\left[\begin{array}{cccc}A& B& C& D\\ E& F& G& H\end{array}\right]`$. The map $`\psi `$ sends the new variables $`a,\mathrm{}f`$ to the appropriate minor, regarded as elements of $`\mathrm{𝚂𝚋𝚊𝚛}`$.

```
    i10 : PP5overBase = QQ[t, a..f];

```

```
    i11 : psi = matrix{{Sbar_"t"}} | substitute(
               exteriorPower(2, matrix{{A, B, C, D}, {E, F, G, H}}), Sbar);

```

```
                     1          7
    o11 : Matrix Sbar  <--- Sbar

```

```
    i12 : fanoOfQt = trim ker map(Sbar, PP5overBase, psi);

```

```
    o12 : Ideal of PP5overBase

```

We next determine the limit as $`t`$ tends to $`0`$.

```
    i13 : fanoOfQ0 = trim substitute(saturate(
              fanoOfQt, t), {t => 0})

```

```
                                            2   2   2
    o13 = ideal (e*f, d*f, d*e, a*e + b*f, d , f , e , c*d - b*e + a*f,

```

```
                                        2    2    2
                 b*d + c*e, a*d - c*f, a  + b  + c )

```

```
    o13 : Ideal of PP5overBase

```

```
    i14 : degree( ideal(d, e, f, a^2+b^2+c^2) / fanoOfQ0)

```

```
    o14 = 2

```

We see that $`F_1(Q_0)`$ is supported on the plane conic $`d,e,f,a^2+b^2+c^2`$ and $`F_1(Q_0)`$ is not reduced — it has multiplicity two.

From section IV.3.2 in Eisenbud and Harris (1999), we know that $`F_1(Q_1)`$ is the union of two conics lying in complementary planes. We verify this as follows:

```
    i15 : fanoOfQ1 = trim substitute(saturate(
              fanoOfQt, t), {t => 1})

```

```
                             2    2    2
    o15 = ideal (a*e + b*f, d  + e  + f , c*d - b*e + a*f, b*d + c*e,

```

```
                             2    2    2                         2    2
                 a*d - c*f, c  + e  + f , b*c + d*e, a*c - d*f, b  - e ,

```

```
                             2    2
                 a*b + e*f, a  - f )

```

```
    o15 : Ideal of PP5overBase

```

```
    i16 : fanoOfQ1 == intersect(ideal(c-d, b+e, a-f, d^2+e^2+f^2),
               ideal(c+d, b-e, a+f, d^2+e^2+f^2))

```

```
    o16 = true

```

Thus, $`F_1(Q_0)`$ is the double conic obtained when the two conics in $`F_1(Q_1)`$ move together.


text:
# Implementation of the Projector Augmented Wave LDA+U Method: Application to the Electronic Structure of NiO

## I Introduction

For many materials, the density functional theory (DFT) in the local-spin-density approximation (LSDA) provides a good description of their ground-state properties. However, problems arise when the DFT-LSDA approach is applied to materials with ions that contain incomplete $`d`$\- or $`f`$-shells, such as transition-metal oxides or heavy fermion systems. For example, most transition-metal oxides are wide-gap antiferromagnetic insulators, and the DFT-LSDA predicts them to be either metals (FeO and CoO) or small-gap semiconductors (MnO and NiO). The failure of the DFT-LSDA can be traced to the mean-field character of the Kohn-Sham equations as well as to the poor description of strong correlation effects within the homogeneous electron gas. The strong correlation effects are responsible for the breakdown of the DFT-LSDA description of the electronic structure of these compounds. In order to provide a better description of these effects, the Mott-Hubbard picture has been introduced.

In the Mott-Hubbard picture of NiO, the $`d`$-$`d`$ Coulomb interaction splits the Ni $`d`$ sub-bands into the so-called lower and upper Hubbard bands. The upper Hubbard band has mostly Ni $`3d^9`$ character, while the top of the valence band is of $`3d^8`$ character, leading to a Mott-Hubbard gap of $`d`$-$`d`$ type. However, O $`1s`$ x-ray absorption as well as x-ray photoemission and bremsstrahlung isochromat spectroscopies on LixNi1-xO have shown that the additional hole has mainly oxygen character. In contrast to the Mott-Hubbard model, the energy band gap caused by the Ni $`3d`$ correlations is therefore of the charge-transfer type between the occupied oxygen $`2p`$ and the Ni $`3d`$ empty states.

On the other hand, localized approaches,  such as the local cluster scheme based on the configuration interaction method or the Anderson impurity model, in which transition-metal ions are treated like an impurity in an oxygen $`2p`$ host, predict a well-defined band gap of 5.0 eV. However, the oxygen $`2p`$ band dispersion observed in angle-resolved photoemission spectroscopy cannot be described by these methods because the lattice effects are neglected.

Several attempts have been made to include the missing correlation effects in DFT-LSDA. The generalized gradient approximation (GGA), which takes into account the radial and angular gradient corrections, can only open a small band gap. The self-interaction correction (SIC) eliminates the spurious interaction of an electron with itself from the conventional DFT-LSDA method. Compared to LSDA, the band gap and the magnetic moments are significantly increased. However, the band gap still is too small, and the SIC-LSDA method predicts a larger energy band gap for NiO than for FeO and CoO, in contradiction to experiment. The crystal-field orbital polarization introduced by Norman  to determine the magnetism and insulating band gap of transition-metal oxides is promising but underestimates both the spin magnetic moment and the band gap.

Another promising approach for correlated materials is the so-called local density approximation (LDA) plus the multi-orbital mean-field Hubbard model (LDA+U) which includes the on-site Coulomb interaction in the LSDA Hamiltonian. After adding the on-site Coulomb interaction to the LSDA Hamiltonian, the potential becomes spin- and orbital-dependent. Because a larger energy cost is associated with fluctuations of the $`d`$-occupancy, the orbital-dependent potential reduces the fluctuations of the $`d`$-occupancy, resulting in a better justification of a mean-field approach. LDA+U, although it is a mean-field approach, has the advantage of describing both the chemical bonding and the electron-electron interaction.

The question regarding the best value for the Coulomb repulsion parameter $`U`$ is, however, still under debate. The $`U`$ parameter for NiO obtained from a constrained LDA calculation is about 7 to 8 eV, and this is the value generally used in LDA+U calculations. A similar value of $`U`$ has been obtained from a constrained LDA calculation for bulk Fe, even though a much smaller value had been expected because of the metallic screening in Fe. The authors argued that the higher value could be an artifact due to the poor screening within the atomic sphere approximation (ASA), and that within a full-potential calculation a much smaller value of less than 4 eV would be expected. In contrast, an unpublished full-potential linear muffin-tin orbital method calculation by Alouani and Wills clearly shows that the value of $`U`$ for bulk Fe is even slightly larger than the ASA value. Furthermore, Bulut, Scalapino, and White showed that the renormalization of the Coulomb interaction depends on the type of model used. As LDA is not a diagrammatic method, it is not known which type of renormalization is the most appropriate for the LDA+U model.

In this work we shed light on this problem by treating the Coulomb repulsion parameter $`U`$ as adjustable parameter, and by investigating how the electronic and optical properties depend on its value. We show that for an intermediate value of $`U=5`$ eV, good agreement with the measured ground-state antiferromagnetic magnetic moment and optical properties is obtained. We also show that the O $`2p`$ character near the top of the valence states is enhanced for a larger value of $`U`$. Our calculation seems to indicate that the nature of the band gap at intermediate $`U`$ is a mixture of charge transfer and Mott-Hubbard type, and that it becomes almost purely of the charge-transfer type for higher values of $`U`$.

Our calculations are based on the projector augmented wave (PAW) method, an efficient all-electron method without shape approximations on the potential or electron density to avoid uncertainties due to the ASA approach. Based on a Car-Parrinello-like formalism, the PAW method allows complex relaxations and dynamical properties in strongly correlated systems to be studied. Our implementation of LDA+U within the PAW method is described in detail. Furthermore we discuss possible extensions of the existing method that will enhance its applicability.

The paper is organized as follows: In Section II we present those aspects of the PAW formalism that are needed for the implementation of the LDA+U method. In Section III we present and discuss the LSDA and LDA+U ground-state properties of NiO, and in Section IV we study the optical properties of NiO, namely the imaginary part of the dielectric function, and compare the results to experiment.

## II Formalism

### A PAW method

The PAW method developed by one of us  combines ideas of the pseudo-potential (PP) and the linear augmented plane wave (LAPW) methods. It is applicable to all elements of the periodic table. The nodal behavior of the wave function is correctly described and, as in the PP method, the forces on the ions are easily expressed.

In the PAW method, the all-electron (AE) crystal wave function is constructed from a pseudo (PS) wave function and atom-like functions localized near the nuclei. The PS wave function $`|\stackrel{~}{\mathrm{\Psi }}`$ coincides with the crystal AE wave function $`|\mathrm{\Psi }`$ in the interstitial region, i.e. outside the atomic regions. Inside the atomic regions $`\mathrm{\Omega }_t`$, called augmentation regions, the wave function is almost atom-like because the effect of the surrounding crystal is small. Therefore, a natural choice is to use solutions $`|\varphi _\mathrm{\Lambda }`$ of Schrödinger’s equation for the isolated atom, the so-called AE partial waves, as a basis set for the augmentation region. Here $`\mathrm{\Lambda }=\{t,\alpha ,\mathrm{},m\}`$ is a global index for the atom $`t`$, the angular momentum $`\mathrm{}`$, the magnetic quantum number $`m`$, and the index $`\alpha `$, the energy for which Schrödinger’s equation is solved.

To link the expansion in atom-like functions near the nuclei to the PS wave function, we introduce a set of auxiliary functions $`|\stackrel{~}{\varphi }_\mathrm{\Lambda }`$, so-called PS partial waves, which are centered on the atom and coincide per construction with the corresponding AE partial waves $`|\varphi _\mathrm{\Lambda }`$ outside their augmentation regions:

$$\varphi _\mathrm{\Lambda }(𝐫)=\stackrel{~}{\varphi }_\mathrm{\Lambda }(𝐫)\text{for}r\mathrm{\Omega }_t.$$
(1)
The coefficients $`c_\mathrm{\Lambda }`$ of the expansions in AE and PS partial waves are chosen such that the PS partial wave expansion $`_\mathrm{\Lambda }|\stackrel{~}{\varphi }_\mathrm{\Lambda }c_\mathrm{\Lambda }`$ cancels out the PS wave function $`|\stackrel{~}{\mathrm{\Psi }}`$ inside the augmentation region. For this purpose we introduce so-called projector functions $`\stackrel{~}{p}|`$ such that

$$\underset{\mathrm{\Lambda }}{}|\stackrel{~}{\varphi }_\mathrm{\Lambda }\stackrel{~}{p}_\mathrm{\Lambda }|=1,$$
(2)
and therefore

$$|\stackrel{~}{\mathrm{\Psi }}=\underset{\mathrm{\Lambda }}{}|\stackrel{~}{\varphi }_\mathrm{\Lambda }\stackrel{~}{p}_\mathrm{\Lambda }|\stackrel{~}{\mathrm{\Psi }}$$
(3)
for the Hilbert space spanned by the PS partial waves $`|\stackrel{~}{\varphi }_\mathrm{\Lambda }`$. Thus we identify the expansion coefficients with $`c_\mathrm{\Lambda }=\stackrel{~}{p}_\mathrm{\Lambda }|\stackrel{~}{\mathrm{\Psi }}`$. Equation (2) results in the biorthogonality condition

$$\stackrel{~}{p}_\mathrm{\Lambda }|\stackrel{~}{\varphi }_\mathrm{\Lambda }^{}=\delta _{\mathrm{\Lambda },\mathrm{\Lambda }^{}}$$
(4)
for the projector functions, which moreover are chosen to be localized within the corresponding augmentation region.

With these conditions, the AE Bloch wave function $`\mathrm{\Psi }(𝐫)`$ can be obtained from the PS wave function $`\stackrel{~}{\mathrm{\Psi }}(𝐫)`$ as

$$\mathrm{\Psi }(𝐫)=\stackrel{~}{\mathrm{\Psi }}(𝐫)+\underset{\mathrm{\Lambda }}{}[\varphi _\mathrm{\Lambda }(𝐫)\stackrel{~}{\varphi }_\mathrm{\Lambda }(𝐫)]\stackrel{~}{p}_\mathrm{\Lambda }|\stackrel{~}{\mathrm{\Psi }}.$$
(5)
The first term represents the PS wave function defined over the entire space, which is equal to the AE wave function in the interstitial region, and which is expanded in plane waves. The second term is the AE partial wave expansion, which describes the correct nodal behavior of the wave function in the augmentation region $`\mathrm{\Omega }_R^t`$ ($`rr_c^t`$). The third term eliminates the spurious contribution of the PS wave function in the augmentation region.

Note that Eq. (3) holds only approximately if the set of PS partial waves is not entirely within the augmentation regions. However, this has the advantage that only those contributions of the PS wave function will be removed that are also replaced by AE partial waves. As a result, the AE wave function converges rapidly with the number of partial waves used and, moreover, it is continuous and differentiable for every truncation of the partial-wave expansion.

Expectation values of any sufficiently local operator $`A`$ are obtained as

$$\mathrm{\Psi }|A|\mathrm{\Psi }=\stackrel{~}{\mathrm{\Psi }}|A|\stackrel{~}{\mathrm{\Psi }}+\underset{\mathrm{\Lambda },\mathrm{\Lambda }^{}}{}\stackrel{~}{\mathrm{\Psi }}|\stackrel{~}{p}_\mathrm{\Lambda }\left(\varphi _\mathrm{\Lambda }|A|\varphi _\mathrm{\Lambda }^{}\stackrel{~}{\varphi }_\mathrm{\Lambda }|A|\stackrel{~}{\varphi }_\mathrm{\Lambda }^{}\right)\stackrel{~}{p}_\mathrm{\Lambda }^{}|\stackrel{~}{\mathrm{\Psi }}.$$
(6)
Note that the double sum is diagonal in the site indices $`t,t^{}`$. This equation is exact for a complete set of PS partial waves and rapidly attains the converged result if incomplete. The PAW method provides the freedom to represent a zero operator in the form

$$0=\stackrel{~}{\mathrm{\Psi }}|B|\stackrel{~}{\mathrm{\Psi }}\underset{\mathrm{\Lambda },\mathrm{\Lambda }^{}}{}\stackrel{~}{\mathrm{\Psi }}|\stackrel{~}{p}_\mathrm{\Lambda }\stackrel{~}{\varphi }_\mathrm{\Lambda }|B|\stackrel{~}{\varphi }_\mathrm{\Lambda }^{}\stackrel{~}{p}_\mathrm{\Lambda }^{}|\stackrel{~}{\mathrm{\Psi }}$$
(7)
by any operator $`B`$ entirely localized within the augmentation region. Equation (7) has the same range of validity as Eq. (6) does and allows a further acceleration of the convergence by using a well-chosen operator $`B`$ and adding the corresponding zero operator to the expression for the expectation value.

### B The LDA+U total energy functional

For transition-metal oxides, the $`d`$-orbitals are well localized and keep a strong atom-like character. Even though LDA provides a good approximation for the average Coulomb energy of the $`d`$-$`d`$ interactions, it fails to describe correctly the strong Coulomb and exchange interaction between electrons in the same $`d`$-shell. The main intention of LDA+U is to identify these atomic orbitals and to describe their electronic interactions as strongly correlated states. The other orbitals are delocalized and considered to be properly described by the LDA. The procedure is to eliminate the averaged LDA energy contribution of these atom-like orbitals from the LDA total energy functional $`E^{\mathrm{LDA}}`$, and to add an orbital- and spin-dependent correction. The total energy within the LDA+U method then has the form

$`E`$ $`=`$ $`E^{\mathrm{LDA}}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{t,\sigma }{}}{\displaystyle \underset{i,j,k,l}{}}\chi _i^t;\chi _k^t|V_{ee}|\chi _j^t;\chi _l^tn_{i,j}^{t,\sigma }n_{k,l}^{t,\sigma }`$ (8)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{t,\sigma }{}}{\displaystyle \underset{i,j,k,l}{}}\left(\chi _i^t;\chi _k^t|V_{ee}|\chi _j^t;\chi _l^t\chi _i^t;\chi _k^t|V_{ee}|\chi _l^t;\chi _j^t\right)n_{i,j}^{t,\sigma }n_{k,l}^{t,\sigma }`$ (9)
$``$ $`{\displaystyle \underset{t}{}}\left[{\displaystyle \frac{1}{2}}U{\displaystyle \underset{\sigma ,\sigma ^{}}{}}N^{t,\sigma }(N^{t,\sigma ^{}}\delta _{\sigma ,\sigma ^{}}){\displaystyle \frac{1}{2}}J{\displaystyle \underset{\sigma }{}}N^{t,\sigma }(N^{t,\sigma }1)\right],`$ (10)
where $`N^{t,\sigma }=_in_{i,i}^{t,\sigma }`$ is the average occupation of the $`d`$-shell for each spin direction as obtained from the $`d`$-orbital occupancies $`n_{i,j}^{t,\sigma }`$. $`U`$ and $`J`$ are the Coulomb self-energy and the exchange parameter, respectively. The expressions $`\chi _i^t;\chi _k^t|V_{ee}|\chi _j^t;\chi _l^t`$ are the four-center matrix elements of the screened Coulomb interaction $`V_{ee}`$. An additional requirement of the LDA+U approach is that the additional energy is applied only to the valence electrons, which are re-optimized while constrained to remain orthogonal to the core states.

### C Orbital occupations

The rationale behind the definition of the orbital occupations is to project the density matrix onto one consisting of a restricted set of localized orbitals $`|\chi _i^t`$, where each orbital index $`i`$ stands for an angular momentum and $`t`$ for an atomic site. Here we focus in particular on the $`d`$-orbitals of Ni.

The projected density matrix operator has the form

$$\overline{\rho }^\sigma =\underset{i,j}{}|\chi _i^tn_{i,j}^{t,\sigma }\chi _j^t|.$$
(11)
The density matrix elements $`n_{i,j}^t`$ in the localized basis set can be obtained from the condition that the expectation values of a set of projection operators $`P_s`$ are optimally reproduced by the localized basis set, i.e. the weighted square deviation

$`F`$ $`=`$ $`{\displaystyle \underset{s}{}}w_s\left[{\displaystyle \underset{n,𝐤}{}}\mathrm{\Psi }_n^𝐤|P_s^t|\mathrm{\Psi }_n^{𝐤,\sigma }{\displaystyle \underset{i,j}{}}n_{i,j}^{t,\sigma }\chi _i^t|P_s^t|\chi _j^t\right]^2`$
of the projections should be minimized.

Here we choose the projection operators acting on the Ni $`d`$-orbitals in analogy to previous implementations as

$$P_{m,m^{}}^t(𝐫,𝐫^{})=\theta _{\mathrm{\Omega }_t}(𝐫)\delta (|𝐫^{}𝐑_𝐭||𝐫𝐑_𝐭|)Y_{d,m}^t(𝐫)Y_{d,m^{}}^{}(𝐫^{}),$$
(12)
where the site index $`t`$ refers to a particular Ni site, and $`|Y_{d,m}^t`$ is the spherical harmonic centered at site $`t`$. The step functions $`\theta _{\mathrm{\Omega }_t}(r)`$ are unity for $`|rR_t|<r_c^t`$ and zero otherwise. $`Y_{d,m}`$ are the spherical harmonics for $`\mathrm{}=2`$. (Note that $`\chi |P|\chi ^{}=𝑑𝐫𝑑𝐫^{}\chi (𝐫)P(𝐫,𝐫^{})\chi ^{}(𝐫^{})`$.) The radius $`r_c^t`$ for Ni corresponds to an atomic sphere radius, and we have chosen $`r_c^t=2.1`$ a0. The localized orbitals $`|\chi _m^t`$ have been chosen to be identical to those of the spheridized, non-spin-polarized atoms. We used an atom in the $`3d^84s^2`$ configuration of Ni. For such a localized Ni $`d`$-orbital $`|\chi _m^t`$ with magnetic quantum number $`m`$, we obtain

$$I=\chi _m^t|P_{m,m}^t|\chi _m^t=\chi _m^t|\theta _{\mathrm{\Omega }_t}(r)|\chi _m^t,$$
(13)
independent of $`m`$ and $`t`$. Thus we obtain the orbital occupations

$$n_{m,m^{}}^{t,\sigma }=\frac{1}{I}\underset{n,𝐤}{}\mathrm{\Psi }_n^{𝐤,\sigma }|P_{m,m^{}}^t|\mathrm{\Psi }_n^{𝐤,\sigma }$$
(14)
from the minimal condition for $`F`$, which in this case are independent of the weights $`w_s`$.

In the PAW method we obtain the orbital occupations directly from the PS wave functions $`|\stackrel{~}{\mathrm{\Psi }}_n^{𝐤,\sigma }`$ as

$$n_{m,m^{}}^{t,\sigma }=\frac{1}{I}\underset{𝐤,n}{}\stackrel{~}{\mathrm{\Psi }}_n^{𝐤,\sigma }|\stackrel{~}{P}_{m,m^{}}^t|\stackrel{~}{\mathrm{\Psi }}_n^{𝐤,\sigma },$$
(15)
using the pseudo-version

$$\stackrel{~}{P}_{m,m^{}}^t=\underset{\mathrm{\Lambda },\mathrm{\Lambda }^{}}{}|\stackrel{~}{p}_\mathrm{\Lambda }\varphi _\mathrm{\Lambda }|P_{m,m^{}}^t|\varphi _\mathrm{\Lambda }^{}\stackrel{~}{p}_\mathrm{\Lambda }^{}|+\mathrm{\Delta }\stackrel{~}{P}_{m,m^{}}^t$$
(16)
of the projection operator $`P_{m,m^{}}^t`$. Consistent with the PAW formalism, the small correction

$$\mathrm{\Delta }\stackrel{~}{P}_{m,m^{}}^t=P_{m,m^{}}^t\underset{\mathrm{\Lambda },\mathrm{\Lambda }^{}}{}|\stackrel{~}{p}_\mathrm{\Lambda }\stackrel{~}{\varphi }_\mathrm{\Lambda }|P_{m,m^{}}^t|\stackrel{~}{\varphi }_\mathrm{\Lambda }^{}\stackrel{~}{p}_\mathrm{\Lambda }^{}|$$
(17)
is ignored in the present calculations because it can be considered the pseudo-version of the zero operator.

Note that the present choice for the projection operator may not be the most elegant form. Here we suggest another version using projector functions $`q_i^t|`$. The construction is analogous that of the projector functions in the PAW method, but the purpose is different: namely, the goal is to decompose a wave function into local orbitals such that

$$\underset{i,t}{}|\chi _i^tq_i^t|\mathrm{\Psi }=|\mathrm{\Psi }$$
(18)
holds for any wave function $`|\mathrm{\Psi }`$ that can be represented exactly as a superposition of the local orbitals $`|\chi _i^t`$. This requirement is fulfilled whenever $`q_i^t|\chi _j^t^{}=\delta _{i,j}\delta _{t,t^{}}`$. The projector functions can be obtained from a set of functions $`|f_i^t`$, which are linearly independent and localized within an augmentation region as

$$\underset{j,t^{}}{}f_i^t|\chi _j^t^{}q_j^t^{}|=f_i^t|$$
(19)
after multiplication with the inverse of $`f_i^t|\chi _j^t^{}`$. If we now use

$$\overline{P}_{i,j}^t=|q_i^tq_j^t|,$$
(20)
we obtain the matrix elements in the representation of localized orbitals directly as

$$n_{i,j}^t=\underset{k,l,t^{\prime \prime },t^{\prime \prime \prime }}{}n_{k,l}^{t^{\prime \prime },t^{\prime \prime \prime }}\chi _k^{t^{\prime \prime }}|q_i^tq_j^t|\chi _l^{t^{\prime \prime \prime }}=\underset{𝐤,n}{}\mathrm{\Psi }_n^𝐤|\overline{P}_{i,j}^t|\mathrm{\Psi }_n^𝐤.$$
(21)
It is now a simple matter to obtain the corresponding pseudo-projection operator

$$\stackrel{~}{\overline{P}}_{i,j}^t=\underset{\mathrm{\Lambda },\mathrm{\Lambda }^{}}{}|\stackrel{~}{p}_\mathrm{\Lambda }\varphi _\mathrm{\Lambda }|q_i^tq_j^t|\varphi _\mathrm{\Lambda }^{}\stackrel{~}{p}_\mathrm{\Lambda }^{}|.$$
(22)
This expression has a number of advantages: (i) It is the most general expression because any operator can be expressed in a fully separable expression, and, unlike Eq. (12), it is not limited to a semi-local form. (ii) The step function in Eq. (12) introduces a potential that is discontinuous, whereas the operator suggested here can be constructed such that the resulting potential is continuous and differentiable, while still being fully localized in a well-defined region. (iii) The separable form of the operator renders the evaluation numerically convenient. In particular, it is conceivable to evaluate $`\mathrm{\Delta }\stackrel{~}{P}`$ explicitly if necessary, whereas the corresponding operation using the semi-local projection operator would be computationally prohibitive. (iv) We can systematically construct projection operators that specifically act on a well-defined portion of the Hilbert space: By enlarging the set of functions $`|\chi _i`$ to include also orbitals that should not be represented in the restricted set, they are excluded by simply setting the corresponding matrix elements of the density matrix $`N_{i,j}`$ to zero.

### D Coulomb and exchange parameters

The four-center integrals used in the expression of the LDA+U total energy are defined as

$$\chi _i^t;\chi _j^t|V_{ee}|\chi _k^t;\chi _l^t=𝑑𝐫_1𝑑𝐫_2\chi _i^t(𝐫_1)\chi _j^t(𝐫_2)v_{ee}(𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐})\chi _k^t(𝐫_1)\chi _l^t(𝐫_2),$$
(23)
where $`v_{ee}(𝐫,𝐫^{})`$ is the screened Coulomb interaction between two electrons.

If we choose localized $`d`$-orbitals that are described by an atomic $`d`$-wave function $`\chi _m^t(𝐫)=\chi _d(|𝐫𝐑_t|)Y_{d,m}(𝐫𝐑_t)`$ with magnetic quantum number $`m`$ (where $`Y_{d,m}`$ are the spherical harmonics for $`\mathrm{}=2`$) and furthermore assume that the static dielectric function $`ϵ`$ is constant in space, we can exploit the multipole expansion of $`1/|𝐫_1𝐫_2|`$:

$$v_{ee}(𝐫_1,𝐫_2)=\frac{1}{ϵ\left|𝐫_1𝐫_2\right|}=\frac{1}{ϵ}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\frac{4\pi }{2\mathrm{}+1}\frac{r_<^{\mathrm{}}}{r_>^{\mathrm{}+1}}\underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}Y_\mathrm{}m(𝐫_\mathrm{𝟏})Y_\mathrm{}m^{}(𝐫_\mathrm{𝟐}).$$
(24)
Here $`r_<`$ and $`r_>`$ denote the smallest and largest values of $`r_1`$ and $`r_2`$, respectively. Under these assumptions we can transform Eq. (23) into

$$\chi _1^t;\chi _3^t|V_{ee}|\chi _2^t;\chi _4^t=\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\frac{4\pi }{2\mathrm{}+1}\underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}\mathrm{}_1,m_1|Y_{\mathrm{},m}|\mathrm{}_2,m_2\mathrm{}_3,m_3|Y_{\mathrm{},m}^{}|\mathrm{}_4,m_4F^{\mathrm{}},$$
(25)
where $`(\mathrm{}_i,m_i)`$ are the angular momenta quantum numbers of $`|\chi _i^t`$, $`\mathrm{},m|Y_{\mathrm{}^{\prime \prime },m^{\prime \prime }}|\mathrm{},m^{}`$ the Gaunt coefficients, and $`F^{\mathrm{}}`$ the so-called screened Slater’s integrals. Because of the special properties of the Gaunt coefficients, only $`F^0,F^2,`$ and $`F^4`$ contribute to the Coulomb integrals:

$`F^{\mathrm{}}={\displaystyle \frac{1}{ϵ}}{\displaystyle _0^{\mathrm{}}}𝑑r_1{\displaystyle _0^{\mathrm{}}}𝑑r_2r_1^2r_2^2\chi _d^2(r_1)\chi _d^2(r_2){\displaystyle \frac{r_<^{\mathrm{}}}{r_>^{\mathrm{}+1}}}.`$ (26)
The parameters $`U`$ and $`J`$ are identified with averages of the Coulomb and exchange integrals, which are related to the Slater integrals $`F^0,F^2,`$ and $`F^4`$ by the properties of Clebsch-Gordan coefficients, Eq. (25):

$`U`$ $`=`$ $`{\displaystyle \frac{1}{(2\mathrm{}+1)^2}}{\displaystyle \underset{m,m^{}}{}}\chi _m^t;\chi _m^{}^t|V_{ee}|\chi _m^t;\chi _m^{}^t=F^0,`$ (27)
$`J`$ $`=`$ $`{\displaystyle \frac{1}{(2\mathrm{})(2\mathrm{}+1)}}{\displaystyle \underset{mm^{},m^{}}{}}\chi _m^t;\chi _m^{}^t|V_{ee}|\chi _m^{}^t;\chi _m^t=(F^2+F^4)/14,`$ (28)
The dielectric constant and therefore the Coulomb and exchange parameters $`U`$ and $`J`$ are not known a priori. Usually, they are obtained from a constrained DFT calculation. Here, we are interested in how sensitively the results depend on the choice of Coulomb parameters, and which Coulomb parameters will provide the best agreement with reality as probed by optical absorption. Therefore, we adopt the general form for the four-center integrals as function of the $`U`$ and $`J`$ suggested by the arguments provided above, and perform calculations for different $`U`$ values, namely $`U=5`$ eV$`\times I^2`$ and $`U=8`$ eV$`\times I^2`$. In the following we suppose that $`I^2=1`$. Because the results are fairly insensitive to the exchange parameter $`J`$, we have adopted $`J=0.95`$ eV$`\times I^2`$ from previous constrained LDA calculations. The third mandatory relation is obtained from the work of DeGroot et al., who determined the ratio $`F^4/F^2`$ for transition-metal oxides to be between 0.62 and 0.63. We therefore adopt a ratio $`F^4/F^2=0.625`$.

### E Hamiltonian

The pseudo-Hamiltonian operator

$$\stackrel{~}{H}_\sigma =\stackrel{~}{H}_\sigma ^{LDA}+\stackrel{~}{H}_\sigma ^U,$$
(29)
which acts on the PS wave functions, is obtained as the derivative of the total energy functional with respect to the two-center pseudo-density matrix operator $`\stackrel{~}{\rho }=_{n,𝐤}|\stackrel{~}{\mathrm{\Psi }}_n^{𝐤,\sigma }f_n^{𝐤,\sigma }\stackrel{~}{\mathrm{\Psi }}_n^{𝐤,\sigma }|`$. The non-LDA contribution of the LDA+U Hamiltonian is then obtained as the product of the derivative $`V_{m,m^{}}^{t,\sigma }`$ of the non-LDA contribution to the total energy and the projection operator $`\stackrel{~}{P}`$, which is the derivative of the occupation with respect to the two-center pseudo-density matrix operator. Thus we obtain

$$\stackrel{~}{H}_\sigma ^U=\underset{t,m,m^{}}{}\frac{1}{I}\stackrel{~}{P}_{m,m^{}}^tV_{m^{},m}^{t,\sigma },$$
(30)
where

$`V_{m_1,m_2}^{t,\sigma }`$ $`=`$ $`{\displaystyle \underset{m_3,m_4}{}}\chi _{m_1}^t;\chi _{m_3}^t|V_{ee}|\chi _{m_2}^t;\chi _{m_4}^tn_{m_3,m_4}^{t,\sigma }`$ (31)
$`+`$ $`{\displaystyle \underset{m_3,m_4}{}}\left[\chi _{m_1}^t;\chi _{m_3}^t|V_{ee}|\chi _{m_2}^t;\chi _{m_4}^t\chi _{m_1}^t;\chi _{m_3}^t|V_{ee}|\chi _{m_4}^t;\chi _{m_2}^t\right]n_{m_3,m_4}^{t,\sigma }`$ (32)
$``$ $`{\displaystyle \underset{\sigma ^{}}{}}\left[U(N^{t,\sigma ^{}}{\displaystyle \frac{1}{2}}\delta _{\sigma ,\sigma ^{}})\delta _{\sigma ,\sigma ^{}}J(N^{t,\sigma ^{}}{\displaystyle \frac{1}{2}})\right]\delta _{m,m^{}}.`$ (33)
The LDA contribution of the pseudo-Hamiltonian has the usual form:

$`\stackrel{~}{H}_\sigma ^{LDA}={\displaystyle \frac{^2}{2}}+\stackrel{~}{v}+{\displaystyle \underset{\mathrm{\Lambda }_1,\mathrm{\Lambda }_2}{}}|\stackrel{~}{p}_{\mathrm{\Lambda }_1}\left[\varphi _{\mathrm{\Lambda }_1}\left|{\displaystyle \frac{^2}{2}}+v^1\right|\varphi _{\mathrm{\Lambda }_2}\stackrel{~}{\varphi }_{\mathrm{\Lambda }_1}\left|{\displaystyle \frac{^2}{2}}+\stackrel{~}{v^1}\right|\stackrel{~}{\varphi }_{\mathrm{\Lambda }_2}\right]\stackrel{~}{p}_{\mathrm{\Lambda }_2}|.`$ (34)
For the LDA+U calculation, we first performed a self-consistent LSDA calculation using the all-electron PAW method, and then used the self-consistent potential to construct the LSDA Hamiltonian for a large number of $`𝐤`$-points in the Brillouin zone. Next, the Hubbard correction is added to the LSDA Hamiltonian as given by Eq. (29), and the new Hamiltonian iterated until the occupation numbers $`n_{m,m^{}}^{t,\sigma }`$ have converged. We did not use the so-called second-variation procedure for self-consistent LDA+U, in which the LSDA potential is updated, because it failed to yield any improvement on calculations similar to ours.

## III Ground state of NiO

### A LSDA ground state

The ground state of NiO has been calculated using the PAW method, and the density of states (DOS) is calculated from the self-consistent PAW potential using the tetrahedron method for the Brillouin-zone integration. Figure 1 presents the atom-resolved DOS in the augmentation region. As can be seen, LSDA produced an antiferromagnetic insulating ground state with a small band gap. The LSDA magnetic moment is 0.95 $`\mu _B`$, and mainly due to the Ni$`_{e_g}`$ band splitting. This value is much smaller than the experimental value (1.64–1.90 $`\mu _B`$, Refs. ). The calculated band gap of about 0.1 eV is also much smaller than the experimental one (3.0–4.0 eV). The most interesting feature of our LSDA DOS is that the $`d`$-states of Ni dominate the region in the vicinity of the band gap, and that the top of the valence state is of Ni$`_{e_g}`$ type for the first and Ni$`_{t_2g}`$ for the second spin. This electronic structure suggests that the band gap is of Mott-Hubbard type. Hence, this LSDA picture of NiO disagrees completely with experiment. It is surprising that the quasi-particle calculation within the so-called GW approximation, performed by Aryasetiawan and Gunnarsson, produced results qualitatively similar to LSDA except for an increased band gap of 6 eV and an increased magnetic moment of 1.6 $`\mu _B`$. The quantitative change is the reduction of the O $`2p`$ bandwidth by almost 1 eV. However, a recent self-consistent model GW calculation by Massida et al., in which the dynamic effects were neglected, produced other results than those of Aryasetiawan and Gunnarsson. It was argued by Massida al. that the results of the former GW calculation are not quite self-consistent, presumably because of the additional non-local ad-hoc potential that is adjusted to the GW calculation in each self-consistent step. The main difference between the two reported GW model calculations is that the latter calculation produced (i) a spreading of the Ni $`d`$-states over the entire valence bandwidth, (ii) a vanishing gap between the O $`2p`$ and Ni $`3d`$ and, most importantly, (iii) an enhancement of O $`2p`$ states at the valence band maximum. The latter effect attributes the origin of the band gap mainly to a charge transfer gap because this gap is now between the O $`2p`$ and Ni 3$`d`$ conduction states. Next, we will show that this latter finding is in agreement with the results of the LDA+U model.

### B LDA+U ground state

We have used our implementation of the LDA+U model to determine the ground-state electronic structure of NiO. Although it is common practice to use the $`U`$ extracted from a constrained LDA calculation, we adopt a different point of view here. As stated in the introduction and in agreement with recent results reported in the literature, we believe that the value of $`U`$ extracted from constrained LDA is not the best possible choice. Therefore we have determined the electronic structure of NiO for an intermediate $`U`$ of 5 eV as well as for a larger value of 8 eV.

Figure 2 shows our calculated LDA+U DOS for $`U=5`$ and 8 eV. The energy band gap is found to be 2.8 and 4.1 eV, respectively. The total antiferromagnetic spin moment is 1.73 and 1.83 $`\mu _B`$, respectively. Our DOS obtained for $`U=8`$ eV agrees well with previous LDA+U calculations. However, the DOS for $`U`$ = 5 eV is in better agreement with the GW model calculation of Massida et al. The top of valence band is reinforced by the O $`2p`$ states, rendering the band gap a mixture of charge-transfer type and Ni $`dd`$-like excitations. In agreement with previous LDA+U calculations, the spin majority Ni $`e_g`$ states are pushed towards lower energies, and the energy difference between $`e_g^{}`$ and $`e_g^{}`$ is about 11 eV for $`U=8`$ eV and 8.6 eV for $`U=5`$ eV. Here again the GW model yields a value of about 9 eV for this splitting, in good agreement with our results for $`U`$ = 5 eV.

## IV Dielectric function

For insulators, the imaginary part of the macroscopic dielectric function is obtained within the random phase approximation (RPA) in the long wave length limit without local-field effects as

$$ϵ_2(\omega )=\underset{q0}{lim}\underset{v,c}{}\underset{𝐤}{}\frac{8\pi ^2}{\mathrm{\Omega }q^2}\left|M_{v,𝐤𝐪}^{c,𝐤}\right|^2f_{v,𝐤}(1f_{c,𝐤})\delta (E_c^𝐤E_v^𝐤\mathrm{}\omega ).$$
(35)
Here $`M_{v,𝐤𝐪}^{c,𝐤}`$ are the interband transition matrix elements, $`f_{v,𝐤}`$ is the zero-temperature Fermi distribution, $`\mathrm{\Omega }`$ is the cell volume, $`c`$ denotes the conduction-band and $`v`$ the valence-band index. In the case of a local potential, the interband transition matrix elements are given by

$$\underset{q0}{lim}M_{v,𝐤𝐪}^{c,𝐤}=\frac{𝐪}{ϵ_{c,𝐤}ϵ_{v,𝐤}}\mathrm{\Psi }_v^𝐤|𝐩|\mathrm{\Psi }_c^𝐤,$$
(36)
where the matrix elements $`\mathrm{\Psi }_v^𝐤|𝐩|\mathrm{\Psi }_c^𝐤`$ are calculated using the PAW crystal wave functions $`\mathrm{\Psi }^𝐤`$ described by Eq. (5):

$$\mathrm{\Psi }_v^𝐤|𝐩|\mathrm{\Psi }_c^𝐤=\stackrel{~}{\mathrm{\Psi }}_v^𝐤|𝐩|\stackrel{~}{\mathrm{\Psi }}_c^𝐤+\underset{\mathrm{\Lambda },\mathrm{\Lambda }^{}}{}\stackrel{~}{\mathrm{\Psi }}_v^𝐤|\stackrel{~}{p}_\mathrm{\Lambda }\left[\varphi _\mathrm{\Lambda }|𝐩|\varphi _\mathrm{\Lambda }^{}\stackrel{~}{\varphi }_\mathrm{\Lambda }|𝐩|\stackrel{~}{\varphi }_\mathrm{\Lambda }^{}\right]\stackrel{~}{p}_\mathrm{\Lambda }^{}|\stackrel{~}{\mathrm{\Psi }}_c^𝐤,$$
(37)
In the most general case, i.e. where the potential is nonlocal as in LDA+U, a nonlocal contribution has to be added to the interband transition matrix elements. The full derivation is given in the Appendix by Eq. (A7). For NiO the nonlocal contribution to the matrix elements is found to be small, i.e. of a few percent.

Figure 3 shows the imaginary part of the dielectric function calculated within the LSDA. The resulting optical spectrum is not in agreement with experiment (see Fig. 4). In particular, the optical gap is considerably underestimated and the first structure has a much higher intensity compared to experiment. Conversely, for $`U`$ = 5 eV, our calculated imaginary part of the dielectric function within the LDA+U is in a better agreement with experiment, as shown in Fig. 4. The optical band gap and the oscillator strength of the first excitation peak are in excellent agreement with experiment. However, at higher photon energies the agreement with experiment is only qualitative, which is not expected owing to the mean-field approximation of this simple model. A much higher value of $`U`$, i.e. 8 eV, produces a much larger optical gap in contrast to experiment. In agreement with our conclusion that a much smaller value of $`U`$ is required to describe NiO, Dudarev and coworkers also found that $`U=6.2`$ eV reproduces the lattice parameter and the measured electron-energy loss spectra. It is too early to draw a definitive conclusion about the excited states of NiO, as our LDA+U model is a mean-field-like model in which excitonic effects are not included. To our knowledge, no calculations of excitonic effects have so far been attempted. Our results represent the first investigation of the low-lying excited states of NiO that considers all the subtleties of chemical bonding and strong electron-electron interaction.

The interband transitions are responsible for the first structure in the optical spectrum of NiO, located between 4.1 and 4.4 eV. We found that 40.2% of the contribution results from the transition from band 15 (second highest occupied band) to band 17 (lowest empty band), 36.2% from the transition 16$``$18, and 15.9% from the interband transition 16$``$17. To analyze the character of the initial and final states of the interband transitions, the band structure dispersion along some of the high-symmetry directions is shown in Fig. 5 together with the DOS of the states that give rise to the first optical peak. The arrows between the parallel bands indicate the interband transitions from the initial to the final state responsible for the first structure in the optical spectrum. Figure 6 shows the charge density plot of the initial and final states $`\mathrm{\Psi }_n^{𝐤,\sigma }`$ for bands 16 (highest occupied band) and 18 (second lowest empty band) at point $`𝐤=(\frac{127}{120},\frac{127}{360},\frac{\sqrt{3}}{90})\frac{\pi }{a}`$ located between the high-symmetry points K and U, where the optical matrix element value is among the largest. It also shows that the initial state is of mixed O $`2p`$ and Ni $`3d`$ character, whereas the final state mainly is of Ni $`e_g`$ character. The optical transitions then are between the O $`2p`$ and the Ni $`3d`$ states, resulting in excitations of the charge-transfer type. Our interpretation differs from that of Fujimori and Minami, who used a configuration interaction within the metal-ligand cluster and claimed that the $`dd`$ charge-transfer transitions are the origin of fundamental edge. The drawback of the cluster calculation is that in reality the oxygen $`2p`$ orbitals are delocalized and therefore not well described in a small cluster. On the other hand, an earlier, band-structure-based interpretation by Messick and coworkers assigned the peak to one electron interband transitions associated with Ni $`3d`$ to the Ni $`4s`$ states. This interpretation is not correct either because the Ni $`4s`$ are far above the top of the valence states (greater than 6 eV), and only quadrupolar interband transitions are permitted between the $`3d`$ and $`4s`$, which substantially reduces the peak intensity.

## V Conclusion

We have presented a new implementation of LDA+U model based on the PAW method, which is an all-electron method without any shape approximation to the potential or the charge density. We tested the method on NiO and obtained results that are in good agreement with previous LDA+U calculations and a recent GW model calculation. In particular, we obtained the correct antiferromagnetic insulating ground state of NiO.

We discussed the results in terms of the strength of the Hubbard interaction $`U`$. The optimum value of $`U`$ has been determined by comparison with the experimental dielectric function as well as with the ground state properties. We observed a large enhancement of the O $`2p`$ character at the top of the valence state, resulting in a more charge-transfer than Ni $`dd`$ LSDA -type band gap. The calculated antiferromagnetic moment is in good agreement with experiment.

Our calculated dielectric function for an intermediate value of $`U`$, namely 5 eV, is in good agreement with experiment. The low-lying, strong structure in the optical spectrum has been assigned to an interband transition from O $`2p`$ states at top of the valence band to the Ni $`e_g`$ states at the conduction-band bottom. Hence the origin of the first optical peak is due to a charge-transfer excitation.

Our calculation is supported by a recent LDA+U calculation by Dudarev and coworkers, who also argue that a much smaller value of $`U`$ than the one obtained from constrained LDA calculation is needed to describe the electron energy loss spectra and the equilibrium lattice parameter. It should be interesting to apply this method to other transition-metal oxides and check the applicability of LSDA+U for producing excitation energies.

## VI Acknowledgment

Part of this work was done during our (O.B. and M.A.) visit to the Ohio State University in the summer of 1998. We thank J. W. Wilkins and J. G. LePage for useful discussions. Supercomputer time was provided by CINES (project gem1101) on the IBM SP2 and by the Université Louis Pasteur on the SGI O2000 supercomputer.

## A Optical transition matrix elements in LDA+U

The interband transition matrix elements for a given Hamiltonian are obtained as follows:

$`M_{v,𝐤𝐪}^{c,𝐤}`$ $`=`$ $`\mathrm{\Psi }_{v,𝐤𝐪}|e^{i\mathrm{𝐪𝐫}}|\mathrm{\Psi }_{c,𝐤}={\displaystyle \frac{\mathrm{\Psi }_{v,𝐤𝐪}|(ϵ_{v,𝐤𝐪}ϵ_{c,𝐤})e^{i\mathrm{𝐪𝐫}}|\mathrm{\Psi }_{c,𝐤}}{ϵ_{v,𝐤𝐪}ϵ_{c,𝐤}}}`$ (A1)
$`=`$ $`{\displaystyle \frac{\mathrm{\Psi }_{v,𝐤𝐪}|He^{i\mathrm{𝐪𝐫}}e^{i\mathrm{𝐪𝐫}}H|\mathrm{\Psi }_{c,𝐤}}{ϵ_{v,𝐤𝐪}ϵ_{c,𝐤}}}={\displaystyle \frac{\mathrm{\Psi }_{v,𝐤𝐪}|[H,e^{i\mathrm{𝐪𝐫}}]_{}|\mathrm{\Psi }_{c,𝐤}}{ϵ_{v,𝐤𝐪}ϵ_{c,𝐤}}}.`$ (A2)
The commutator involving of the LDA contribution $`H^{\mathrm{LDA}}`$ to the Hamiltonian $`H=H^{\mathrm{LDA}}+H^U`$ is

$$[H^{\mathrm{LDA}},e^{i\mathrm{𝐪𝐫}}]_{}=\frac{1}{2}[^2,e^{i\mathrm{𝐪𝐫}}]_{}=\frac{1}{2}(2i𝐪+q^2)=\mathrm{𝐪𝐩}+O(q)^2.$$
(A3)
The quadratic and higher-order terms in $`𝐪`$ can be ignored in the long-wavelength limit appropriate for optical transitions. The commutator involving the non-LDA part $`H^U`$ is obtained as

$$[H^U,e^{i\mathrm{𝐪𝐫}}]_{}=i𝐪[H^U,𝐫]_{}+O(q)^2=i𝐪\underset{t,m,m^{}}{}\frac{V_{m,m^{},\sigma }^t}{I}[P_{m^{},m}^t,𝐫]_{}+O(q)^2$$
(A4)
Next we use the relation which holds for the special form of the projector operator presented in Eq. (12)

$$P_{m^{\prime \prime },m^{}}^t^{}|\varphi _{t,\mathrm{},m,\alpha }=\theta _{\mathrm{\Omega }_t}(𝐫)|\varphi _{t,\mathrm{},m^{\prime \prime },\alpha }\delta _{t,t^{}}\delta _{\mathrm{},2}\delta _{m,m^{}},$$
(A5)
i.e. the projection operator changes the magnetic quantum number of a specific $`d`$-like partial wave from $`m`$ to $`m^{\prime \prime }`$ and truncates it beyond the atomic sphere $`\mathrm{\Omega }_t`$. Hence we obtain

$$\varphi _\mathrm{\Lambda }|[P_{m^{},m}^t,𝐫]_{}|\varphi _\mathrm{\Lambda }^{}^t=P_{m^{},m}^t\varphi _\mathrm{\Lambda }|\theta _{\mathrm{\Omega }_t}(𝐫)𝐫|\varphi _\mathrm{\Lambda }^{}\varphi _\mathrm{\Lambda }|\theta _{\mathrm{\Omega }_t}(𝐫)𝐫|P_{m,m^{}}^t\varphi _\mathrm{\Lambda }^{}.$$
(A6)
Finally, we obtain the expression of the matrix elements for the dipole transition, with the PAW LDA+U formalism:

$`\underset{q0}{lim}M_{v,𝐤𝐪}^{c,𝐤}`$ $`=`$ $`{\displaystyle \frac{𝐪}{(ϵ_{v,𝐤𝐪}ϵ_{c,𝐤})}}\{\stackrel{~}{\mathrm{\Psi }}_{v,𝐤}|𝐩|\stackrel{~}{\mathrm{\Psi }}_{c,𝐤}+{\displaystyle \underset{\mathrm{\Lambda },\mathrm{\Lambda }^{}}{}}\delta _{t,t^{}}\stackrel{~}{\mathrm{\Psi }}_{v,𝐤}|\stackrel{~}{p}_\mathrm{\Lambda }[\varphi _\mathrm{\Lambda }|𝐩|\varphi _\mathrm{\Lambda }^{}\stackrel{~}{\varphi }_\mathrm{\Lambda }|𝐩|\stackrel{~}{\varphi }_\mathrm{\Lambda }^{}`$ (A7)
$`+`$ $`i{\displaystyle \underset{m,m^{}}{}}{\displaystyle \frac{V_{m,m^{}}}{I}}(P_{m,m^{}}^t\varphi _\mathrm{\Lambda }|\theta _{\mathrm{\Omega }_t}(𝐫)𝐫|\varphi _\mathrm{\Lambda }^{}\varphi _\mathrm{\Lambda }|\theta _{\mathrm{\Omega }_t}(𝐫)𝐫|P_{m,m^{}}^t\varphi _\mathrm{\Lambda }^{})]\stackrel{~}{p}_\mathrm{\Lambda }^{}|\stackrel{~}{\mathrm{\Psi }}_{c,𝐤}\}.`$ (A8)
The difference between wave functions $`|\stackrel{~}{\mathrm{\Psi }}_{v,𝐤}`$ and $`|\stackrel{~}{\mathrm{\Psi }}_{v,𝐤𝐪}`$ has been ignored because it only contributes to terms that are proportional to $`q^2`$, which are ignored in the long-wavelength limit.


text:
# Key aspects in the photometric monitoring of 𝛿 Sct stars

## 1 Introduction

In the last two decades a great observational effort has been undertaken to solve the light curves of $`\delta `$ Sct stars. These complicated variables, located on or just above the Zero–Age–Main–Sequence, constitute the lower limit of the classical instability strip. The presence of convective zones makes these stars very interesting from the evolutionary point of view, since various stages, from pre–main sequences to giant, can be observed.

Photometric monitoring is the most practiced approach to study their properties and several stars have been deeply investigated; for a review, see Poretti (2000). However, the spectroscopic approach tells us that many modes that are not photometrically detectable are actually excited (Mantegazza 2000).

The observed frequencies are between 60 and 405 $`\mu `$Hz and multiperiodicity is very common; rotation acts as an important factor in the increase of the number of excited modes. Only a few stars show a monoperiodic behaviour above the current limit of the detectable amplitude from ground, i.e. $``$1 mmag. Of course, the number of detected modes is strongly dependent on the accuracy, the number and the sampling of the measurements.

The observed modes are in the domain of pressure ($`p`$) modes; there is no observational evidence that gravity ($`g`$) modes are excited in $`\delta `$ Sct stars, even if some cases suggest they might be. At the moment, $`g`$-modes seem to be present only in the $`\gamma `$ Dor stars, which in turn do not show $`p`$–modes.

How the photometric results can be used to identify modes (i.e. to classify the oscillation in terms of quantum numbers $`n`$, $`\mathrm{}`$ and $`m`$) is discussed by Garrido (2000).

## 2 The state–of–art in the study of $`\delta `$ Sct stars

Several approaches have been used in the past to attack the problem of the detection of excited modes in their light curves. The debate on the nature of the variations, multiperiodic or irregular, occupied the years from 1974 to 1980. At that time the observations consisted in time–series collected over a few nights, a practice that had been borrowed from the study of high amplitude $`\delta `$ Sct stars. In the years that followed, a number of teams observed more carefully $`\delta `$ Sct stars and clearly evidenced their multiperiodicity, detecting stable frequencies and obtaining very satisfactory fits. It appeared also evident that nonradial modes are excited. The number of frequencies is so high that the aliases at 1 cd-1 become a serious problem in the mode detection. For that, multisite campaigns were planned and the best results now available have been obtained from these ground–based observations.

The effort made by the observers these last years allowed us to detect a large number of frequencies in many stars: we will recall here the cases of FG Vir, XX Pyx, 4 CVn and BH Psc. By means of multisite campaigns we are able to detect terms with amplitude less than 1 mmag. The mode identification techniques are based on both the phase shifts and amplitude ratios of the light and colour curves and the synergic approach performed by considering spectroscopic curves. Their full exploitation should guarantee an important role for our researches in stellar physics.

Figure 1 shows the frequency content of three well-studied pulsators in which all the terms with a half-amplitude larger than 1 mmag are detected. It appears that their content is completely different: high frequencies only for XX Pyx, low frequencies only for 4 CVn, two groups of intermediate frequency values for FG Vir. The combination terms are also reported in the upper part of each panel; 4 CVn displays both $`f_if_j`$ and $`f_i+f_j`$ terms, located right and left of the bunch of independent frequencies. FG Vir displays only three combination modes in the same region of the independent frequencies. XX Pyx displays only a 2$`f`$ harmonic, at 76.2 cd-1.

The obvious conclusion is that $`\delta `$ Sct stars are very complicated pulsators and that a general recipe cannot be applied to predict the range of the excited modes or the importance of combination terms. However, a curious similarity in the frequency content was noted by Mantegazza & Poretti (1999) comparing 4 CVn to HD 2724: at least 7 frequencies have almost the same value.

## 3 The amplitude variations

The variations of the amplitudes have been detected when a $`\delta `$ Sct star has been re–observed in order to investigate some unclear facts: a comparison of the mode amplitudes in the two different seasons often evidences strong differences. Even if the data sampling can affect the reliability of these conclusions, some results seem to be very much in favour of a physical variability (Breger 2000). These variations can be seen as a further complication, even if we can argue that the disappearance or the damping of some terms can enhance or make discernible other modes, increasing the number of known frequencies.

Amplitude variations are observed in a large variety of stars, both multiperiodic (XX Pyx, Handler et al. 1997; 4 CVn, Breger 2000; 44 Tau, Poretti et al. 1992) and monoperiodic (28 And, Rodríguez et al. 1998; BF Phe, Poretti et al. 1996). The observations of mode growth in the light curve of a relatively simple pulsator as V663 Cas can clarify what happens in much more complicated cases: some modes can be damped and then re–excited (Poretti et al. 1996). Some observational facts favour this explanation even if the model of a beating between two close frequencies with similar, constant amplitude cannot be ruled out.

BH Psc is another example. It was observed by our group in 1991 (on two occasions), 1994 and 1995 at the European Southern Observatory, Chile. The analysis of the first two campaigns showed a very complex light variability resulting from the superimposition of more than 10 pulsation modes with frequencies between 5 and 12 cd-1 and semi–amplitudes between 17 and 3 mmag (Mantegazza, Poretti, & Zerbi 1995). The fit left a high r.m.s. residual 2.3 times greater than the one measured between the two comparison stars. A second photometric campaign was carried out in October and November 1994; we hoped to reveal more terms and to check the stability of their amplitudes. The analysis of the new data allowed us to single out 13 frequencies (Mantegazza, Poretti, & Bossi 1996). They are concentrated in the region 5–11.5 cd-1: this distribution is slightly larger than the one observed for 4 CVn, but they can be considered very similar. However, more terms should be present, with an amplitude below 3 mmag, since about 15% of the variance could not be explained with the detected terms and the noise. Moreover, the standard deviation around the mean value was considerably higher in 1991 than in 1994 (26 mmag against 18 mmag); considerations about the amplitudes of the modes confirmed that the pulsation energy was lower in 1994 than in 1991. Moreover, comparing the 1991 September and 1991 October data it is quite evident that we can detect such variations on a baseline of 1 month or less (see Tab. 1, where the results of a preliminary analysis are shown).

## 4 Scientific topics to be matched by MONS

The precision in the measurements achievable by MONS will of course allow us to detect terms with amplitudes much smaller than 1 mmag. This in turn will generate a much more detailed knowledge of the frequency content of $`\delta `$ Sct variables. We could then obviously match one of the requests of the asteroseismology theory, i.e. increasing the number of detected terms. This is quite obviuos. However, the values of the amplitudes of these modes are also an important parameter that can help in the theoretical models. The observational evidence is that these amplitudes are variable. Can MONS then give us an important clue as regard this aspect ?

To reach the goal, observing twice the same field can be very useful. Doing so, we could:

1. gain in frequency resolution since, as shown in Fig. 1, there are many close terms. It should also be noted that multisite campaigns span more than 50 d and that a single 30-50 d MONS run will not have a sufficient frequency resolution to solve already known close peaks. Owing to the expected large number of peaks, interference between different terms can seriously affect the quality of the results: an unappropriate frequency resolution could be considered a weakness of the MONS results;
2. detect amplitude variations in a very clear way. The absence of alias terms combined with the appropriate frequency resolution will allow a clarification on the nature of the amplitude variations, i.e. beating or intrinsic variations.

## 5 A plausible solution

In the current configuration, can MONS satisfy these two requirements ?

Among the $`\delta `$ Sct stars, FG Vir is by far the most studied: not only an excellent observational basis is available, but also more and more thorough theoretical models are proposed on the basis of the observed frequencies (Breger et al. 1999). FG Vir ($`\alpha `$ 12$`{}_{}{}^{\mathrm{h}}14_{}^{\mathrm{m}}15^\mathrm{s}`$, $`\delta 5^{}43^{}00^{\prime \prime }`$, 2000.0) is located very close to $`\beta `$ Vir ($`\alpha `$ 11$`{}_{}{}^{\mathrm{h}}50_{}^{\mathrm{m}}42^\mathrm{s}`$, $`\delta +1^{}45^{}53^{\prime \prime }`$, 2000.0), a F9 V star candidate for searching solar–like oscillations. Then, FG Vir ($`V`$=6.1) can be measured not only as a primary target with the Main Camera, but also with the Forward–pointing star tracker (BST or SI-2) when $`\beta `$ Vir is measured in the Main Camera. From this lucky coincidence, we can obtain two runs on a well–studied $`\delta `$ Sct star, greatly increasing the number of detected frequencies and saving frequency resolution. Moreover, in both cases, it will be possible to measure II Vir, another $`\delta `$ Sct star close to FG Vir.

Is that the only result ? Certainly not. There would be one more bonus. Looking through the backward–pointing star tracker (FST or SI-1) we could see, in both cases, an interesting field, located at 176 from FG and $`\beta `$ Vir. Three intriguing objects are located in this field: namely, BH Psc ($`\alpha `$ 23$`{}_{}{}^{\mathrm{h}}59_{}^{\mathrm{m}}31^\mathrm{s}`$, $`\delta 2^{}50^{}37^{\prime \prime }`$, 2000.0), i.e. a $`\delta `$ Sct showing strong amplitude variations and HD 224639 and HD 224945, two $`\gamma `$ Dor variables. The requirements about frequency resolution and amplitude variations also apply to $`\gamma `$ Dor variables.

As a conclusion, it looks possible for MONS to perform an observational study of $`\delta `$ Sct stars in a very powerful way (increase in the number of detected terms, two–colour photometry, sufficiently long baseline allowing frequency resolution and amplitude monitoring) just considering some scheduling precautions and putting an adequate target in the Main Camera.