text:
# DUKE-TH-02-219 Unexpected Results in the Chiral Limit with Staggered Fermions

## 1 Motivation

Consider a two dimensional $`U(1)`$ lattice gauge theory with staggered fermions. This model has been studied extensively on the lattice for over two decades . In the limit when the lattice gauge coupling goes to zero, the model is believed to describe the continuum two flavor Schwinger model; the fermions are confined and the low energy physical particles are mesons. Further, in the chiral limit three massless and one massive pseudo-scalar mesons emerge like in QCD . However, since in two dimensions continuous chiral symmetries remain unbroken, the chiral condensate vanishes in the chiral limit. The actual prediction is

$$\overline{\psi }\psi m^{1/3},$$
(1)
as discussed in . This means that the chiral susceptibility $`\chi =\frac{}{m}\overline{\psi }\psi `$ diverges in the chiral limit reflecting the presence of massless excitations. This is consistent with the fact that $`M_\pi m^{2/3}`$, i.e., the pion mass vanishes at $`m=0`$. If one computes the susceptibility at $`m=0`$ in a finite box of size $`L\times L`$, one expects

$$\chi L,$$
(2)
i.e., it diverges linearly with $`L`$ .

When the lattice gauge coupling is non-zero, lattice artifacts in the staggered fermion formulation break the chiral symmetry of the two flavor Schwinger model explicitly down to a $`U(1)`$ subgroup. As far as we know, no one has completely analyzed the effects of lattice artifacts on the particle spectrum. Do massless particles remain in the spectrum? The conventional wisdom from higher dimensions is that there should be one massless pion, since the remnant $`U(1)`$ chiral symmetry is expected to break spontaneously at strong couplings. Of course the Mermin-Wagner theorem forbids the breaking of a continuous symmetry in two dimensions, although non-interacting massless bosons can emerge. On the other hand a $`U(1)`$ (or equivalently an $`O(2)`$) symmetry is special. In such a case bosons can interact through topological excitations as discovered by Kosterlitz and Thouless . Thus, the lattice model could either be in a massive or a massless phase depending on the effective couplings of the low energy effective theory describing the bosonic excitations. In the massless phase there are predictions for the behavior of the chiral condensate and the susceptibility based on universality . One expects

$$\overline{\psi }\psi m^{\eta /(4\eta )}.$$
(3)
where $`0\eta 0.25`$. Again the susceptibility $`\chi `$ diverges in the chiral limit reflecting the presence of massless excitations. In this case the finite size scaling formula for the susceptibility is given by

$$\chi L^{2\eta }.$$
(4)
It is amusing that if we set $`\eta =1`$ in eqs. (3) and (4) we recover eqs. (1) and (2).

Most results from earlier work on the lattice Schwinger model favor the view point that there is one massless pion at finite lattice spacings. This is based on the observation that the pion mass decreases with the fermion mass like in the continuum Schwinger model. However, on closer examination the lattice model shows deviations, which become larger at stronger couplings as one might expect . On the other hand, no one has ever found scaling relations expected from universality in a massless phase of an $`O(2)`$ model in two dimensions. In particular no one has been able to confirm relations similar to eqs. (3) or (4). Thus, inspite of the large amount of literature on the subject, questions related to the chiral limit of the lattice Schwinger model with staggered fermions still remain unanswered.

The essential difficulty is the absence of a reliable numerical approach to study interacting fermionic field theories close to the chiral limit. For this reason most previous studies alluded to above, have focused on calculations away from the chiral limit and have used extrapolation techniques to predict the chiral limit. As we will see this can be misleading. In the last few years fermion algorithms based on cluster updates have emerged, which do not suffer from problems that conventional algorithms face near the chiral limit . One can now work directly in the chiral limit using the new approach, a luxury not available with earlier methods. Recently, a fermion cluster algorithm has confirmed the scaling predictions near a Kosterlitz-Thouless transition in a Hubbard type model with unmatched precision . Interestingly, these new ideas can also be applied to study the chiral limit of the lattice Schwinger model at strong couplings . The algorithm is based on a novel loop representation of the model. Although, it may be possible to extend the method to weaker couplings, in this article we concentrate on the strong coupling limit. We focus on the question of whether there are massless excitations at strong couplings by looking for a divergence in the chiral susceptibility as predicted by eq.(4).

## 2 The Model

The two dimensional $`U(1)`$ lattice gauge theory with staggered fermions is described by the action

$`S`$ $`=`$ $`{\displaystyle \frac{1}{g^2}}{\displaystyle \underset{x,\mu ,\nu ,\mu \nu }{}}\mathrm{Re}[U_\mu (x)U_\nu (x+\widehat{\mu })U_\mu ^{}(x+\widehat{\nu })U_\nu ^{}(x)]`$ (5)
$`+{\displaystyle \underset{x,\mu }{}}\overline{\psi }_x\eta _\mu (x)\left(U_\mu (x)\psi _{x+\widehat{\mu }}U_\mu ^{}(x\widehat{\mu })\psi _{x\widehat{\mu }}\right)`$
$`+m{\displaystyle \underset{x}{}}\overline{\psi }_x\psi _x,`$
where $`x(x_1,x_2)`$ represents a lattice site on an $`L\times L`$ lattice, $`\mu ,\nu =1,2`$ represent the two directions and $`\widehat{\mu },\widehat{\nu }`$ are the unit vector in the $`\mu `$ and the $`\nu `$ direction respectively. The phase factors $`\eta _1(x)=1`$ and $`\eta _2(x)=(1)^{x_1}`$ are the staggered fermion phase factors. The gauge fields $`U_\mu (x)`$ are elements of $`U(1)`$ phase factors and the fermion fields $`\psi (x)`$ and $`\overline{\psi }(x)`$ are Grassmann numbers.

In the strong coupling limit ($`g=\mathrm{}`$) it is possible to integrate over the gauge fields and obtain

$`Z(m)`$ $`=`$ $`{\displaystyle [d\psi ][d\overline{\psi }]\underset{x}{}\left(1m\overline{\psi }_x\psi _x\right)}`$ (6)
$`{\displaystyle \underset{x,\mu }{}}\left(1+\overline{\psi }_x\psi _x\overline{\psi }_{x+\mu }\psi _{x+\mu }\right),`$
which can be written as a sum over weights of configurations of monomers and dimers . Each configuration consists of $`n_x=0,1`$ monomers on site $`x`$ and $`b_{x,\mu }=0,1`$ dimers on the bond connecting $`x`$ and $`x+\widehat{\mu }`$. In order for the Grassmann integration to give a non-zero result we need

$$n_x+b_{x,1}+b_{x,2}+b_{x\widehat{1},1}+b_{x\widehat{2},2}=\mathrm{\hspace{0.33em}1}.$$
(7)
Assuming this constraint the partition function can be written as

$$Z(m)=\underset{[n,b]}{}\mathrm{exp}\left[\mathrm{log}(m)\underset{x}{}n_x\right]$$
(8)
A typical configuration is shown in figure 1.

When $`m=0`$ no monomers are allowed and the partition function is given by the number of closely-packed-dimer (CPD) configurations on the lattice. Such configurations are interesting even in condensed matter physics and have played an important role in the study of the 2-d Ising model . The chiral symmetry of staggered fermions is manifest at finite volumes through the fact that the chiral condensate vanishes since it is impossible to find a CPD configuration with just one monomer. The chiral susceptibility on the other hand is non-zero and is a useful observable in the chiral limit. It is equal to the total number of CPD configurations with two monomers, where one of the monomers is constrained to be at a fixed position, divided by the partition function.

## 3 Loop Representation

It is easy to construct a local Metropolis algorithm for the monomer-dimer model when $`m0`$. The algorithm is based on an update which either breaks a dimer into two monomers or vice-versa. This algorithm works reasonably well for $`m0.01`$. However, the algorithm fails in the chiral limit since no monomers are allowed when $`m=0`$. In fact in the chiral limit it is easy to find configurations where simple local updates do not lead to another allowed configuration. Perhaps for this reason, as far as we know, no one has successfully studied the chiral limit.

To construct an algorithm that is applicable in the chiral limit we extend CPD configurations to configurations of loops made up of bonds which include the original or “filled” dimers (represented here by “solid” bonds) and “empty” dimers (represented by “dashed” bonds) such that the partition function can be expressed as a sum over weights of new loop configurations. Figure 2 shows the rules of one such extension. If we ignore monomers each shaded plaquette of the configuration of Fig. 1 carries one of the seven plaquette configurations given on the left side of equations in Fig. 2. On the other hand, the right side of the equations represent the new configurations and their weights. Figure 3 gives an example of a loop configuration.

In the absence of monomers it is easy to check that all constraints are satisfied if each loop is made up of a repeating sequence of filled and empty dimers. This means that there are two allowed configurations associated with each loop. The usefulness of the loop variable is that a dimer system can be updated by “flipping” a loop where filled dimers are emptied and vice versa. The Metropolis acceptance of such a loop flip is found to be reasonable although large loops are not flipped as often as small loops. When monomers are allowed each loop consists of an even number of them. Further, when monomers are present in a loop then it contains a unique pattern of filled and empty dimers and a flip is not allowed. However, close to the chiral limit loops containing monomers are negligible.

The algorithm to measure the susceptibility in the chiral limit is straight forward in the loop representation. Typically, we choose a point at random on the lattice and traverse the loop it is on, starting in the direction of a filled bond. As one goes around the loop, the bonds are flipped and the change in the weight of the configuration is noted. Interestingly, every time a filled bond is emptied one gets a configuration that contributes to the susceptibility. Since its weight relative to the original configuration is known at the time of the flip, it is recorded as a part of the measurement. When the whole loop is flipped one knows the change in the weight of the configuration and a Metropolis accept-reject step can be performed. If the configuration is accepted then one gets a new configuration. Otherwise the loop is flipped back to the original configuration. In any case, the sum over all the weight changes recorded while flipping the filled bonds in the loop divided by one plus the total weight change due to the loop flip is taken as a contribution to the susceptibility during that update. It is possible to show that this algorithm is ergodic.

## 4 Results

Since the remnant $`U(1)`$ chiral symmetry of staggered fermions can break spontaneously in higher dimensions, the chiral condensate can approach a constant in the chiral limit. However, in two dimensions a continuous symmetry cannot break and it is almost guaranteed that

$$\overline{\psi }\psi m^\delta .$$
(9)
with $`0<\delta 1`$. Comparing with eqs. (1) and (3) we see that $`\delta =1/3`$ in the continuum two flavor Schwinger model and $`\delta 1/15`$ is expected from $`O(2)`$ universality in two dimensions. It is important to remember that the $`U(1)`$ chiral symmetry of staggered fermions is non-anomalous. Thus, a non-zero expectation value of $`\overline{\psi }\psi `$ in the chiral limit would definitely imply a spontaneous breaking of the $`U(1)`$ symmetry.

Using the local Metropolis algorithm we calculated the chiral condensate for five masses in the region $`0.01m0.1`$ for $`L=16,24,32`$ and $`64`$. Below $`m=0.01`$ the algorithm slows down and is not very efficient. Our results are shown in figure 4. The solid lines are linear fits of the data to the finite size scaling formula $`\overline{\psi }\psi =A+B/L^2`$ at a fixed value of $`m`$. The constant $`A`$ then yields the condensate at infinite volume for a fixed mass. Figure 5 shows the condensate as a function of the mass for $`L=16`$ and for infinite $`L`$ obtained from fits shown in figure 4. At $`L=16`$ we see that the condensate goes linearly to zero as expected due to finite volume effects111 In a finite box the partition function $`Z(m)`$ is a polynomial of the variable $`m^2`$. The condensate is proportional to the first derivative of $`Z(m)`$ with respect to $`m`$ and hence vanishes linearly with $`m`$.. On the other hand the infinite volume results fit beautifully to a power law of the form $`\overline{\psi }\psi =0.455(6)m^{0.108(5)}`$. Since we are in the strong coupling limit it is not surprising that the power does not match the continuum two flavor result of eq. (1). However, it neither matches the predictions of the two dimensional $`O(2)`$ model based on universality (see eq.(3)).

Before attempting to understand the unexpected value of the power it is useful to confirm that the power law behavior is valid all the way to the chiral limit. A power law for the chiral condensate implies that the chiral susceptibility will diverge at $`m=0`$ in the thermodynamic limit. This divergence typically manifests itself in a finite size scaling of the form

$$\chi =\frac{}{m}\overline{\psi }\psi L^\gamma ,$$
(10)
consistent with eqs. (2) and (4). This behavior is difficult to study with conventional algorithms. However, the loop cluster algorithm discussed in the earlier section is designed exactly for this purpose. In figure 6 we plot our results for the chiral susceptibility obtained with the new algorithm at $`m=0`$ as a function of the lattice size.

Notice that for $`L=16`$ the value of the susceptibility is consistent with $`15`$, the slope of the condensate at $`m=0`$ shown in figure 5.

As a function of $`L`$ the susceptibility behaves differently in different regions. The data is described well by a power law in the region $`8L20`$ with $`\gamma =1.56`$, while a straight line fits the data in the region $`20L40`$. For $`L>40`$ the susceptibility appears to be approaching a constant. If we compare the power law behavior of the susceptibility in the region $`8L20`$ with eq. (4) we find that $`\eta 0.44`$ which is inconsistent with the predictions of $`O(2)`$ universality in which one would expect $`0\eta 0.25`$. However, it explains the strange power law behavior of the infinite volume chiral condensate we encountered in figure 5. If we substitute $`\eta =0.44`$ in eq. (3) we find that $`\overline{\psi }\psi m^{0.12}`$, which is quite close to the $`L=\mathrm{}`$ results of figure 5 in the region $`0.01m0.1`$. Since the power law behavior of the susceptibility does not extend to larger volumes, it is clear that the power law behavior of the infinite volume chiral condensate shown in figure 5 cannot hold all the way to the chiral limit. In fact the susceptibility at $`m=0`$ saturates at large volumes. Equivalently, the first derivative of the chiral condensate with respect to the mass at $`m=0`$ reaches a constant as the volume becomes large222Although we cannot rule out a mild divergence of the susceptibility in figure 6, we think it is unlikely.. This means the infinite volume chiral condensate will also vanish linearly at $`m=0`$. Typically, divergence of a susceptibility such as the one defined in eq. (10), is related to the presence of massless particles in the theory. Our results indicate that there are no such particles.

## 5 Discussion

The results that we have obtained are somewhat surprising. The conventional wisdom is that strong couplings break chiral symmetry at least in three or more dimensions. This suggests that the low energy effective model that describes the $`U(1)`$ chiral dynamics of the lattice Schwinger model with staggered fermions is most likely in the low temperature (massless) phase at strong couplings. This reasoning appears to be in contradiction with our findings since we do not find the expected divergence described by eq. (4) which in turn implies that there are no massless particles at strong couplings.

Is it possible we have missed something? Of course, we rely heavily on our algorithm. It is rather new and may have weaknesses like long auto-correlation times that we have not yet appreciated. In order to alleviate such fears we have compared our Monte Carlo results with exact calculations of the susceptibility at $`m=0`$ on small lattices. Table 1 shows this comparison.

Clearly, we can reproduce the exact results with great precision at least on small lattices. We have not seen any pathologies in the simulations at larger volumes except for the fact that the fluctuations increase which require us to increase the statistics proportionally. Based on this we are prejudiced to believe in our estimate of the errors within a factor of two or three.

Assuming our results are correct, we find that at strong couplings the lattice model is in the high temperature (massive) phase of the low energy effective model. We know from experience in higher dimensions that chiral symmetry breaking effects become weaker at smaller couplings. This, coupled with the fact that there is no deconfinement transition in two dimensions, suggests that it is unlikely that there is a phase transition to a massless phase at weaker couplings. Our results then imply that there are no massless excitations for all $`g0`$ (or equivalently at finite lattice spacings); the non-singlet pseudo-scalar meson (which would be expected to become the Goldstone boson in higher dimensions ) is actually massive at all non-zero couplings. This conclusion may come as a surprise to some, but it does not contradict any known physics as far as we can tell. In particular it does not contradict our expectations that in the zero coupling limit one must recover the continuum two flavor Schwinger model. It is well known that two other non-singlet pseudo-scalar mesons acquire a mass due to lattice artifacts in the staggered fermion formulation. Our data is suggesting that all mesons become massive at non-zero couplings contrary to expectations from higher dimensions where the Goldstone boson will remain massless at finite lattice spacings. In the lattice Schwinger model this boson is still perhaps the lightest. The gauge coupling controls the lattice spacing in the model; typically one uses the mass of the iso-singlet meson to set the lattice spacing. This means that in order to recover the correct continuum limit the singlet and the non-singlet meson masses must become smaller in lattice units as one approaches weaker couplings but their ratio must diverge.

The current work can be extended in many directions. There are results in the large $`N_c`$ and large $`d`$ limits at strong couplings which show that chiral symmetry is indeed broken . Here we have shown that at $`N_c=1`$ and $`d=2`$ chiral symmetry remains unbroken and the theory is in the massive phase. It would be interesting to find the value of $`N_c`$ for $`d=2`$ at which the theory will move into a phase with massless excitations. There is a lot of evidence from lattice simulations and chiral extrapolations that shows lattice QCD with staggered fermions to be in the chirally broken phase. However, as we have seen in the present context, chiral extrapolations can sometimes be misleading. Fortunately, cluster algorithms of the type used here can be developed for any $`N_c`$ and $`d`$ so that these questions can be answered directly in the chiral limit. One can start with the results presented in . It would be useful to confirm that chiral symmetry is indeed broken at $`N_c=3`$ and $`d=4`$. On the algorithmic side, it is exciting that we have a new method to explore the chiral limit in certain gauge theories. Loop variables play an important role in this method. Such variables have already been discovered in a limited class of models . Perhaps it is possible to extend them to other interesting models.

## Acknowledgment

I would like to thank C. Gattringer, C.B. Lang, J.C. Osborn, J.J.M. Verbaarschot and U.-J. Wiese for helpful comments. This work was supported in part by US Department of Energy funds under the grant DE-FG02-96ER40945. The computations were performed on BRAHMA, a Pentium based Beowulf cluster funded by grants from the Intel Corporation.


format: text

text:
# References

COMMISSIONS 27 AND 42 OF THE IAU

INFORMATION BULLETIN ON VARIABLE STARS

Number xxxx

Konkoly Observatory Budapest xx March 2002 HU ISSN 0374 – 0676

V2540 Oph (Nova Oph 2002): LARGE-AMPLITUDE SLOW NOVA

WITH STRONG POST-OUTBURST OSCILLATIONS

KATO, TAICHI1, YAMAOKA, HITOSHI2, ISHIOKA, RYOKO1

1 Dept. of Astronomy, Kyoto University, Kyoto 606-8502, Japan,
e-mail: (tkato,ishioka)@kusastro.kyoto-u.ac.jp

2 Faculty of Science, Kyushu University, Fukuoka 810-8560, Japan, e-mail: yamaoka@rc.kyushu-u.ac.jp

V2540 Oph (Nova Oph 2002) was independently discovered by Katsumi Haseda and Yuji Nakamura at magnitude 9.0 on 2002 January 24 (Haseda et al. (2002)). Retter et al. (2002) detected emission lines of hydrogen and Fe II, indicating that the object is an Fe II class nova caught in the early decline stage. Later examination of photographs revealed that the nova was already at magnitude 8.9 on 2002 January 19 (Seki et al. (2002)).

Since the detection of the outburst, the nova has been intensively monitored by a number of observers. Figure 1 shows the light curve constructed from visual, CCD $`V`$-band and photovisual observations reported to the VSNET Collaboration (http://www.kusastro.kyoto-u.ac.jp/vsnet/). The nova showed strong post-maximum oscillations up to 1 mag. The large-amplitude early stage oscillations resemble those observed in V1178 Sco = Nova Sco 2001) and V4361 Sgr = Nova Sgr 1996 (Kato, Fujii (2001)).

Superimposed on these oscillations, the nova showed a steady fade at 0.033 mag d-1. \[This rate was determined using the data between 2452294 and 2452341, during which the general trend of the fading can be approximated by a single decline rate. The last part of the light curve, when the nova underwent a long-lasting brightening, was not used in this analysis. If we incorporate the last part of the light curve, the average decline rate becomes 0.013 mag d-1, which may more severely constrain the following discussion\]. By applying the recently calibrated relation (della Valle, Livio (1995)) of absolute maximum magnitude vs rate-of-decline (MMRD) in classical novae, we obtain the expected absolute $`V`$-band maximum magnitude of $`M_\mathrm{V}`$=$``$7.0$`\pm `$0.5. We performed the accurate astrometry with the images obtained by Kyoto 0.30-m telescope taken on Mar. 8.23 UT, which revealed the position of the nova as: R.A. = 17h 37m 34s.385 $`\pm `$ 0s.017, Decl. = $``$16 23 18′′.19 $`\pm `$ 0′′.18 (equinox 2000.0, using 59 UCAC1 reference stars). This position is marginally consistent of the reported position by K. Kadota, who measured Haseda’s discovery films (Haseda et al. (2002)). No corresponding object was found on DSS and 2MASS scans within 2′′.5 of the nova, setting an upper limit of the prenova magnitude of $``$21.

This indicates that the lower limit of the outburst amplitude is $``$12.5 (by adopting the observed maximum magnitude of 8.5), which is unusually large for a slow nova with a decay rate of 0.033 mag d-1. By using the above expected absolute $`V`$-band maximum magnitude, we can set an upper limit of $`M_\mathrm{V}`$$``$5.5 for the nova progenitor. This magnitude is extremely faint for known prenova magnitudes and other novalike cataclysmic variables (Warner (1986, 1987)). \[Available observations suggest that the true maximum of the nova msut have been missed. By considering this, both the decline rate and the outburst amplitude could be larger than the values in this discussion. However, we consider this effect will not severely affect the conclusion, because 1) the MMRD-relation (della Valle, Livio (1995)) is known to be relatively flat (i.e. little depends on the decline rate) around the decline rate in question; a brighter maximum will therefore tend to pose a more stringent upper limit for the prenova), and 2) the reported spectrum (Retter et al. (2002)) suggests that the object was caught during an early decay stage.\]

Such a faint prenova magnitude would require a small mass-transfer rate, a small dimension of the disk, or a high inclination. Because V2540 Oph is a slow nova, the low mass-transfer rate is a rather unlikely explanation. We propose that the nova should have either a short orbital period or a high inclination. We know an excellent example of such a nova: V2214 Oph (Nova Oph 1988), which is believed to be a magnetic nova inside the period gap (Baptista et al. (1993)). V2540 Oph resembles V2214 Oph in many aspects: large outburst amplitude, slow rate of decline, and the presence of prominent oscillations. Since the characteristic double-wave orbital modulations were already present during the decay stage of V2214 Oph (Baptista et al. (1993)), we strongly encourage observers to detect orbital signatures in V2540 Oph.

We are grateful to all observers who reported vital observations to VSNET. This work is partly supported by a grant-in aid (13640239) from the Japanese Ministry of Education, Culture, Sports, Science and Technology. This research has made use of the 2MASS scan at NASA/IPAC Infrared Science Archive, the Digitized Sky Survey producted by STScI, and the VizieR catalogue access tool.


text:
# Chaotic dynamics in a storage-ring Free Electron Laser

## 1 Introduction

Since its advent, the laser has evolved as a unique tool both for fundamental investigations and important applications in optical science and technology . In providing coherent and intrinsically stable emission, lasers result fundamentally different from conventional light sources, which are characterized by a signal that is composite of random and uncorrelated emissions. Anyway, the view of the laser emission as an ordered and time-invariant process does not provide a complete description of the whole picture. During the last twenty years, profound mathematical discoveries have revolutionized the understanding of nonlinear science. Lasers, in particular, have been found to exhibit a rich variety of nonlinear dynamical behaviors, including unstable or even chaotic regimes. As in the case of other nonlinear dynamical systems, the transition from stable to chaotic dynamics is obtained by varying a control parameter, such like the losses of the cavity , its optical length  or the gain of the amplification process . The transitions are found to follow well-defined paths, regardless of the peculiar characteristics of the considered system. These universal “signatures” motivated both experimentalists and theoreticians in the search for physical systems exhibiting such phenomena, and in the investigation of their intrinsic similarity.
Concerning lasers, first experimental studies have been carried on a $`\text{CO}_2`$ laser , . It has been found that the modulation at a frequency $`f`$ of the chosen control parameter may result not only in oscillations with frequency $`nf`$, being $`n`$ an integer number, but also in a response at the sub-harmonics $`f/n`$. When the order $`n`$ of the sub-harmonics grows indefinitely, the response of the laser becomes irregular, though the system remains deterministic. A quite impressive similitude has been remarked between the “multi-furcation” diagram of the laser and that of the logistic map . In order to reproduce these experimental results, a simple theoretical model has been used, resting on the dynamical interplay between the electromagnetic field within the optical cavity, and the variables of the material employed as “active medium”. Theoretical work has been also done for demonstrating that the chaotic behavior of a laser can be stabilized either by using a self-controlling feedback procedure ,  or by a further modulation of a control parameter ,  .
Experimental evidences of a close link with deterministic chaos has also been given in the pioneering work  for a kind of non-conventional laser source: a storage-ring Free Electron Laser (SRFEL).
After the first operation of a FEL in the infrared spectral range on a linear accelerator in 1977 , a second FEL was installed on the storage ring ACO and provided the first visible radiation in 1983 . Present Storage Ring based FELs still supply the shortest FEL wavelength in the oscillator configuration (on ELETTRA ) and in the harmonic generation scheme (on Super-ACO ). User applications performed since 1993 on the Super-ACO FEL have demonstrated the very good quality offered by such sources in terms of tunability, high average power, rather short pulse duration and high degree of coherence . These characteristics make, for example, the FEL very suitable for two colours experiments performed in combination with the naturally synchronized synchrotron radiation .
A SRFEL, the principle of which is shown in Figure 1, is a coherent source of radiation in which the active medium consists of an ultra-relativistic electron beam, moving in a periodic magnetic undulator. The alternating magnetic field forces the electrons to move along sine-like trajectories and, consequently, to emit radiation, known as spontaneous emission. Once stored in the optical cavity, the radiation is amplified to the detriment of the kinetic energy of the electrons. This process generally leads to the increase of the rms value of the particles energy (the so-called electron beam energy spread) and, as a consequence, to the reduction of the amplification gain, until this latter reaches the level of the cavity losses. Since it origins from synchrotron radiation, the SRFEL output presents a micro-temporal structure with a period of the order of hundred of nanoseconds, which is determined both by the longitudinal dimension of the electron bunch and by the beam revolution period. On a larger (millisecond) temporal scale, the SRFEL dynamics depends strongly on the longitudinal overlap between the electron bunch(es) and the laser pulses at each pass inside the optical cavity (detuning condition). As it is shown in Figure 2 for the particular case of Super-ACO, a given detuning leads to a cumulative delay between the electrons and the laser pulses: the laser intensity may then appear “cw” (for a weak or strong detuning) or show a stable pulsed behavior (for intermediate detuning) , , . The narrow “cw” zone of the detuning curve (few ps around the perfect synchronism) is generally the most interesting for user applications. In fact, when in this zone, the laser is characterized by the maximum average power and the signal is the most close to the Fourier limit . In order to keep the laser-electron beam synchronism and avoid the jittering, which could determine a migration towards one of the unstable, pulsed zones of the detuning curve, efficient feedback systems have been implemented on the Super-ACO  and UVSOR FELs. SRFELs are complex, strongly coupled dynamical systems. The strong laser-electron beam coupling origins from the fact that, unlike a LINAC based Free Electron Laser, where the beam is renewed after each interaction, electrons are re-circulated. As a consequence, at every light-electron beam energy exchange, the system keeps memory of previous interactions.
The work that will be presented in this paper starts from the results obtained in 1990 on the Super-ACO FEL . In the latter reference, it was experimentally shown that a periodic modulation at a given frequency $`f`$ of the laser-electron beam detuning may lead to a period doubling of the laser intensity (i.e. to a laser response at a frequency $`f/2`$) or even to chaos for more important modulation amplitudes. Reference  contains also a first attempt to define the conditions (i.e. the amplitude and the frequency of the modulation) for which the bifurcation most likely occurs. The experimental data were discussed in connection with a simplified model based on a set of phenomenological rate equations, accounting for the coupled evolution of the laser intensity, of the electron-beam energy spread and of the gain of the amplification process. Such a model allowed to reproduce the basic features of the laser intensity evolution (in particular the saturation mechanism) only close to the perfect synchronism between the laser pulse and the electron bunch. This was a limiting factor for developing a comprehensive theoretical picture. In fact, the fundamental role played by the detuning has been only later understood: in  such effect has been shown to be responsible for the behavior of the laser intensity (“cw” or pulsed) on the millisecond scale and to produce a further gain reduction over that induced by the increase of the electron-beam energy spread. The simple model presented in reference  was then improved by including the detuning effect on the whole laser intensity distribution. This allowed to find a qualitative agreement between experiments and numerical results, for the case of the Super-ACO FEL .
It is worth to stress that the loss of the laser-electron beam longitudinal overlapping may be induced by different phenomena such like a vibration of the cavity mirrors at the line frequency  or a modulation of the electron beam energy .
In this paper an improved version of the model presented in  is used to get a deeper insight (with respect to that obtained in ) into the physics of deterministic chaos in a SRFEL. Numerical simulations are compared to a set of experiments performed on the Super-ACO FEL. In Section 2 the model is presented and shown to be able to reproduce quantitatively the observed features of the detuned FEL dynamics. Section 3 discusses the modifications induced to the dynamics by assuming an externally modulated detuning. Bifurcations, as well multi-furcation, are shown to occur for both the intensity and centroid position of the laser distribution. In Section 4 experimental results are presented. Finally, Section 5 contains concluding remarks and perspectives.

## 2 Theoretical model

The longitudinal dynamics of a SRFEL can generally be described by a system of rate equations accounting for the coupled evolution of the electromagnetic field and of the longitudinal parameters of the electron bunch .
The temporal profile of the laser intensity, $`y_n`$, is updated at each pass, $`n`$, inside the optical cavity according to:

$$y_{n+1}(\tau )=R^2y_n(\tau ϵ)\left[1+g_n(\tau )\right]+i_s(\tau ),$$
(1)
where $`\tau `$ is the temporal position with respect to the centroid of the electron bunch distribution; $`R`$ is the mirror reflectivity; the detuning parameter $`ϵ`$ is the difference between the electrons revolution period (divided by the number of bunches) and the period of the photons inside the cavity; $`i_s`$ accounts for the spontaneous emission of the optical klystron111 The case of SRFELs implemented on an optical klystron is here considered , . An optical klystron consists of two undulators separated by a dispersive section, (i.e. a strong magnetic field) favoring the interference between the emission of the two undulators.. The FEL gain $`g_n(\tau )`$ is given by:

$$g_n(\tau )=g_i\frac{\sigma _0}{\sigma _n}\mathrm{exp}\left[\frac{\sigma _n^2\sigma _0^2}{2\sigma _0^2}\right]\mathrm{exp}\left[\frac{\tau ^2}{2\sigma _{\tau ,n}^2}\right]$$
(2)
where $`g_i`$ and $`\sigma _0`$ are the initial (laser-off) peak gain and beam energy spread, while $`\sigma _n`$ and $`\sigma _{\tau ,n}`$ are the energy spread and the bunch length after the n$`th`$ light-electron beam interaction. The first exponential in the right-hand side of equation (2) accounts for the modulation rate of the optical-klystron spectrum222The laser-off peak gain has been optimized by assuming $`N+N_d=1/(4\pi \sigma _0)`$ in the expression of the modulation rate, where $`N`$ is the periods’ number of the undulators of the optical klystron and $`N_d`$ is the interference order due to its dispersive section ., while the second one reproduces the temporal profile of the electron bunch distribution. The bunch distribution is therefore assumed to keep its “natural” Gaussian profile under the action of the laser onset. This hypothesis entails that the interaction of the electron beam with the ring environment , ,  is neglected. This important point will be further discussed with particular concern to the case of Super ACO.
Defining $`g_{n,0}`$ as the peak gain after the n$`th`$ interaction, $`g_n(\tau )`$ can be written in the form:

$$g_n(\tau )=g_{n,0}\mathrm{exp}\left[\frac{\tau ^2}{2\sigma _{\tau ,n}^2}\right].$$
(3)
Figure 3 shows a schematic layout of the light-electron beam interaction in presence of a finite detuning $`ϵ`$.
The evolution of the normalized laser-induced energy spread $`\mathrm{\Sigma }_n=(\sigma _n^2\sigma _0^2)/(\sigma _e^2\sigma _0^2)`$ is given by:

$$\mathrm{\Sigma }_{n+1}=\mathrm{\Sigma }_n+\frac{2\mathrm{\Delta }T}{\tau _s}(I_n\mathrm{\Sigma }_n)$$
(4)
where $`\sigma _e`$ is the equilibrium value (i.e. that reached at the laser saturation) of the energy spread at the perfect tuning and $`\mathrm{\Delta }T`$ is the bunching period of the laser inside the optical cavity; $`I_n`$ is the normalized laser intensity defined as $`I_n=(1/I_e)_{\mathrm{}}^{\mathrm{}}y_n(\tau )𝑑\tau `$ (being $`I_e`$ the equilibrium value) and $`\tau _s`$ stands for the synchrotron damping time. Assuming that the saturation is achieved when the peak gain is equal to the cavity losses, $`P`$, the following relation holds333As function of the cavity losses, the mirror reflectivity $`R`$ is given by $`\sqrt{1P}`$.:

$$P=g_i\frac{\sigma _0}{\sigma _e}\mathrm{exp}\left[\frac{\sigma _e^2\sigma _0^2}{2\sigma _0^2}\right].$$
(5)
By inserting equation (5) in equation (3) and recalling the definition of $`\mathrm{\Sigma }_n`$, a closed expression for the peak gain is obtained:

$$g_{n,0}=g_i\frac{\sigma _0}{\sigma _n}\left[\frac{P}{g_i}\right]^{\mathrm{\Sigma }_n}\left(\frac{\sigma _e}{\sigma _0}\right)^{\mathrm{\Sigma }_n}.$$
(6)
Note that in the derivation of eq. (6) the variation of the bunch length during the FEL interaction has been taken into account. In this respect, the present model represents an improvement of the one proposed in , where the additional assumption $`\sigma _n=\sigma _0`$ in the definition of $`g_{n,0}`$ is made. A systematic numerical comparison allowed to show that this simplification may alter the period of the laser intensity in the pulsed regime, that is, as it will be shown in the next Sections, a parameter of paramount importance for the study the laser response to an external modulation.
Figure 2 shows a typical detuning curve obtained for the case of the Super-ACO FEL operated with only the main radio-frequency (RF) cavity. The parameters characterizing this configuration, which is the one of concern for the experiments reported in this paper, are listed in Table 1.

It is worth to mention that Super ACO can be also operated making use of an additional RF (harmonic) cavity , which has been installed with the aim of improving the FEL performances. Its effect results in a shortening of the electron bunch. The consequent increase of the electron density leads to the enhancement of the gain of the amplification process. On the other hand, the interaction of the electron beam with the ring environment (i.e. the metallic wall of the ring vacuum chamber) is also reinforced. This leads to a degradation of the electron beam quality (i.e. “anomalous” increase of the laser-off bunch length vs. current and deformation of temporal beam distribution) limiting in part the beneficial effect of electron-density increase . As regards the shape of the electron beam temporal distribution, the interaction with the ring vacuum chamber induces a perturbation of the “natural” Gaussian profile of the electron beam: the higher the cavity voltage, the stronger becomes the head-tail effect  deforming the electron distribution. It is worth to stress that, due to relative high value of the Super-ACO vacuum chamber impedance , the anomalous bunch lengthening and the perturbation of the beam profile are significant even when the harmonic cavity is passive. For the purpose of this paper, these perturbations can however be neglected in first approximation.
A quasi-symmetric beam profile reflects in a detuning curve that is almost symmetric for positive and negative detuning amounts (see Figure 2). The structure of the detuning curve has been studied by making use of the theoretical model presented above. The extension of the central (“cw”) and lateral (pulsed) zones, which are the “playground” of the experiments that will be discussed in the next Sections, has been found to be well reproduced (see caption of Figure 4). Figure 4 shows the regimes of the laser intensity for different detuning amounts. Again, the numerical results have been found in quantitative agreement with experiments. Concerning, in particular, the pulsed regime, a careful experimental and numerical analyses  pointed out that the period of the laser-intensity oscillations is not a constant (as it has been implicitly assumed in ) but it is instead an almost linear function of the detuning amount.
The effect of an imperfect light-electron beam longitudinal overlapping also influences the equilibrium position of the laser centroid. As shown in Figure 5 , the laser position is quite sensitive to small detuning amounts, while it attains an almost constant value as the side zones of the detuning curve are approached.
Numerical simulations are in agreement with experiments (see Figure 6) and allow to detail the interesting behavior of the laser centroid in the close-to-zero detuning region.

## 3 Numerical simulations for modulated detuning

Consider now the modification induced to the system dynamics by the cumulative effect of an external periodic modulation applied to the detuning $`ϵ`$, namely:

$$ϵ=a\mathrm{sin}(2\pi ft)+b,$$
(7)
where $`t`$ represents the time elapsed, integer multiple of $`\mathrm{\Delta }T`$; $`f`$ is the frequency of the oscillation, while the two amplitudes $`a`$ and $`b`$ control the maximum detuning amount. The oscillations are centered around $`b`$, being $`a`$ the maximum elongation. Tuning the value of $`b`$, allows to explore the effects induced by the modulation in different zones. Simulations have been performed based on the model previously introduced, the main motivation being a systematic comparison with preliminary results of experiments carried on in SuperACO. This point is addressed in the next Section. This paragraph is instead devoted to a more generally-oriented analysis, concerning the peculiar phenomena, that have been numerically found.
As already discussed in the previous Section, several different regimes are produced in a SRFEL, for constant detuning amounts. For small $`ϵ`$ the laser displays a “cw” intensity, while larger values of $`ϵ`$ drive the system into stable oscillations. This “natural” pulsed regime is of paramount importance for the present study. The right boundary of the central, symmetric, “cw” zone will be labeled as $`ϵ_{th}`$.
Consider first the case for which $`b`$ is set to zero. Hence, the average modulation is also zero. The external forcing, defined by equation (7), leads to a stable ($`1T`$) response only for small values of $`a`$, which prevent excursion of the modulated detuning outside the boundaries of the “cw” zone. When this excursion occurs, the modulation is generally found to induce a chaotic response of the laser intensity.
Quite different scenario is produced when $`b`$ is set to some small value different from zero. This choice corresponds to introducing a slight asymmetry which, surprisingly enough, is responsible of a significant regularization of the signal444 Note that this case can be considered more realistic than the previous one. In fact, even when the laser is operated as close as possible to the perfect tuning, a slight, uncontrollable, shift (corresponding to a $`b`$ value of the order of one or two fs in the case of Super ACO) has to be generally expected.. In Figure 7 the laser intensity is displayed as function of time, the upper panel being the modulation of $`ϵ`$. Simulations refer to a fixed value of the frequency while the amplitude is varied. The transition from a $`1T`$ to a $`2T`$ clearly occurs. It is worth stressing that the latter regime is observed for $`a>ϵ_{th}`$.
In Figure 8 the laser intensity is plotted as function of time. Again, simulations refer to different values of $`a`$, while $`b`$ and $`f`$ are maintained constant. For small values of $`a`$, the laser response is locked to the frequency of the external modulation, thus displaying a $`1T`$ regime. For larger values the laser intensity passes through a chaotic region and, finally, attains a stable $`3T`$ regime.
In Figure 9 the position of laser centroid is represented, as function of time, for the same choice of parameters as in Figure 7. Oscillations are displayed and a transition from a $`1T`$ to $`2T`$ regimes is observed. This is, indeed, a quite general result: each time a bifurcation of the laser intensity takes place, an analogous response is observed for the position of the centroid.
Summing up, a wide number of distinct behaviors are recovered and, within them, bifurcations of both intensity and centroid position of the laser, depending on the values of the parameters $`a`$ and $`f`$. For $`a<ϵ_{th}`$ , i.e. when the oscillation is confined in the central “cw” zone of the detuning curve (few ps), the system is locked to the frequency imposed by the modulation and laser intensity exhibits regular oscillations of frequency $`f`$. As soon as the value of $`a`$ exceeds the threshold $`ϵ_{th}`$, a cascade can occur. Therefore, one is led to conclude that excursions in the pulsed regime are a necessary requirement for the modulations to produce more complex behaviors. It is, in fact, believed that the observed periodic structures could be the result of the combined effect of the external and natural modulations. In addition, the system seems to be more sensitive to the natural frequencies (which, for the case of Super-ACO, lie between 400 Hz and 630 Hz), even though multi-furcations have been detected outside this privileged range. These observations mark an important point of distinction with similar studies on conventional lasers, being the pulsed regime an intrinsic characteristic of the SRFEL.
Another interesting feature regards the value, $`\widehat{a}`$, of the modulation amplitude for which the transition from $`1T`$ regime to more structured behaviors (i.e. $`2T`$, $`3T`$, …) occurs. Numerically it has been observed that increasing the frequency , the value of $`\widehat{a}`$ significantly decreases. This phenomenon, has been systematically observed in a wide range of frequencies (i.e. 200 Hz to 900 Hz). For larger values of $`f`$, the system is insensitive to the rapid external modulation, which therefore becomes ineffective.
Further, numerical simulations allowed to detect narrow windows, for which more complex periodic structures (i.e. $`4T`$, $`5T`$, $`6T`$, $`8T`$…) are displayed. Nevertheless, the existence of universal paths towards deterministic chaos remains an open issue.
Similar analysis have been carried out for larger values of $`b`$, such that the modulation is stably centered in one of the pulsed zones (see Figure 10). This condition is achieved for small enough values of the amplitude $`a`$. Numerical results produce remarkable features, which will be commented in the following. As already pointed out in Section 2, the frequency of the oscillation in the pulsed zone was shown to increase linearly with $`ϵ`$ . Assume $`f_p`$ to be the frequency of the natural oscillation which is found for $`ϵ=b`$, being $`a=0`$. When the modulation of $`ϵ`$ is switched on and the condition:

$$\frac{n}{f}=\frac{m}{f_p}n,mN,n<m,\frac{n}{m}Q$$
(8)
is fulfilled, the laser shows oscillations of frequency $`f/n`$, characterized by $`m`$ peaks, spaced with period $`2\pi /f_p`$ (see Figure 10, which refers to $`n=2`$ and $`m=3`$). Deviations are observed when condition (8) is no longer satisfied (for instance, by progressively changing the value of $`b`$). Note, however, that such behavior occurs for discrete values of $`f`$, and are not maintained over a finite set of frequencies. Hence, its characterization in term of multi-furcation is, somehow, stretched. Further experimental investigations will be addressed to validate these numerical findings.
Finally, an extensive campaign of simulations has also been performed varying the frequency $`f`$ of the modulation, while keeping $`a`$ fixed. A small value has been again assigned to the parameter $`b`$, in such a way that the center of the oscillations falls in the region of “cw” laser. A typical results is reported in Figure 11, where a cascade $`1T2T3T`$ is shown to occur. It is again worth to stress the crucial role that seems to be played by the pulsed zone.

## 4 Experimental results and their interpretation

Experiments have been performed by modifying the longitudinal overlap between the laser pulse and the electron bunch. This is done by modulating the radio-frequency in such a way that the variation of the detuning $`ϵ`$, in the form specified by equation (7), is achieved. Different results have been obtained by varying either the modulation amplitude (keeping $`f`$ constant) or the modulation frequency (for fixed value of the amplitude of the modulation).
In Figure 12 the evolution of the laser intensity is investigated, when adopting the first procedure: $`f`$ is set to a value close to the largest natural frequency. For small amplitudes of the modulation, a noisy $`1T`$ regime is observed. A transition toward a $`2T`$ response is produced when $`a`$ is larger than $`ϵ_{th}`$, in qualitative agreement with the general picture outlined in the previous Section. Larger values of $`a`$ induce a chaotic response, before a stable $`3T`$ regime is attained (see Figure 8, for a qualitative comparison).
In Figure 13, a $`1T2T`$ transition is observed when varying the modulation frequency. Again, the bifurcation has been found to occur for an amplitude of the modulation large enough to drive the laser in the pulsed zone of the detuning curve. The transition is observed when $`f`$ approaches the lowest natural frequency. This result is in a good agreement with the numerical simulations shown in Figure 11, which suggest that a further bifurcation $`2T3T`$ could have been found, for larger frequencies. In addition, the frequency for which the transition is shown to occur, increases when the amplitude of the modulation is reduced (see Figure 13). This result confirms the prediction of the theoretical analysis.
The response of the FEL to an external modulation of the laser-electron beam detuning has been also investigated by means of a double sweep streak camera. In particular, in Figure 14 the occurrence of a period doubling of the laser centroid is detected.

## 5 Conclusions

A complete analysis of the dynamics of a SRFEL in presence of a longitudinal laser-electron beam detuning has been performed. Numerical simulations, based on the model discussed in , and experimental results have been compared and shown to agree quantitatively, for a constant value of $`ϵ`$.
Further, the effect of an external modulated $`ϵ`$ has been considered. The system has been shown to display bifurcations, multi-furcations as well chaotic regimes for both the laser intensity and the position of its centroid, when either the amplitude or frequency of the modulation are tuned. Numerical analysis and preliminary experiments carried on for SuperACO, agree in this respect. Moreover, a detailed series of simulations, over a wide range of values of the modulation parameters, allowed to point out the crucial role played by the naturally pulsed zones of the SRFEL detuning curve. This observation, confirmed by a careful analysis of the experimental data, seems to indicate that the link with analogous investigations for conventional lasers is not trivial, the pulsed regime being an intrinsic characteristic of the SRFEL.
In addition, conventional laser are shown to produce a complete cascade towards deterministic chaos, similar to that of the logistic map. In the case of SRFEL, experiments proved clearly the existence of $`2T`$ and $`3T`$ regimes, as well of chaotic behaviors. Nevertheless, a progressive increase of the periodicity of the structures was not observed. Numerically, the interstitial regions laying between stable signals and chaotic regimes have been object of a detailed analysis. More complex periodic structures have been clearly detected ($`4T`$, $`5T`$, $`6T`$, $`8T`$…), even though, the attempt of identifying recurrent paths has, so far, failed. However, it is worthwhile stressing that the latter are shown to hold for a narrow range of tunability of the parameters. Hence, the transition from a low periodicity signal to chaos is generally sharp.
Finally, a complete understanding of the laser response in presence of an externally modulated detuning could allow to improve the performances of the FEL. It has in fact been proven in - that the chaotic behavior of a conventional laser can be successfully stabilized. It is planned to investigate deeper this possibility, being interested in the development of new self-controlling procedure.

## Authorship

The theoretical part of this work, described in section 2 and 3 above, was done by G.D.N. and D.F. These authors contributed equally to this part of the work. The experimental part, described in section 4 above, was done by C.B. and M.E.C.

## Acknowledgments

We thank M. Billardon, P. Glorieux and D. Dangoisse for stimulating discussions. We gratefully acknowledge the help of D. Garzella, D. Nutarelli S. Randoux, R. Roux and B. Visentin for carrying on experiments. We also thank E. Aurell for a careful reading of the manuscript. The work of D.F. was supported by the Swedish Research Council through grant NFR F 650 -19981250. Financial support to the work of G.D.N. was provided by the TMR contract N ERB 4061 PL 97-0102.

## Figures


text:
# The Band-Gap Problem in Semiconductors Revisited: Effects of Core States and Many-Body Self-Consistency

Present address: \]Department of Physics, University of California, Davis, California 95616.

## Abstract

A novel picture of the quasiparticle (QP) gap in prototype semiconductors Si and Ge emerges from an analysis based on all-electron, self-consistent, GW calculations. The deep-core electrons are shown to play a key role via the exchange diagram —if this effect is neglected, Si becomes a semimetal. Contrary to current lore, the Ge 3$`d`$ semicore states (e.g., their polarization) have no impact on the GW gap. Self-consistency improves the calculated gaps —a first clear-cut success story for the Baym-Kadanoff method in the study of real-materials spectroscopy; it also has a significant impact on the QP lifetimes. Our results embody a new paradigm for ab initio QP theory.

The modern “band-gap problem” originated with the realization that density-functional theory (DFT) Kohn , implemented in the local-density approximation (LDA), failed drastically in the description of the fundamental excitation gap of semiconductors and insulators. A significant step forward was achieved in the mid-eighties, when the first ab initio calculations of quasiparticle (QP) states were performed Hybertsen ; Godby  within Hedin’s GW approximation (GWA) Hedin . At the present time, it is nearly-universally accepted Hybertsen ; Godby ; Hedin ; Aulbur ; Shirley ; Rohlfing ; Rojas ; Hamada ; Arnaud  that the GWA yields QP gaps in semiconductors and insulators to within 0.1 eV of experiment —which is the level of accuracy required in the study of transport in these materials Shirley .

In this Letter we uncover a novel picture of the physical ingredients underlying the observed QP gap in Si and Ge. Central elements of this picture are the impact of the core electrons on the many-body problem for the states at the gap, and the role of self-consistency. The GW schemes alluded to above turn out to benefit from “cancellation of errors” involving the neglect of both effects. Our results illustrate the practical importance of the Baym-Kadanoff conserving method Baym  for the study of excitations in real materials.

These conclusions are arrived at by eliminating approximations which are routinely introduced in the implementation of Hedin’s scheme. First, the usual GW work invokes the pseudopotential (PS) approximation —by which the core states are effectively eliminated from the gap problem. However, PS theory does not guarantee that a “partitioning” of the electrons into two groups may lead to an accurate description of the dynamical self-energy of the valence states —which, according to DFT, is a non-linear functional of the total density. Semicore states pose a special challenge Shirley ; Rohlfing ; significantly, on the basis of a phenomenological model, the indirect nature of the Ge gap has been assigned to an effect of the polarization of the 3$`d`$ states Shirley .

Second, in most GW calculations the Dyson equation (DEq) is not solved to self-consistency. However, it has been shown, for a Hubbard-type model, that this practice leads to a genuine violation of charge conservation Schindlmayr . Still, from the available self-consistent solutions of the DEq Groot95 ; Holm98 ; Schone  it has been inferred that, while self-consistency, at the GW level, is a must in total-energy calculations Holm99 , the same is “to be avoided” in the study of spectroscopy Aulbur ; Groot95 ; Holm98 ; Groot96 . This state of affairs is unsatisfactory, since self-consistency is a necessary condition for the fulfillment of all the conservation laws Baym  —and, thus, for a proper theory of transport Faleev .

In our GW calculations all the electrons are taken into account in the evaluation of the valence-electron self-energy. Remarkably, the deep core states are found to play a significant role in the QP gap problem via the core-valence exchange diagram. An additional surprise is that the (presumably important) shallow Ge 3$`d`$ semicore states have no effect on the GW gap. Self-consistency at the GW level does improve the QP gaps; it also impacts the QP lifetimes. Our results for the Si gap and the indirect Ge gap al L agree with experiment very well. Other aspects of our calculated QP band structures provide signatures of physics beyond the GWA.

We recall that the exact self-energy is “$`\mathrm{\Phi }`$-derivable,” Luttinger  $`\mathrm{\Sigma }\left[G\right]\left(11^{}\right)=\delta \mathrm{\Phi }/\delta G\left(1^{}1\right)`$, where $`\mathrm{\Phi }`$ is the Luttinger-Ward “free-energy” functional; our notation stresses the fact that $`\mathrm{\Sigma }`$ is a functional of the dressed Green’s function $`G`$. Now, a $`\mathrm{\Phi }`$-derivable $`\mathrm{\Sigma }`$, obtained on the basis of an approximate $`\mathrm{\Phi }`$-functional, coupled with a self-consistent solution of the DEq, $`G^1(1,1^{})=G_0^1(1,1^{})\left(\mathrm{\Sigma }\left[G\right](1,1^{})\left(\mu \mu _0\right)\delta \left(11^{}\right)\right)`$, ensures that $`G`$ fulfills the conservation laws exactly Baym . \[Here $`1,1^{}`$ denote space-time points; Matsubara times $`\tau `$ are defined for $`0\tau \beta \mathrm{}`$. $`\mu `$, the chemical potential for the correlated system with a fixed number of electrons, is obtained self-consistently with $`G`$; $`\mu `$ differs appreciably ($``$1eV) Ku2000  from its counterpart $`\mu _0`$ for the reference one-electron system whose Green’s function is $`G_0`$.\] The GWA Hedin  is defined by the $`\mathrm{\Phi }`$-functional

$$\text{},$$
(1)
where the particle-hole bubbles are made up of $`G`$’s (not $`G_0`$’s), and the dashes represent the Coulomb interaction $`v`$. Functional differentiation of $`\mathrm{\Phi }_{GW}`$ yields $`\mathrm{\Sigma }_{GW}\left[G\right]=\mathrm{\Sigma }_H\left[G\right]+\mathrm{\Sigma }_{xc}\left[G\right]`$, where $`\mathrm{\Sigma }_H\left[G\right]`$ is the Hartree term, and the exchange-correlation (XC) term is of the Hedin form $`\mathrm{\Sigma }_{xc}\left[G\right](1,1^{})=G(1,1^{})W\left[G\right](1,1^{})`$, where $`W\left[G\right](1,1^{})`$ is the screened interaction Hedin .

We work in the basis of the Kohn-Sham (KS) states $`\varphi _{k,j}\left(x\right)`$ provided by the full-potential, linearized augmented-plane-wave (FLAPW) method Blaha ; here $`k`$ is a wave vector in the Brillouin zone, and $`j`$ is a band index. Adopting (without lack of generality) the KS system as the reference one-particle system, described within the LDA, the DEq can be written as diag

$`G_{k,j}\left(\tau \tau ^{}\right)=G_{k,j}^{LDA}\left(\tau \tau ^{}\right)+G_{k,j}^{LDA}\left(\tau \overline{\tau }_1\right)`$
$`\times \left[\mathrm{\Sigma }_{k,j}\left(\overline{\tau }_1\overline{\tau }_2\right)\left(V_{k,j}^{LDA}+\left(\mu \mu _0\right)\right)\delta \left(\overline{\tau }_1\overline{\tau }_2\right)\right]`$
$`\times G_{k,j}\left(\overline{\tau }_2\tau ^{}\right),`$ (2)
where $`G_{k,j}^{LDA}\left(\tau \right)=\frac{1}{\mathrm{}}e^{\epsilon _{k,j}\tau /\mathrm{}}\left(\theta \left(\tau \right)n_F\left(\epsilon _{k,j}\right)\right)`$, the KS eigenvalues $`\epsilon _{k,j}`$ being measured from $`\mu _0`$, $`V^{LDA}`$ is the KS potential without the nuclear contribution, and summation over variables with a bar on top is understood.

We solve Eq. (2) and the integral equation for $`W\left[G\right](1,1^{})`$ (in the latter case, in reciprocal space) on the $`\tau `$–axis. Our approach is ideally suited for a self-consistent evaluation of $`G`$ and should prove valuable in calculations beyond the GWA Ku2000 . A novel feature of our scheme —which is devoid of the cutoff effects encountered in $`\omega `$-axis formulations Schone — is the use of a non-uniform “power mesh” (PM) Ku2000  which, as outlined in Fig. 1, accounts for the nature of both $`G`$ and the particle-hole bubble —which are strongly peaked at the ends of the interval $`0\tau \beta \mathrm{}`$, being flat in between; moreover, the PM allows us to perform high-order interpolation (scaling $``$linearly) to generate all functions on the dense, uniform mesh required in the evaluation of the products entering the $`\tau `$-integrals Ku2000 . From the solution of Eq. (2) we evaluate the spectral function $`A_{k,j}\left(\omega \right)=\frac{1}{\pi }\text{Im}G_{k,j}\left(\omega \right)`$ via analytic continuation of $`G`$ onto the real-$`\omega `$ axis Ku99 . For the $`j`$-bands of interest, $`A_{k,j}\left(\omega \right)`$ shows a well-defined peak, whose $`k\omega `$ dependence yields the QP band structure —see Fig. 1.

The convergence of our results was tested by varying the parameters involved in: k-space sampling (we used 5x5x5 and 8x8x8 meshes), number of bands (14 and 24), number of (reciprocal lattice-) K-vectors used in the evaluation of $`\mathrm{\Sigma }_{xc}`$ (9, 27, 51, and 65), temperature (4000, 2000, 1000, and 300 K), and PM mesh (p5u20, p6u10, and p6u20). The most demanding parameter, the number of K-vectors kept in the valence exchange term, is associated with a monotonic opening of the Si gap. (The core contribution converges even slower; thus, a 6D real-space integral is evaluated instead.) In the case of Si we estimate that the absolute gap is converged to $``$0.1 eV from below; for Ge the convergence is even better, from both directions.

Table 1 provides the framework for a discussion of our results for the absolute QP gap (located at $``$80% $`\mathrm{\Gamma }`$X), the direct gap at $`\mathrm{\Gamma }`$, and the occupied bandwidth of Si. Our results comprise two levels of implementations of the GWA: $`(i)`$ A self-consistent solution of Eq. (2), corresponding to the evaluation of the self-energy $`\mathrm{\Sigma }_{GW}\left[G_{GW}\right]`$ (3rd row);  this calculation represents a numerical realization of a conserving approximation Baym . (ii) A non-conserving calculation (4th row) corresponding to the use of the self-energy $`\mathrm{\Sigma }_{GW}\left[G_{LDA}\right]`$ —the standard “$`G_{LDA}W_{LDA}`$” approximation Hybertsen ; Godby ; Hedin ; Aulbur ; Shirley ; Rohlfing ; Rojas ; Hamada ; Arnaud . Clearly, the gaps obtained from $`\mathrm{\Sigma }_{GW}\left[G_{GW}\right]`$ are in good agreement with experiment Landolt ; the mechanisms behind this agreement turn out to be quite unexpected —and instructive.

Indeed, on the basis of several numerical tests, we uncovered, first of all, the role of the deep-core electrons. In one set of calculations we suppressed their contribution to the exchange self-energy —2nd term in Eq. (1)— for the states at the gap exchange . At the $`\mathrm{\Sigma }_{GW}\left[G_{LDA}\right]`$ level, the absolute gap ($``$0.85 eV, see Table 1) is then reduced by 0.9 eV —i.e., the gap is closed, the ensuing QP band structure of Si corresponding to a semi-metal. While the size of this effect is surprising, its physics is easy to visualize: $`(i)`$The core electrons shield the attractive field of the nuclei (via the 1st term in Eq. (1), thereby raising the energy of the valence and conduction states; (ii) the exchange process partially compensates for this effect; (iii) the states across the gap have different amplitudes in the core region; these amplitudes control the strength of the exchange. It is the larger lowering of the energy of the QP states below the gap, relative to those above it, due to the non-local core-valence exchange process, that leads to this novel all-electron effect core .

Further insight into the role of the core electrons is obtained by simulating their contribution to $`\mathrm{\Sigma }_{xc}\left[G_{LDA}\right]`$ as $`\mathrm{\Sigma }_{xc}^{\text{from core}}V_{xc}^{LDA}\left[n_{tot}\right]V_{xc}^{LDA}\left[n_{val}\right]`$, where $`n_{tot}`$ and $`n_{val}`$ are the total and valence densities, and $`V_{xc}^{LDA}`$ is the XC contribution to $`V^{LDA}`$. This uncontrolled “LDA recipe” yields a spurious 0.15 eV additional opening of the Si gap Ku2000 . The significance of this test is that it accounts (Table 1) for the difference between the approximate FLAPW-based $`\mathrm{\Sigma }_{GW}\left[G_{LDA}\right]`$-level result of Hamada, Hwang, and Freeman Hamada  (1.01 eV) and our own (0.85 eV) —Ref. Hamada  relies on this LDA ansatz for the core contribution to the valence-electron self-energy, as the PS-based GW schemes implicitly do.

From Table 1 we draw a second key message: contrary to current wisdom Hedin ; Groot95 ; Holm98 ; Schone ; Groot96 , self-consistency does improve the quality of the calculated GW gaps of Si; cf. rows 1, 3, and 4. Indeed, the additional opening of the gaps obtained from $`\mathrm{\Sigma }_{GW}\left[G_{GW}\right]`$ brings them closer to their experimental values. This effect is traced to the dressing of both $`W`$ and $`G`$. Indeed, the dynamical screening built into $`W`$, which greatly reduces the exaggerated Hartree-Fock (HF) gap, is weakened for the dressed $`W`$; the dressing of $`G`$ widens the gap as well, a trend most easily visualized within self-consistent HF Ku2000 . We stress that the success of our conserving scheme is intimately related to the fact that we have carried out a full all-electron calculation exaggeration .

Interestingly, self-consistency also has a significant effect on the QP lifetimes, given by the inverse of the full-width at half-maximum (FWHM) of the QP peak in $`A_{\stackrel{}{k},j}\left(\omega \right)`$. As seen in Fig. 1d) for the hole state at the midpoint of the Si band, the FWHM (solid line) is significantly reduced, relative to its non-self-consistent counterpart (dashes). This effect —which has been ignored in the rapidly-growing literature on “hot carrier” lifetimes Echenique — recognizes two sources. First, fully dressed QP’s scatter less frequently off the Fermi sea than in the $`\mathrm{\Sigma }_{GW}\left[G_{LDA}\right]`$ case. Second, in the latter, non-self-consistent case, the gap edges recognized by $`G_{LDA}`$ and Im$`\mathrm{\Sigma }_{GW}\left[G_{LDA}\right]`$ are different. This mismatch in the gap edge, $`\mathrm{\Delta }`$, exaggerates the FWHM according to $`(\epsilon +\mathrm{\Delta })^2\epsilon ^2=2\mathrm{\Delta }\epsilon +\mathrm{\Delta }^2`$, where $`\epsilon `$ is the QP energy measured from each gap edge, and $`\mathrm{\Delta }`$\[$`\left(\mu \mu _0\right)`$ plus the shift of the respective edge relative to LDA\] Ku2000 . (N.B.: the QP states obtained from $`\mathrm{\Sigma }_{GW}\left[G_{LDA}\right]`$ have an unphysical finite lifetime, $`\mathrm{\Delta }^2`$, at the gap edges!)

From Table 1 it also follows that PS-based GW schemes carry a built-in error which is comparable with the LDA gap —cf. the representative PS-based result (1.29 eV Hybertsen ; 7th row) with our corresponding all-electron result (0.85 eV; 4th row). In addition to the non-local many-body effect of the core elucidated above, there is the impact on the matrix elements of $`\mathrm{\Sigma }_{GW}`$ of the removal of the oscillations of the valence wave functions at the atomic sites Ku2000 . Interestingly, the projector-augmented-wave result of Arnaud and Alouani Arnaud  (6th row) provides an independent measure of the error due to the latter source. We stress that, at the $`\mathrm{\Sigma }_{GW}\left[G_{LDA}\right]`$ level, the uncontrolled PS approximation Hybertsen ; Godby ; Hedin ; Aulbur ; Shirley ; Rohlfing ; Rojas  is masked by the neglect of self-consistency —which introduces an opposite-sign error, and results in apparent agreement with the experimental gap exaggeration .

Ge is isoelectronic with Si, and has the same diamond crystal structure. It is then reassuring that an analysis along the above lines confirms that our conclusions concerning the impact of the deep core states and many-body self-consistency apply in the case of Ge as well. Note, in particular, the 1st column of Table 2: The additional opening of the indirect gap at L obtained from $`\mathrm{\Sigma }_{GW}\left[G_{GW}\right]`$ (3rd row) —relative to the value obtained from $`\mathrm{\Sigma }_{GW}\left[G_{LDA}\right]`$; 4th row— leads to excellent agreement with experiment Landolt .

A significant issue in Ge is the role of the 3$`d`$ semicore states (which lie $``$25eV below the gap). In fact, in the standard PS-based GW approach, the correct topology of the QP band structure along $`\mathrm{\Gamma }`$L —which is automatically produced in our results— is obtained only upon introducing the 3$`d`$’s in the gap problem via a phenomenological core-polarization potential (CPP) model Shirley . Our all-electron approach places us in a position to address this issue. Shown in Fig. 2 is the density-response function Ku99  of Ge for a small wave vector, evaluated with and without the contribution from the 3$`d`$’s. Evidently, the screening is virtually identical in both calculations for energies up to $``$10 eV, by virtue of the weak transition strength of the 3$`d`$’s. Thus, the impact of the CPP model Shirley  is not physically justified. Moreover, if the 3$`d`$’s are excluded from both $`W`$ and $`G`$, our all-electron QP band structure remains unchanged (5th row). We conclude that the 3d’s play no role within the GWA.

Now, although our calculated indirect gap at L agrees with experiment very well, other aspects of our GW results for Ge are less satisfactory. As illustrated in Table 2 (3rd column), the empty states near X lie slightly below the lowest empty state at L —contrary to experiment Landolt . This ordering, which is reversed in the absence of screening —i.e., if $`\mathrm{\Sigma }`$ is evaluated within HF Ku2000 — provides a signature of the limitations of the GWA (note also our results for the Ge direct gap and the Si and Ge bandwidths). Thus, we expect that the localized Ge 3$`d`$ states may induce short-range correlation effects in $`\mathrm{\Sigma }`$; within the conserving method, such effects are to be included by adding an appropriate $`\mathrm{\Phi }`$-functional to Eq. (1) Eguiluz .

In summary, our results embody a new paradigm for ab initio QP theory, as we have demonstrated the non-trivial role played by the deep core states, and many-body self-consistency, in the QP gap problem; self-consistency —and thus, the fulfillment of the conservation laws— was also shown to impact the QP lifetimes. These effects, whose inclusion leads to excellent agreement with experiment for the Si absolute gap, and for the Ge indirect gap al L, are masked in the standard GW schemes. The Ge 3$`d`$ semicore states were found to play no role in the GW gap; their impact is likely to come via effects beyond the GWA. Our results indicate that a fundamental description of the entire valence QP band structure of Si and Ge within 0.1eV requires the inclusion of mechanisms beyond the GWA.

###### Acknowledgements.

This work was supported by NSF and NERSC.