text:
# The 2003 eclipse of EE Cep is coming

## 1 Introduction

EE Cep (BD+55°2693) is an extraordinary long-period eclipsing binary showing very large and complex changes of the shape and the depth of the minima (e.g. Baldinelli & Ghedini 1977). Quite recently Mikołajewski & Graczyk (MG (1999) \- hereafter MG) have suggested that the reason for these changes is the presence of a non-stellar eclipsing body in the system – a precessing, opaque and dark disc.

This star was identified as a variable star by Romano (roman56 (1956)) in 1952 and it was confirmed by Weber (weber56 (1956)) who reported observations of the previous minimum in 1947. The colour changes during the 1952 minimum were very small and the star was originally classified by Romano as an R CrB variable. Perova (pero57 (1957)) reported about 500 photographic observations of EE Cep between 1899-1908 and 1933-1955 but no brightness changes were recorded. Romano & Perissinotto (roman66 (1966)) have suggested that it may be an eclipsing binary but they did not give any approximate ephemeris. Its binary nature was confirmed later by Meinunger (mein73 (1973)) who analysed over 500 photographic plates from Sonneberg Observatory. The period of the star is almost 2050 days and only the primary minimum can be detected. Meinunger has not detected any changes of the brightness between the 1958 and 1964 minima. The brightness of the star outside of the eclipse phases is constant within the observational errors: $`V_{\mathrm{max}}=10.78`$ (Meinunger & Pfau mein81 (1981), MG). The primary component is apparently a normal B5 III star (Herbig her60 (1960), Baldinelli et al. bal81 (1981)). The emission visible in $`H_\alpha `$ and $`H_\beta `$ lines arises most likely from an invisible companion (Brückner bryk76 (1976)). The absorption spectrum of the companion is never seen.

The early model of EE Cep presented by Meinunger (1976) assumed that the companion should be a cool star of M type, much less luminous but larger than the primary B5 III star, surrounded by a strongly extended envelope. The pulsations of the cool star would change its diameter and thus lead to the changes of the minima duration and depth. The envelope would cause an atmospheric eclipse in the beginning and in the end of the minimum, producing typical wings in the light curves. However, such a simple model cannot explain some of the features observed in EE Cep, as was argued by MG and in this paper.

## 2 Photometric behaviour

### 2.1 A review of past eclipses

The typical minimum of EE Cep has large wings, the duration of the eclipse is $`D40`$ days or a phase interval $`\mathrm{\Delta }P=0.02`$, and the depth of the eclipse is moderate $`\mathrm{\Delta }B1.5`$ mag. Such minima were observed in 1975, 1980 and 1997. But some eclipses observed up to date deviated strongly from this picture. We have collected observations of all previously observed eclipses of EE Cep in two roughly similar photometric wavelength bands: 1) $`B`$ and $`m_{\mathrm{pg}}`$; 2) $`V`$, $`m_{\mathrm{pv}}`$ and $`m_{\mathrm{vis}}`$. The references used for the photometry are listed in Table1, whereas the data are presented in Figs. 1 and 2. Part of the historical data were published only as diagrams and they needed to be read from these figures. All observations were phased using the ephemeris given by MG.

The 1947 minimum ($`E=1`$) was the first one detected and its depth was 1$`\stackrel{m}{.}`$5 in $`m_{\mathrm{pg}}`$. There is no information about the duration and shape of the eclipse and it is not shown in Fig. 1. The 1952 minimum ($`E=0`$) was deep (at least 1.7 mag in $`m_{\mathrm{pg}}`$) and has an asymmetric shape. During the next eclipse, in 1958 ($`E=1`$), there are only two observational points and two upper limits $`13`$$`\stackrel{m}{.}`$in $`m_{\mathrm{pg}}`$ at the bottom of the minimum (Romano & Perissinotto 1966), therefore we also omitted this minimum in Fig. 1. Nevertheless, it could be very similar, or maybe a little deeper, to the following one in 1964 ($`E=2`$). The depth of these two minima is at least 2 mag. The duration of the 1964 eclipse was only 30 days.

During the 1969 ($`E=3`$) eclipse, the minimum was, in contrast to the previous one ($`E=2`$), extremely shallow ($`\mathrm{\Delta }m_{\mathrm{pg}}=0.6`$ mag) and a flat bottom phase was clearly observed. Also, the 1969 eclipse was probably the longest one – its total duration was about 60 days. The 1975 ($`E=4`$) and 1980 ($`E=5`$) eclipses were more or less typical, while the 1986 ($`E=6`$) minimum was again shallow and strongly asymmetric. The 1992 ($`E=7`$) eclipse was observed visually by three amateur astronomers from AAVSO and seven $`VR`$ measurements with a CCD camera have been obtained by J. Borovicka from Ondrejov Observatory in the Czech Republic (Halbach hal92 (1992)). Unfortunately, there is a zero-point offset of +0$`\stackrel{m}{.}`$16 between Borovicka’s $`V`$ magnitude scale and Meinungers’s (mein76 (1976)) and MG’s scales. Also, the visual light curves strongly suffer from differences in zero-point scales. Thus we have shifted visual estimations to the corrected CCD observations, although three points still lie out of the light curve (Fig. 2). This minimum was shallow ($`\mathrm{\Delta }m_{vis}0.8`$) and probably the shortest one – its total duration was about 22 days. The last observed 1997 ($`E=8`$) eclipse was, in contrary, deeper than usual. However, its shape and the total duration were similar to those of 1975 and 1980.

### 2.2 The characteristic phases of the minimum

In general we can note that 1) the shape of the eclipse may change from the deep and narrow minimum (like in 1964) to a shallow and very wide minimum (1969, 1986), 2) a flat-bottom phase may occur which is characteristic for annular eclipses, 3) wide wings are often visible and 4) there is a persistent asymmetry of the light curve between the descending and the ascending branches of the minima.

Although the light curve of EE Cep shows these changes, it is possible to single out the most typical shape of the minimum. MG suggested that two eclipses from 1975 and 1980 are most typical. Basing on the data from these two eclipses they construct the representative (”template”) eclipse shown in Fig. 3. There are some typical features of the minimum which are also observed in other eclipses of EE Cep. The most important is the presence of six characteristic moments during an eclipse. Moments 1st, 2nd, 3rd and 4th correspond to contacts observed in normal eclipsing binaries, while moments $`1^{}`$st and $`4^+`$th denote the start and the end of the phase of semi-atmospheric eclipse, respectively. Another feature is the asymmetry of the light curve between the ascending and descending branches of the minimum. Such asymmetry is clearly visible in eclipses of epochs 4, 5 and 8 and may be supposed to occur in the eclipse of epoch 6. The asymmetry seems to be caused by the presence of the sloping-bottom phase (between 2nd and 3rd contacts – Fig. 3) which influences the shape of the ascending branch. This picture is supported by the fact that during eclipse of epoch 3 four distinct contacts were observed and the bottom was flat. Thus the part of the light curve between the 2nd and 3rd moments is not part of the egress.

The characteristic shape and phases of the eclipse are independent of the photometric band. Over a broad spectral range, from $`U`$ until $`I`$, the eclipse seems to be almost grey (MG). During three eclipses (E = 4, 5, 6) the depth in the $`B`$ band was on average $`0.07`$ mag deeper than in the $`V`$ band, possibly due to a selective extinction of the eclipsing body matter.

### 2.3 Depth – duration relation

Table 1 presents the total duration and the depth of EE Cep’s eclipses estimated from data in Figs. 1 and 2. Figure 4 presents the relation between the duration and the depth of eclipses in $`B/m_{\mathrm{pg}}`$. For epochs 6 and 7 we have only $`V/vis`$ estimates so we applied a correction of $`\mathrm{\Delta }m=+0.07`$ mag to the depth of these two eclipses. It is worth noting that there is a quite distinct correlation between these two parameters – longer eclipses correspond to shallower minima. Also, from the inspection of Figs 1 and 2 one can note that the shallowest minima have the most asymmetric wings. The minimum of Epoch 7 seems to be an exception to the rule. However, the estimates of $`\mathrm{\Delta }P`$ and $`\mathrm{\Delta }V`$ for this eclipse are based mainly on inconsistent visual observations, this may reflect the presence of a more complex relation than simple linear dependence.

The value of $`\mathrm{\Delta }B/m_{\mathrm{pg}}`$ depths versus the number of epoch is presented in Fig. 5. Its evolution is rather complex, but a possible period of nine epochs ($`50`$ years) can be suggested. Because it is almost equal to the timespan of all the historical data, new observations are needed to confirm its existence. The duration of the eclipses, which is a more difficult parameter to measure, do not show such clear time evolution.

## 3 The model of disc eclipses

The model presented by MG assumes that the invisible component is an optically thick, dark disc inclined to the orbital plane – a “cousin” of the disc in $`ϵ`$ Aur system (e.g. Carroll et al. 1991). The spin axis precession of the disc changes two angles: 1) inclination of the disc to the line of sight (which causes changes of the projected area of the disc) and 2) inclination of the projected disc to the orbital movement direction (tilt). As a result we observe varying depth of the minima as it is shown in Fig. 4 in MG. Another kind of a precession – the nodal regression – was recently proposed also for two relatively close eclipsing binaries: SS Lac (Torres 2001) and V907 Sco (Lacy et al. 1999) as a explanation of a varying eclipse depth. In both cases there are precession motion (timescale of hundred years) of the whole short-period (several days) binary orbit caused by an interaction with a third body moving on a wider orbit with a period of hundreds days. In the case of the long period binary EE Cep, however, we observe only the precession of the disc itself.

During our quantitative analysis and after doing some simple numerical tests, we have found that the most probable reason for the observed asymmetries was not the asymmetry of the disc itself but the tilt of the projected disc to the direction of motion. The tilt produces typical asymmetry of the minima as the result of the sloping-bottom between inner contacts (Fig. 6). The wide wings between $`1^{}`$st and 1st, and between 4th and $`4^+`$th may be caused by the outer, semitransparent part of the disc. The “true” eclipse is caused by the opaque, inner part: the four characteristic moments 1st–4th may be easily understood as contacts between the outer and inner parts of the tilted disc with the eclipsed star – see Fig. 3 and 6. During the extreme, rare, edge-on and small tilt configuration the eclipses should have a flat bottom phase as was observed in 1969. The possible period of nine epochs ($``$ 50 years) may reflect precession of the disc.

## 4 Spectroscopic behaviour

There are very few spectroscopic observations of EE Cep described in the literature. Herbig (her60 (1960)) gave its spectral type as B5:ne$`\beta `$. Brückner (bryk76 (1976)) noted that the emission profile in $`H_\alpha `$ is not typical of Be stars and during the eclipse of epoch 4 there was a considerable strengthening of this emission in respect to the continuum. However, his observations were done in very low resolution (200 Å/mm). Baldinelli et al. (bal81 (1981)) did not record any changes in the spectrum during the following eclipse (E=5), in spite of the slightly higher resolution (60 Å/mm). Moreover, they noted broad absorption in the Balmer series from $`H_\alpha `$ to $`H_8`$ with strong, broad ($`H_\alpha `$) and narrow ($`H_\beta `$) emission lines superimposed. This spectrum was estimated as a typical Be star of B5 III spectral type, which agrees with the $`BV`$, $`UB`$ colours analysis (MG). Thus, the connection between Balmer emission and eclipsed (Be star) or eclipsing (disc) body is not clear.

We have secured several spectra covering the $`H_\alpha `$ and the sodium doublet Na I regions of EE Cep. The spectra were taken at NAO Rozen in Bulgaria on a 2-m telescope equipped with a coude-spectrograph + LN2-cooled CCD camera at resolution 0.4 Å. Unexpected large changes of the $`H_\alpha `$ profile can be noted outside the eclipse phases on a time scale of several months to years (Fig. 7). The $`H_\alpha `$ profiles are strong and broad ($`600`$ km/s), show one or two peaked emission with variable asymmetries. Two profiles obtained during two consecutive days are almost identical, however a red-shifted bulge around 60 km/s seems to be stronger in the later spectrum. The sodium doublet D1D2 shows strong single-component, probably interstellar absorptions. However, its heliocentric radial velocity probably changed from $`V_\mathrm{r}=36`$ km/s in September 1998 to $`V_\mathrm{r}=24`$ km/s in November 2002. Also the equivalent width of D1D2 doublet seems to be variable, decreasing about 10% between September 1998 and November 2002. The mean value of the D1 equivalent width $`EW=0.57`$ Å corresponds to $`E_{\mathrm{B}\mathrm{V}}=0.46`$ using the calibration given by Munari & Zwitter (mun97 (1997)), which is fairly consistent with the reddening derived from $`UBV`$ photometry (MG).

## 5 The next eclipse: 2003 ($`E=9`$)

The moment of the mid-eclipse according to the elements given in MG is JD 2452794 (3 June 2003). As shown in Fig. 4 the eclipse can start even about 30 days before mid-eclipse. The reason for it is the presence of the wide wing during ingress which can be usually observed, thus the photometric observations should start at least 40 days before the expected mid-eclipse. The brightness of the star should decrease by at least 1 mag in $`B,V`$ filters, but the observers should keep in mind that this decrease may be much larger: up to 2 mag! Also, a significantly shallower (as epochs 3 and 7) minimum can occur. If the periodicity of nine epochs suggested by Fig. 5 is real we can expect an eclipse similar to that of E=0, i.e. narrow and deep. Near-far infrared photometry during the eclipses will be very important as it would give the answer to what is the temperature of the companion. Also, good quality $`B,V`$ or $`R`$ light curves from the whole eclipse will be valuable to make numerical tests of the disc model.

The spectroscopic observations appear to be very important. We should expect the discovery of a shell spectrum produced by the gaseous, semitransparent part of the disc like in $`ϵ`$ Aur (Ferluga & Mangiacapra 1991). Moreover, such observations would give us an additional clue to the answer of the origin of the variable $`H_\alpha `$ line. The systematic observations would be very useful to obtain the spectroscopic orbit of the visible B5 component.

The eclipsing body in EE Cep system seems to be unique. The confirmation (or rejection) of the disc hypothesis can be attained only by detailed photometric and spectroscopic study of the incoming eclipses.

###### Acknowledgements.

This paper was supported by KBN Grant No. 5 P03D 003 20 and BAN/PAN Joint Research Project No.28. We would like to thank: Ms Karolina Wojtkowska for her help in obtaining some rare literature, dr Alojzy Burnicki for his help in translation of some German publications and dr Boud Roukema for the very careful reading of the manuscript and for the English corrections in the text.


format: text

text:
# Phase Transitions of a Bose Gas in an Optical Lattice

## 1 Introduction

In this paper an interacting Bose gas in an optical lattice is studied. It is based on some previously published ideas for the hard-core Bose gas and its generalization to soluble models. The application to Bose-Einstein condensation, depletion of the condensate and the Mott insulating state in these systems is discussed.

A real Bose gas, as has been studied in the case of trapped cold alkali atoms , is always subject to interparticle interaction. In a good approximation this interaction is a hard core, characterized by a scattering length $`a`$. Depending on the type of atoms and external conditions (e.g. absence or presence of an external magnetic field), the scattering length can vary over a wide range. Typical values of $`a`$ observed in experiments are $`a5`$nm in 87Rb  and up to $`a500`$nm in a gas of 85Rb atoms near the Feshbach resonance . Negative scattering lengths are also known but will not be considered here. The effect of the hard-core interaction can be characterized by the interaction parameter

$$I=\frac{\mathrm{scattering}\mathrm{length}}{\mathrm{typical}\mathrm{interparticle}\mathrm{distance}}.$$

If the typical interparticle distance in a Bose gas is large (e.g., $`200`$nm in a very dilute gas that undergoes a Bose-Einstein condensation ), $`I`$ is small ($`I0.01`$), and the interaction can be treated as a weak perturbation. This is the regime where the Gross-Pitaevskii physics can be used .

If the Bose gas is transferred to an optical lattice the characteristic length scale is given by the distance of the potential minima, the lattice constant $`a_l`$. Typical experimental values are $`a_l250\mathrm{}500`$nm. Then the interaction is $`I0.01\mathrm{}1`$.

¿From the theoretical standpoint the simplest case is the non-interacting Bose gas ($`I=0`$), where no qualitative effect of the lattice on the Bose-Einstein condensation is expected. Indeed, Bose-Einstein condensation in a non-interacting Bose gas is not a phase transition, characterized by spontaneous symmetry breaking, but a macroscopic occupation of the groundstate. A grand-canonical ensemble of non-interacting bosons at inverse temperature $`\beta `$, with quantum numbers $`l=0,1,\mathrm{}`$ and energy $`ϵ_l0`$ is given by the partition function

$$Z_0=Tre^{\beta \widehat{H}}=\underset{l0}{}\underset{n_l0}{}(\zeta e^{ϵ_l})^{\beta n_l}.$$
(1)
The fugacity $`\zeta `$ is chosen such that $`ϵ_0=0`$ and $`ϵ_l>0`$ for $`l>0`$. It is assumed that the groundstate is unique. Then $`\zeta >1`$ is not allowed here because this would lead to a divergent series in Eq. (1) for $`l=0`$. The average number of particles with quantum number $`l`$ is

$$N_l=n_l=\frac{1}{\beta }\frac{}{ϵ_l}\mathrm{log}Z=\frac{1}{e^{\beta ϵ_l}\zeta ^\beta 1}.$$
(2)
This quantity decreases with increasing $`ϵ_l`$ and has its maximum at $`l=0`$. The probablity of finding a particle in the groundstate is

$$\frac{N_0}{N}=1\frac{1}{N}\underset{l1}{}N_l.$$
(3)
The Bose-Einstein condensate (BEC) is the macroscopically occupied groundstate, i.e., $`lim_N\mathrm{}N_0/N>0`$. According to the result (2) this can only happen if $`\zeta `$ goes to 1 when $`N`$ goes to infinity. Using the density of states $`\rho (ϵ)`$ at energy $`ϵ`$ and keeping the ratio $`N/V`$ ($`V`$: volume of the Bose gas) fixed when $`N\mathrm{}`$, the second term in Eq. (3) becomes the integral

$$\frac{V}{N\beta }_0^{\mathrm{}}\frac{\rho (x/\beta )}{e^x1}𝑑x.$$

There exists a BEC if the integral is smaller than 1. The latter is the case if (i) the density of states $`\rho (ϵ)`$ cuts-off the pole of the integrand (this is usually the case for three-dimensional lattices) and if (ii) $`V/(N\beta )`$ is sufficiently small. Thus the presence of a lattice does not prohibit the formation of a BEC in a non-interacting Bose gas at sufficiently low temperatures. More interesting physics is expected for $`I1`$ when the lattice constant is of the same order as the scattering length because the interaction can block the tunneling processes and lead to an insulating state without phase coherence and with a locally fixed number of bosons, the Mott insulator.

## 2 Random-Walk Representation

In order to discuss the interaction between bosons in an optical lattice it is convenient to derive a random-walk representation for the lattice Bose gas. This is related to the idea that the average over particle fluctuations in an equilibrium system is equivalent to the average of random walks. For this purpose one can start from non-interacting bosons, then develop the weakly (Bose-Hubbard-like) and finally the strongly (hard-core) interacting gas.

### 2.1 Non-interacting Bose Gas

Evaluation of the geometric series in the partition function of Eq. (1) gives

$$Z_0=\underset{l}{}\frac{1}{1\zeta ^\beta e^{\beta ϵ_l}}.$$

Assuming that $`\beta `$ is an integer, the identity

$$1\zeta ^\beta e^{\beta ϵ_l}=det(\mathrm{𝟏}\zeta w_l)$$

with the $`\beta \times \beta `$ matrix

$$w_l=\left(\begin{array}{ccccc}0& u_l& 0& \mathrm{}& 0\\ 0& 0& u_l& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& 0\\ 0& \mathrm{}& 0& 0& u_l\\ u_l& 0& \mathrm{}& 0& 0\end{array}\right)(u_l=e^{ϵ_l})$$
(4)
and the diagonal $`\beta 𝒩\times \beta 𝒩`$ matrix $`𝐰=diag(w_1,w_2,\mathrm{},w_𝒩)`$ can be used to write the partition function as

$$Z_0=\frac{1}{det(\mathrm{𝟏}\zeta 𝐰)}.$$
(5)
The matrix elements in the last equation are defined with respect to $`x=(l,t)`$ ($`t=1,2,\mathrm{},\beta `$, $`l=0,1,\mathrm{},𝒩`$). The diagonal matrix $`𝐰`$ can be replaced by any other matrix to describe a non-interacting Bose gas. For instance, tunneling between sites $`r,r^{}`$ of the three-dimensional optical lattice gives with $`x=(r,t)`$ the matrix elements

$$w_{x,x^{}}=\{\begin{array}{cc}J\delta _{t^{},t+1}/6\hfill & \text{for }r,r^{}\text{ nearest neighbors}\hfill \\ (1J)\delta _{t^{},t+1}\hfill & \text{for }r^{}=r\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}.$$
(6)
This choice implies that a particle either stays in a potential well of the optical lattice at $`r`$ with probability $`1J`$ ($`0J1`$) or tunnels from a well at $`r`$ to one of the six neighboring wells with probability $`J/6`$. An expansion of $`Z_0`$ in powers of $`\zeta `$

$$Z_0=\underset{n0}{}\frac{(1)^n}{n!}\left[\underset{l1}{}\frac{\zeta ^{\beta l}}{l}\frac{1}{\beta }Tr(𝐰^{\beta l})\right]^n$$
(7)
has an interpretation in terms of random walks: the matrix elements $`w_{x,x^{}}`$ can be considered as (random) steps in positive time direction between the lattice sites at $`r`$ and $`r^{}`$. Then the trace $`Tr(𝐰^{\beta l})`$ consists of terms which start at a site $`r`$, go $`l\beta `$ steps and eventually return to $`r`$. This is a consequence of the fact that $`𝐰`$ has non-zero matrix elements only with $`t,t+1`$. For $`l>1`$ the walk is sent back to $`t=1`$ after $`\beta `$ steps due to the periodic boundary conditions defined by the matrix (4). A realization of random walks is shown in Fig.1. The number $`l1`$ can be associated with the winding number of the corresponding walk.

### 2.2 Bose Gas with Intermediate Interaction

A weakly-interacting Bose Gas is described by the Bose-Hubbard model . This can be derived from the random-walk expansion of Eq. (7) by adding an interaction between the random walks. Since the interaction is local and repulsive, the statistical weight of intersecting walks must be reduced in comparison with their non-intersecting counterparts. This can be implemented in the random-walk representation of Eq. (7) by introducing a phase factor to the tunneling matrix element

$$w_{x,x^{}}w_{x,x^{}}e^{2i\phi _x}v_{x,x^{}}.$$

and using for the partition function of the interacting system

$$Z_I=det(\mathrm{𝟏}\zeta 𝐯)^1_\phi $$

with the Gaussian distribution of $`\phi _x`$ with zero mean, and $`\phi _x^2_\phi =U`$. The expansion of $`Z_I`$ in powers of $`\zeta `$ yields terms of the type

$$\underset{x^{}}{}(𝐯_{x,x^{}})^{l_x^{}}_\phi =\underset{x^{}}{}(𝐰_{x,x^{}})^le^{U(_x^{}l_x^{})^2}.$$

For a singly occupied $`x`$ (i.e. $`_x^{}l_x^{}=1`$) the averaging gives a rescaling of the fugacity $`\zeta \zeta e^U`$. In the case of a single well (i.e., $`J=0`$ in expression (4)) the partition function can be easily evaluated and gives with $`\zeta =e^\mu `$ ($`\mu `$ being the chemical potential)

$$Z_I=\left(1\zeta ^\beta \underset{t}{}e^{2i\phi _{t,t+1}}\right)^1_\phi =\underset{n0}{}e^{\beta U(n\mu /U)n}.$$

The summation over $`n0`$, here formally obtained from a geometric series, corresponds with the summation over particle numbers, as one can see by comparing with Eq. (1). For large $`\beta `$ only the term with minimal value of $`(n\mu /U)n`$ contributes significantly. This is a (discrete) saddle point approximation:

$$\underset{n0}{}e^{\beta U(n\mu /U)n}\underset{n0}{\mathrm{max}}e^{\beta U(n\mu /U)n}.$$

The value of $`n`$ that gives the maximum is also the average number of particles. To evaluate this quantity one has to determine the groundstate energy

$$E_0(\mu /U)=\underset{n0}{\mathrm{min}}U(n\mu /U)n$$

for a given value of $`\mu /U`$. The energy $`E(\mu /U,n)=U(n\mu /U)n`$ as a function of $`\mu /U`$ is linear. This implies for the groundstate energy a piecewise linear function with slope $`n`$ for $`2n1\mu /U2n+1`$ (cf. Fig. 2). If one includes tunneling processes, e.g. using an expansion in small $`J`$, the phase diagram would give the lobe structure, discussed for the Bose-Hubbard model .

### 2.3 Hard-core Bose Gas ($`I=1`$)

In the case of hard-core interaction it is not just that random walks have to pay a penalty when they intersect but intersections are prohibited completely. This is a constraint to the random walks which can be included explicitly by using the functional integral representation. Starting from $`Z_0`$

$$Z_0=\frac{1}{det(\mathrm{𝟏}\zeta 𝐰)}=\mathrm{exp}\left(\underset{x,x^{}}{}\mathrm{\Phi }_x(\mathrm{𝟏}\zeta 𝐰)_{x,x^{}}\mathrm{\Phi }_x^{}^{}\right)\underset{x}{}d\mathrm{\Phi }_xd\mathrm{\Phi }_x^{}/\pi ,$$

the complex field $`\mathrm{\Phi }_x`$ is replaced by a nilpotent field: $`\mathrm{\Phi }_x\eta _x^1`$, $`\mathrm{\Phi }_x^{}\eta _x^2`$, where $`(\eta _x^j)^l=0`$ for $`l>1`$. With the integral over the nilpotent field

$$_{HC}\underset{x\mathrm{\Lambda }^{}\mathrm{\Lambda }}{}\underset{\sigma j_x}{}\eta _x^\sigma =\{\begin{array}{cc}1\hfill & \text{if }\mathrm{\Lambda }^{}=\mathrm{\Lambda }\text{}j_x=1,2\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array},$$

where $`\mathrm{\Lambda }`$ is the full lattice spanned by $`x`$, the partition function of the hard-core Bose gas reads

$$Z_{HC}=_{HC}\mathrm{exp}\left[\underset{x,x^{}}{}\eta _x^1(\mathrm{𝟏}+\zeta 𝐰)_{x,x^{}}\eta _x^{}^2\right].$$

The expansion of $`Z_{HC}`$ in powers of $`\zeta `$ produces the non-intersecting random walks thanks to the nilpotent condition of the $`\eta `$’s.

Physical interesting quantities like the density of particles $`n`$ are expressed directly in terms of $`Z_{HC}`$:

$$n=\frac{\zeta }{\beta 𝒩}\frac{\mathrm{log}Z_{HC}}{\zeta }.$$
(8)
The groundstate of a Bose gas can exhibit off-diagonal long-range order that is characterized by a non-vanishing condensate density $`n_0`$, obtained from a correlation function that is related to $`Z_{HC}`$:

$$\underset{|r|\mathrm{}}{lim}\frac{1}{Z_{HC}}\frac{1}{\beta }\underset{t}{}_{HC}\eta _0^1\eta _{(r,t)}^2\mathrm{exp}\left[\underset{x,x^{}}{}\eta _x^1(\mathrm{𝟏}+\zeta 𝐰)_{x,x^{}}\eta _x^{}^2\right]=n_0>0.$$

## 3 Fermion Representation of Hard-Core Bosons

The nilpotent variables $`\eta _x^j`$ can be replaced by products of Grassmann variables as $`\eta _x^\sigma (1)^\sigma \psi _x^\sigma \overline{\psi }_x^\sigma `$. This gives with the usual Grassmann integration  the partition function

$$Z_{HC}=e^{S_{\mathrm{ferm}}}\underset{x,\sigma }{}d\psi _x^\sigma d\overline{\psi }_x^\sigma $$

with

$$S_{\mathrm{ferm}}=\underset{x}{}(\psi _x^1\psi _x^2+\overline{\psi }_x^1\overline{\psi }_x^2)\zeta \underset{x,x^{}}{}w_{x,x^{}}\psi _x^1\overline{\psi }_x^1\overline{\psi }_x^{}^2\psi _x^{}^2.$$

The advantage of this representation is that a quadratic form of Grassmannians in the exponent can be integrated out, leading to a determinant. Of course, in the case under consideration $`S_{ferm}`$ is not a quadratic form. However, the non-quadratic part can be expressed by a coupling of a quadratic form to an Ising spin $`S_{x,x^{}}=\pm 1`$:

$$\mathrm{exp}\left(\zeta \underset{x,x^{}}{}w_{x,x^{}}\psi _x^1\overline{\psi }_x^1\overline{\psi }_x^{}^2\psi _x^{}^2\right)$$

$$=2^{𝒩\beta }\underset{\{S_{x,x^{}}=\pm 1\}}{}\mathrm{exp}\left(\underset{x,x^{}}{}\sqrt{\zeta w_{x,x^{}}}S_{x,x^{}}(\psi _x^1\psi _x^{}^2+\overline{\psi }_x^1\overline{\psi }_x^{}^2)\right),$$

where $`𝒩`$ is the number of lattice sites. Now the Grassmann integration can be performed and gives

$$Z_{HC}=2^{𝒩\beta }\underset{\{S_{x,x^{}}=\pm 1\}}{}det(\mathrm{𝟏}+𝐔)^2$$

with the matrix elements

$$U_{x,x^{}}=\sqrt{\zeta w_{x,x^{}}}S_{x,x^{}}.$$

The density of bosons $`n`$ is the expectation value

$$n=\frac{1}{2𝒩\beta }\frac{\underset{\{S_{x,x^{}}\}}{}det(\mathrm{𝟏}+𝐔)^2Tr((\mathrm{𝟏}+𝐔)^1𝐔)}{_{\{S_{x,x^{}}\}}det(\mathrm{𝟏}+𝐔)^2}$$

$$\frac{1}{2𝒩\beta }Tr((\mathrm{𝟏}+𝐔)^1𝐔)_S,$$

with respect to the distribution of Ising spins

$$P(𝐔)=\frac{det(\mathrm{𝟏}+𝐔)^2}{_{\{S_{x,x^{}}\}}det(\mathrm{𝟏}+𝐔)^2}.$$

This representation is a good starting point for a Monte Carlo simulation of the hard-core Bose gas .

In one dimension the partition function can be evaluated exactly, since periodic boundary conditions are satisfied for each random walk individually (i.e. the winding number is zero) due to the hard-core condition. Then the determinant is just a sum over all products along possible random walks with periodic boundary conditions in $`t`$:

$$det(\mathrm{𝟏}+𝐔)=1+\underset{\pi id}{}\underset{r}{}g_{r,\pi }$$

with

$$g_{r,\pi }=\underset{t}{}U_{x,\pi (x)}.$$

Each $`g_{r,\pi }`$ contains a product of Ising spins. Now one takes the square of the determinant and the sum with respect to the configurations of the Ising spins. This gives

$$\underset{\{S_{x,x^{}}=\pm 1\}}{}(1+\underset{\pi id}{}\underset{r}{}g_{r,\pi })^2=2^{𝒩\beta }+\underset{\pi id}{}\underset{\{S_{x,x^{}}=\pm 1\}}{}\underset{r}{}g_{r,\pi }^2=2^{𝒩\beta }(1+\underset{\pi id}{}\underset{r}{}g_{r,\pi }^2),$$

since $`g_{r,\pi }^2`$ does not depend on the Ising spins:

$$g_{r,\pi }^2=\underset{t}{}U_{x,\pi (x)}^2=\zeta ^\beta \underset{t}{}w_{x,\pi (x)}.$$

Therefore, the partition function can also be written as a determinat

$$Z_{HC}=det(\mathrm{𝟏}+\zeta 𝐰).$$

For the occupation of the quantum number $`l`$ follows

$$N_l=\frac{\zeta ^\beta e^{\beta ϵ_l}}{1+\zeta ^\beta e^{\beta ϵ_l}}\{\begin{array}{cc}1\hfill & \text{for }\zeta e^{ϵ_l}>1\hfill \\ 0\hfill & \text{for }\zeta e^{ϵ_l}<1\hfill \end{array}$$

for $`\beta \mathrm{}`$. Thus a quantum state is either empty or singly occupied, a consequence of the hard-core interaction. This is also the property of a non-interacting fermion gas. It reflects the equivalence of the hard-core Bose gas and the free fermion gas in one dimension .

## 4 Generalization to $`N`$ States and $`N`$ Colors

Generalizations of the hard-core Bose gas have been proposed that lead to exactly soluble limits: the $`N`$ state bosons (NSB)  and the $`N`$ color bosons (NCB) . In the first case a boson can occupy one of $`N`$ different states at each lattice site, where these states have the same energy. The hard-core interaction applies only to bosons in the same state. Therefore, bosons can avoid the interaction during the tunneling process by choosing an unoccupied state. If $`N`$ is large the NSB are effectively weakly interacting. Nevertheless, it turns out that the interaction plays an important role even in the limit $`N\mathrm{}`$ . The other generalization uses colored bosons: the hard-core interaction applies only to bosons with the same color. In this situation the interaction is effectively strong for all values of $`N`$. In particular, it will be shown that even in the limit $`N\mathrm{}`$ it can destroy the Bose-Einstein condensate and create a Mott insulator, an effect not possible in the limit $`N\mathrm{}`$ of the NSB.

### 4.1 Effective Field Theory of NSB

The extension of $`S_{ferm}`$ to the NSB is given by

$$S_{NSB}=\underset{x}{}\underset{\alpha =1}{\overset{N}{}}(\psi _x^{1,\alpha }\psi _x^{2,\alpha }+\overline{\psi }_x^{1,\alpha }\overline{\psi }_x^{2,\alpha })+\frac{\zeta }{N}\underset{x,x^{}}{}w_{x,x^{}}\underset{\alpha ,\beta =1}{\overset{N}{}}\psi _x^{1,\alpha }\overline{\psi }_x^{1,\alpha }\overline{\psi }_x^{}^{2,\beta }\psi _x^{}^{2,\beta }.$$

The Hubbard-Stratonovich transformation

$$\underset{\alpha }{}\psi _x^{1,\alpha }\psi _x^{2,\alpha }\phi _x,\chi _x$$

leads to the partition function

$$Z_{NSB}=e^{NS_{NSB}^{}}𝒟[\phi ,\chi ]$$

with the action

$$S_{NSB}^{}=(\phi ,(\mathrm{𝟏}+𝐰)^1\phi ^{})+(\chi ,\chi ^{})\underset{x}{}\mathrm{log}[\zeta ^1+(\phi _x+i\chi _x)(\phi _x^{}+i\chi _x^{})].$$
(9)
$`S_{NSB}^{}`$ does not depend on $`N`$ and has a $`U(1)`$ symmetry: it is invariant under the transformation $`\phi _xe^{i\eta }\phi _x`$, $`\chi _xe^{i\eta }\chi _x`$, with a global phase $`\eta `$.

### 4.2 Effective Field Theory of NCB

In the case of NCB the action $`S_{ferm}`$ is extended to

$$S_{NCB}=\underset{x}{}\underset{\alpha =1}{\overset{N}{}}(\psi _x^{1,\alpha }\psi _x^{2,\alpha }+\overline{\psi }_x^{1,\alpha }\overline{\psi }_x^{2,\alpha })+\frac{\zeta }{N}\underset{x,x^{}}{}w_{x,x^{}}\underset{\alpha ,\beta =1}{\overset{N}{}}\psi _x^{1,\alpha }\overline{\psi }_x^{1,\beta }\overline{\psi }_x^{}^{2,\beta }\psi _x^{}^{2,\alpha }.$$

The Hubbard-Stratonovich transformation

$$\underset{\alpha }{}(\psi _x^{1,\alpha }\psi _x^{}^{2,\alpha }+\overline{\psi }_x^{1,\alpha }\overline{\psi }_x^{}^{2,\alpha })u_{x,x^{}}$$

leads to the partition function

$$Z_{NCB}=\mathrm{exp}\left\{N\underset{x,x^{}}{}\frac{(u_{x,x^{}})^2}{\zeta w_{x,x^{}}}+2N\mathrm{log}det(\mathrm{𝟏}+𝐮)\right\}𝒟[𝐮].$$
(10)
With the choice of $`𝐰`$ given in Eq. (6) the action $`S`$ can also be expressed, after a decomposition of the matrix field $`𝐮=𝐮_{}+𝐮_{}`$ ($`𝐮_{}`$ is diagonal in space), as

$$S_{NCB}^{}=\frac{1}{\zeta (1J)}Tr(𝐮_{}𝐮_{}^T)+\frac{6}{\zeta J}Tr(𝐮_{}𝐮_{}^T)2\mathrm{log}det(\mathrm{𝟏}+𝐮_{}+𝐮_{})$$

This is invariant under the orthogonal transformation

$$𝐮\mathrm{𝐎𝐮𝐎}^T.$$

## 5 $`N\mathrm{}`$ Limits

The new representations of the partition functions depend on $`N`$ only through a prefactor in front of the $`N`$-independent action. Therefore, the functional integration can be performed for large values of $`N`$ by the saddle-point (SP) approximation. That means we have to solve the extremal condition of the action $`S`$: $`\delta S=0`$. Systematic corrections to the SP solution are obtained in a $`1/N`$ expansion.

### 5.1 Results for the NSB

There are two SP solutions, one is the trivial solution $`\phi =\chi =0`$ which is valid for $`\zeta <1`$ and a non-trivial solution with $`|\chi |=|\phi |/2=\sqrt{11/\zeta }`$ which is valid for $`\zeta >1`$. This implies for total density of bosons at zero temperature

$$n(\zeta )\{\begin{array}{cc}0\hfill & \text{for }\zeta <1\hfill \\ 1\zeta ^1\hfill & \text{for }\zeta >1\hfill \end{array}.$$

Moreover, all bosons are in the condensate in the $`N\mathrm{}`$ limit. However, the $`1/N`$ corrections indicate a depletion of the condensate, as shown in Fig. 3. This phenomenon reflects a partial destruction of the condensate due to interaction. It can be understood as a precursor for a new groundstate in the interacting Bose gas, visible in the NCB as the Mott insulator (s. next section). The critical point $`\zeta =1`$ of the Bose-Einstein condensation does not depend on the tunneling rate $`J`$ and agrees with $`J=0`$ for all $`N`$. It should be noticed that the limit $`N\mathrm{}`$ of the NSB leads to the Gross-Pitaevskii equation , at least in a close vicinity of the critical point (dilute regime).

Quasiparticles are available from $`1/N`$ corrections. In the condensate (i.e. $`\zeta 1`$) they have a dispersion on large scales with respect to the three-dimensional wavevector $`k`$ as

$$ϵ(k)=\sqrt{\frac{J(\zeta 1)}{3}k^2+(J/6)^2k^4}.$$

Due to the spontaneously broken $`U(1)`$ symmetry of the NSB there is a Goldstone mode $`\varphi `$ for $`\zeta >1`$ which satisfies the equation

$$\left(c^2k^2\omega ^2\right)\varphi (k,\omega )=0,$$

where $`c=\sqrt{J(\zeta 1)/3}`$ is the sound velocity.

### 5.2 Results for the NCB

For the SP equation of the action $`S_{NCB}`$ one can distinguish symmetric and symmetry-breaking solutions. Examples of the former are the trivial solution $`𝐮_{}=𝐮_{}=0`$ (empty phase) and $`𝐮_{}0`$, $`𝐮_{}=0`$ with identical matrix elements

$$u_{x,x^{}}=\{\begin{array}{cc}u_{}\hfill & \text{for }r^{}=r\text{}t^{}=t+1\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}.$$

This ansatz leads to the mean-field action

$$\overline{S}_{NCB}=𝒩\beta \left[\frac{u_{}^2}{\zeta (1J)}\frac{2}{\beta }\mathrm{log}(1+u_{}^\beta )\right]$$

which is minimized by the solution of the SP equation

$$u_{}^2=\zeta (1J)\frac{u_{}^\beta }{1+u_{}^\beta }.$$

For zero temperature the value of $`u_{}`$ jumps from 0 to 1 at $`\zeta =1/(1J)`$. The density of bosons can be calculated with Eq. (8) from the mean-field action and gives $`n=\mathrm{\Theta }(\zeta (1J)1)`$. This is the behavior of the system without tunneling. The jump should be intercepted by an intermediate phase, most likely the Bose-Einstein condensate. The off-diagonal long-range order of the NCB in the large-$`N`$ limit is related to a symmetry-breaking SP solution with $`𝐮_{}0`$. This represents a Bose-Einstein condensate with a Goldstone mode. It is anticipated that this solution agrees qualitatively with the result of the NSB, discussed in the previous section. The resulting phase diagram, obtained from the $`N\mathrm{}`$ limits of the NSB and the NCB, is shown in Fig. 4.

## 6 Conclusions

The hard-core Bose gas at zero temperature has been studied in terms of the soluble $`N\mathrm{}`$ limits of NSB and NCB. NCB describe three different phases of the hard-core Bose gas, characterized by a matrix-like order parameter $`𝐮=𝐮_{}+𝐮_{}`$. There is an empty phase ($`𝐮_{}=𝐮_{}=0`$), a Bose-Einstein condensate ($`𝐮_{}=0,𝐮_{}0`$), and a Mott insulator ($`𝐮_{}0,𝐮_{}=0`$). In contrast to the NCB the complex scalar order parameter of the NSB describes only an empty phase and a Bose-Einstein condensate but not the Mott insulator. Thus, the NSB has only one symmetric SP solution (empty phase), whereas the NCB has at least two (empty phase and Mott insulator). However, the depletion of the condensate in $`o(1/N)`$ of the NSB indicates the partial destruction of the latter. The main result of this study is the effective action of the NCB of the matrix field $`𝐮`$

$$S_{NCB}^{}=\frac{1}{\zeta (1J)}Tr(𝐮_{}𝐮_{}^T)+\frac{6}{\zeta J}Tr(𝐮_{}𝐮_{}^T)2\mathrm{log}det(\mathrm{𝟏}+𝐮_{}+𝐮_{}),$$

where its minima represent different phases of the strongly interacting Bose gas.

Acknowledgement

The author is grateful to Ch. Moseley for useful discussions.


text:
# Diluted Magnetic Semiconductor at Finite Temperature

## Abstract

We studied the diluted magnetic semiconductor by the self-consistent Green’s function approach, which treats the spin-wave kinematics appropriately at finite temperatures. Our approach leads to a simple formula for the critical temperature in a wide range of parameter space. In addition, the magnetization curve versus temperature in some regimes is concave, which is dramatically different from the usual convex shape. Finally, we discuss the possibility of generalizing the current theory to include the realistic band structure, electronic correlations and disorders in a systematic way.

Diluted magnetic semiconductor (DMS) has attracted intense attentionsOhno98a ; Prinz98  for its potential applications in spintronics devices. Robust ferromagnetic order has been observed in (Ga1-xMnx)As up to 110 KOhno96 ; Ohno99 . Magnetically doped wide bandgap semiconductors and oxides such as GaN, ZnO, TiO2 even exhibit ferromagnetism at room temperatureMatsumoto01 ; Theodoropoulou01 ; Lee02 , although the magnetization is less robust than the doped III-V semiconductors. While lots of efforts are focused on the search of optimal materials with enhanced critical temperaturesDatta90 , it remains a challenging task to describe the coupled localized moments and the itinerant carriers in an appropriate way.

Part of the difficulty lies in the fact that the spatial fluctuations of the ferromagnetic order is large at finite temperature and the usual Weiss mean-field theory does not work. To account for the spatial fluctuations, numerical approachesDietl00 ; Schliemann01 ; Sham01 ; Sun02 , such as local density functional approximation and Monte Carlo simulations, are quite helpful in estimation of various thermodynamical properties. However, it is rather difficult to study the electronic transport by these numerical approaches. Analytical approachesAkai98 ; Konig00 ; Litvinov01 ; Berciu01  provide partial descriptions in several particular limits by treating the impurity spin semiclassically or ignoring the interactions and the kinematic constraints among spin waves. At finite temperatures, the average spin-wave density is large so that these approximations are no longer appropriate. It is therefore desirable to develop a spin-wave theory which works at finite temperature.

In this Letter, we adapt the self-consistent Green’s function approach to describe the fluctuating spin correlations at finite temperature. Since the kinematic constraint of spin waves are treated exactly in the equation of motion, this method can be applied to a wide range of parameter space, even very close to the critical temperatureVogt01 ; Yang01 . Indeed, the critical temperature is determined by the simple formula,

$$k_BT_c=\frac{S+1}{3}\left[\frac{1}{cV}\underset{p}{}\frac{\alpha _c}{\mathrm{\Omega }_c(p)}\right]^1,$$
(1)
where $`c`$ is the density of impurity spins. The spin-wave dispersion $`\mathrm{\Omega }_c(p)`$ and the impurity spin polarization $`\alpha _c`$ are both determined self-consistently near the critical temperature $`TT_c`$. The trends of the critical temperature upon the change of parameters are studied in detail later.

We model the DMS by the Hamiltonian, containing only the kinetic energy of itinerant carriers and the exchange interaction between the itinerant and the localized impurity spins,

$$H=H_0+Jd^3r𝑺(r)𝒔(r),$$
(2)
where $`J>0`$ is the strength of exchange interaction. The impurity spin density is $`𝑺(r)=_I\delta ^{(3)}(rR_I)𝑺_I`$ while the itinerant spin density is $`𝒔(r)=\psi ^{}(r)(𝝈/2)\psi (r)`$. The band structure of the itinerant carriers is described by $`H_0`$, which depends on the host semiconductors. Since our emphasis here is how to cope with spatially fluctuations appropriately, the dispersion is taken as the simplest parabolic band, $`H_0=p^2/2m^{}`$. Generalization to more realistic but complex band structures, such as the Luttinger model, can be achieved straightforwardly.

Since the impurity spins are randomly doped into the host semiconductor, their positions are random. If the disorder is strong, the itinerant electrons are localized and the percolation approachLitvinov01 ; Berciu01  would be more appropriate. However, we are interested in the metallic regime where itinerant carriers are delocalized. The disorder also plays a crucial role in smoothing out the impurity spin density $`𝑺(r)`$. For instance, the magnitude of the spin density is smeared after coarse-graining,

$$𝑺^2(r)_{R_I}c^2S(S+1),$$
(3)
where $`c`$ is the impurity spin density and the average is taken over random locations of the impurity spins. Therefore, the presence of weak disorder allows a field-theory description in the continuous limit.

The dynamics of the impurity spins is described by the thermal Green’s function,

$`D(r,\tau )`$ $``$ $`S^+(r,\tau );S^{}(0,0)`$ (4)
$``$ $`\mathrm{\Theta }(\tau )[S^+(r,\tau ),S^{}(0,0)].`$
The equation of motion for $`D(r,\tau )`$ would involve more Green’s functions of higher orders. The exact solution then involves the Green’s functions of all orders and is not feasible in general. However, within mean-field approximation, the higher-order Green’s functions can be decomposed into simpler ones and eventually a self-consistent solution is possible. For the spin-wave propagator $`D(q,i\nu _n)`$ in momentum space, the mean-field decomposition simplifies the equation of motion,

$$i\nu _nD(q,i\nu _n)=1+Js_z+\frac{JS_z}{c^{}V}\underset{k}{}F(k,k+q,i\nu _n),$$
(5)
where $`c^{}`$ is the density of itinerant carriers. Notice that it only involves one additional Green’s function $`F(k,k+q,i\nu _n)\psi _{}^{}(k)\psi _{}(k+q);S^{}(0,0)`$. Applying the same trick again, the equation of motion for $`F(k,k+q,i\nu _n)`$ is

$$F(k,k+q,i\nu _n)=\frac{Jc^{}}{2}\frac{f_{}(k)f_{}(k+p)}{i\nu _n+\mathrm{\Delta }+ϵ_kϵ_{k+p}}D(q,i\nu _n),$$
(6)
where $`f_,(k)=[e^{\beta (ϵ_k\pm \mathrm{\Delta }/2\mu )}+1]^1`$ is the Fermi distribution for itinerant carriers with different spins. The Zeeman gap $`\mathrm{\Delta }JS_z`$ in the electronic band structure is due to the ferromagnetic order of the impurity spin and has to be determined self-consistently.

From Eqs. (5) and (6), we can solve for the spin-wave propagator. Upon the Wick rotation, $`i\nu _n\mathrm{\Omega }+i\eta `$, the spin-wave dispersion is identified as the simple pole in the Greens’ function,

$$\mathrm{\Omega }Js_z\frac{J\mathrm{\Delta }}{2V}\underset{k}{}\frac{f_{}(ϵ_k)f_{}(ϵ_{k+p})}{\mathrm{\Omega }+\mathrm{\Delta }+ϵ_kϵ_{k+p}+i\eta }=0.$$
(7)
The average spin densities $`s_z`$ and $`S_z`$ in the above equation remain unknown and need to be determined self-consistently.

The average itinerant spin density $`s_z`$ is just the difference of the spin densities in majority and minority bands, split by the Zeeman gap $`\mathrm{\Delta }`$. The relation between the average localized spin density $`S_z`$ and the spin-wave dispersion $`\mathrm{\Omega }(k)`$ is more subtle due to the non-trivial spin kinematic constraint. From Callen’s formula,Callen63

$$\frac{1}{c}S_z=Sn_{sw}+\frac{(2S+1)n_{sw}^{2S+1}}{(1+n_{sw})^{2S+1}n_{sw}^{2S+1}}.$$
(8)
Here $`n_{sw}=(1/cV)_k[e^{\beta \mathrm{\Omega }(k)}1]^1`$ is the average number of spin waves over all momenta. At low temperatures, the average spin-wave number is small and the last term in Eq. (8) can be safely ignored. The magnetization show the $`T^{3/2}`$ behavior as in the independent spin-wave approximation. As the temperature approaches the critical regime, the spin-wave density becomes large and the kinematic constraint becomes important. Solving Eqs. (7) and (8), both the spin-wave dispersion $`\mathrm{\Omega }(k)`$ and the averaged impurity spin density $`S_z`$ are obtained self-consistently.

The dispersion from Eq. (7) has two branches due to the presence of both itinerant and localized spins. Because the gapless spin-wave fluctuations dominates, the optical branch can be safely ignored. Besides, the spectral weight of the optical mode is also small due to the dilute density of itinerant carriers. In Fig. 1, the lower branch of spin-wave dispersion at zero temperature is shown, corresponding to the Goldstone excitation of the ferromagnetic order. Our numerical results show that the imaginary part of the spin-wave self energy, due to the presence of Stöner continuum (the optical branch), is negligible in all regimes and the dispersion is sharp. Thus, it is a reasonable approximation to assume the spectral weight is carried by the gapless spin-wave excitations only.

The magnetization curves at different spin densities $`c^{},c`$ are plotted in Figs. 3 and 3. Throughout this Letter, the exchange coupling and the effective mass are fixed at typical values $`J=0.15`$ eV nm3 and $`m^{}=0.5m_e`$Ohno99 . First of all, we study the temperature dependence of the magnetization at different impurity spin densities $`c`$, while the ratio of itinerant and localized spin densities fixed at $`c^{}/c=0.1`$. It is quite interesting to notice that the magnetization drops more dramatically near the critical regime when the impurity spin density is large. As shown in Fig. 3, the critical temperature, computed by Eq. 1, increases monotonically as the densities increaseSun03 . This monotonic increase in the critical temperature is qualitatively (not quantitatively) the same as the Weiss mean-field theory and also agree with the experiments.

On the other hand, if the impurity spin density is hold constant, say $`c=1`$ nm-3, both the magnetization curve and the trend of the Curie temperatures show interesting behaviors, which deviate from the Weiss mean-field theory. Starting from extremely dilute density of the itinerant carriers, the magnetization curve is concave as shown in Fig. 3, in contrast to the usual convex curvature in Weiss mean-field theory. This concave shape resembles the magnetization curve in percolation theory near the boundary of the metal-insulator transition and is often used to be an indication of localization in the presence of disorderLitvinov01 ; Berciu01 . However, our results show that a concave magnetization curve is not necessarily tied up to the localization tendency. We emphasize that the disorder of impurity spin is included in our approach only through the coarse-graining procedure and the transport of the itinerant carriers is assumed ballistic here. So it is rather surprising, by including the spin-wave fluctuations appropriately at finite temperatures, the magnetization curve is concave at dilute densities.

In the regime where the magnetization curve is concave, the critical temperature increases with the carrier concentration $`c^{}`$. Beyond an optimal density, $`c^{}0.01`$ nm-3 (the impurity spin density is fixed at $`c=1`$ nm-3), the critical temperature reaches the maximum and starts to fall back as in the inset of Fig. 3. The suppression of the critical temperature is mainly due to the oscillatory RKKY interaction at higher concentrations. The effective coupling between impurity spins is frustrated and not longer purely ferromagnetic. The existence of an optimal density is totally missed in the Weiss mean-field theory, where both the spatially varying fluctuations and the quantum frustrations are ignored.

Beyond the optimal density, not only the critical temperature falls back, the shape of the magnetization also undergoes an interesting change. There are two significant features. One is the magnetization curve becomes convex again. Another feature is that the magnetization has a very sharp decrease near the critical temperature, making it look almost like the first-order phase transition. To make sure that the sharp drop is not an artifact of numerical errors, the data points near this regime are chosen very closely to ensure we still have enough points in the steep regime. In particular, for $`c^{}=0.05`$, it is spectacular that the magnetization drops 50% within 1K. It is not clear at this point what is the underlying mechanism behind this dramatic suppression of magnetization. A more sophisticated theory starting from the critical point might be able to address this interesting dive of magnetization near the critical temperature.

Since the global trend of the critical temperature is important, it is plotted in Fig. 4. Comparing with Fig. 2 in Ref.Dietl00 , the self-consistent Green’s function approach produces a more complicated landscape. In the Weiss mean-filed theory, the critical temperature is

$$T_c^{\mathrm{MF}}=\frac{\chi _P}{(g^{}\mu _B/2)^2}\frac{S(S+1)J^2N}{12},$$
(9)
where $`g^{}`$ is the $`g`$-factor of the carriers and $`\chi _P`$ is their Pauli susceptibility, which is proportional to the effective band mass. The above formula would produce a profile with monotonically increasing $`T_c`$ as the densities become larger. It is clear from Fig. 4 that the profile of the Curie temperature has a ridge, roughly along the curve $`c100c^{}`$ with the given parameters in this Letter, and is qualitatively different from previous studies. A closer check would find that our approach also produces quantitative differences. Comparing with the estimates of Curie temperatures, previously studied by one of the authors in Refs. Schliemann01  and Konig00 , the present approach gives a lower $`T_c`$ in all densities regimes. So it seems that the inclusion of spin kinematics at finite temperatures is not only important to get the right trend of the critical temperatures, but also crucial in estimating their values.

Finally, we address the important aspects of physics which are left out so far. To make quantitative comparison with the experiments, it is crucial to adapt the realistic band structure of the host semiconductors, for instance, the six-band Luttinger model for GaAs. The inclusion of the more complex band structure would not cause formidable messes in calculating the Green’s function self-consistently. While most of the conclusions drawn from the simple parabolic band should remain valid, a more realistic band structure is desirable for making quantitative predictions. The electronic correlation can also be included at the mean-field level in a systematic way. Since Coulomb repulsive interaction stabilizes the ferromagnetic phase, the Curie temperature is expected to be higher. To include the disorder is more subtle. In the diffusive regime, one can replace the ballistic propagator of the itinerant carriers by the diffusive one to account for the disorder effects.

In conclusion, we employ the self-consistent Green’s function approach to study the DMS and derive a general formula for the Curie temperature. In addition, we demonstrate the interesting crossover of the magnetization curve from concave to convex, which is not driven by the localization effects.

We thanks Chung-Yu Mou and Ming-Fong Yang for fruitful discussions and the support of National Science Council in Taiwan. The hospitality of the National Center for Theoretical Science, where the work was initiated, is greatly acknowledged.


text:
# Suppression of 2D superconductivity by the magnetic field: quantum corrections vs superconductor-insulator transition.

## Abstract

Magnetotransport of superconducting Nd2-xCexCuO4+y (NdCeCuO) films is studied in the temperature interval 0.3–30 K. The microscopic theory of the quantum corrections to conductivity, both in the Cooper and in the diffusion channels, qualitatively describes the main features of the experiment including the negative magnetoresistance in the high field limit. Comparison with the model of the field-induced superconductor–insulator transition is included and a crossover between these two theoretical approaches is discussed.

The superconductor–insulator transition (SIT) is an example of the quantum phase transitions Sondhi  which constitutes drastic change of the ground state of the system at zero temperature with varying a parameter. The field was pioneered by A. Goldman et al. in the 1989 Bi  who obtained the transition from insulating to superconductive state in the thin Bi film with the change of its thickness. Later, Fisher Fisher  suggested existence of magnetic field-induced SIT in two-dimensional (2D) systems and Hebard and Paalanen demonstrated Hebard ; Paal  such a transition in amorphous InOx films. Numerous results obtained in several other materials by different groups Kapitulnik ; Okuma ; Tanda ; NCCO  were also interpreted within the framework of the field-induced SIT. Main arguments in favor of this interpretation were negative derivative of resistance $`R/T`$ in the fields above the critical and existence of a finite-size scaling, i.e. existence of some critical region on the ($`T,B`$)–plane where the behavior of the system was governed by competition of the quantum phase transition correlation length $`\xi (BB_c)^\nu `$ and thermal length $`L_TT^{1/z}`$ with $`z`$ and $`\nu `$ being constants called the critical exponents. All relevant quantities in this region are supposed to be universal functions $`f`$ of ratio of the lengths which can be written in the form of scaling variable $`(BB_c)/T^{1/z\nu }`$. For the resistivity in two dimensions $`R_{\mathrm{}}`$ this dependence takes form Fisher

$$R_{\mathrm{}}(B,T)=R_cf[(BB_c)/T^{1/z\nu }],$$
(1)
where $`R_c`$ is a constant of the order of $`h/4e^26.5`$ k$`\mathrm{\Omega }`$. It is called the critical resistance.

In the analysis of the experiments Kapitulnik ; Okuma ; Tanda ; NCCO  the negative derivative $`R/T`$ was rated as an indicator of the insulating state. However, that is not enough: the characteristic of any insulator is the exponential temperature dependence of the resistance. This was demonstrated only in InOx films ourInO . The growth of the resistance with decreasing temperature on the non-superconducting side of the field-induced transition in the experiments with MoGe Kapitulnik , MoSi Okuma  and NdCeCuO Tanda ; NCCO  was minuscule, about ten percent at its best. It reminded more a metal with quantum corrections to its conductivity than an insulator. Usually, the authors do not dwell on the issue, considering weak localization-like behavior to be telltale sign of insulator — as, according to scaling hypothesis 4band , there is no non-superconducting delocalized state at zero temperature in 2D and weak localization is expected to transform sooner or later into strong. However, this crossover might be postponed to extremely low temperature which would never be achieved in practice.

There exists one more sign of SIT. According to the boson–vortex duality model Sondhi ; Fisher , the insulating state which appears as the result of SIT is rather specific; it contains pair correlations between the localized electrons as the remnant of the superconducting pairing. Such insulator is called the Bose-insulator Paal  and the correlated electrons are called localized electron pairs. These correlations should be destroyed by strong magnetic field leading to increase of the carrier mobility, to the negative magnetoresistance Shep1  and even to a reentrant insulator–normal-metal transition ourInO . The negative magnetoresistance was observed in MoSi OkumaNMR  and NdCeCuO NCCO . But it was much weaker than in InO, just the same as the growth of the resistance with decreasing temperature discussed above.

When comparing the whole set of data in InO Hebard ; Paal ; ourInO ; ourInO2  with those in MoGe Kapitulnik , MoSi Okuma  and NdCeCuO NCCO , one can’t help impression that they have many similar features though of different scales of magnitude. At the same time, it was shown in a set of InOx films with various oxygen content $`x`$ that in low-resistivity films a transition to the metallic state substitutes SIT, the rate of the temperature dependence scales down and the whole pattern of curves approaches that of the usual superconducting transition ourInO ; ourInO2 . The main idea of this paper follows from this observation. It is to compare experimental set of data of a “small-scale” type with the theory of the superconducting transition in dirty limit and, keeping in mind its features related to SIT, to build a bridge between SIT and thermodynamic superconducting transition.

Experiment was performed on the 1000Å-thick films of Nd2-xCexCuO4+y (NdCeCuO) obtained by laser ablation with CuO2 planes parallel to the plane of the film. Films were not superconductive as-grown. In order to obtain superconductivity they were annealed at 720C in the flowing 4He gas for several hours. As we aimed to study vicinity of the SIT, we were not trying to reach maximal $`T_c`$ of this material, but were paying attention for smoothness and width of the zero-field transition. A sample was chosen with zero-field transition temperature $`T_{c0}=11.8\pm 0.4`$ K (found by fitting of the superconducting fluctuation contribution to the conductivity above $`T_{c0}`$) and the transition width $`\mathrm{\Delta }T2`$ K.

The resistivity was measured in the ab-plane by the 4-terminal technique. Both current and potential probes were attached on the surface of the films by silver paste. Distance between potential probes corresponded to one square. Magnetic field was applied perpendicular to the film plane (along c-axis). Data, both as a function of field at constant temperature and as a function of temperature at constant field, were obtained though only the latter will be presented below. The upper panel of Fig.1 presents an overview of the impact of the field on $`R(T)`$ dependence and the lower one zooms in on the region of interest, i.e. on the low temperature and high field region.

On the right axis of Fig. 1 the resistance reduced per one CuO2 plane per square is denoted. As NdCeCuO is highly anisotropic anisotropy , it is reasonable to assume the film to be a stack of 2D conducting CuO2 planes with interplane spacing 6Å, quasi-independent and connected in parallel. This is supported by observations of 2D character of quantum interference corrections Ponomarev2D  and magnetoresistance Kussmaul2D . Later on we continue discussion in terms of this variable, disregarding full resistivity and actual thickness of the film. As one could see from Fig. 1, the value of the resistance per layer stays quite far from the quantum resistance $`h/4e^2`$ expected for the SIT.

The data are quite typical for the material (cf., for example, Ichikawa et al. NCCO ). In the low field region, the transition is shifted to the lower temperature as the field increases while the shape of the transition is preserved relatively well. Above 2 T, the transition broadens drastically and eventually disappears; at about 3.5 T, the $`dR/dT`$ changes its sign. At higher fields, above 5 T, the resistance starts to decrease with the increasing field; it follows from the crossing of the 5 T and 7 T curves that a region of the negative magnetoresistance exists below 0.8 K and at $`B>5`$ T.

The set of curves $`R(T)`$ on the lower panel of Fig. 1 is similar to those obtained in Kapitulnik ; Okuma ; Tanda ; NCCO  which had been regarded as a field-induced SIT. Low-field curves (which bend down) may be supposed to reach zero resistivity at zero temperature and to become superconductor, high field curves (which bend up) may be supposed to diverge toward zero temperature and become insulator. In between, there is a curve which is almost horizontal; it manifests itself as common crossing point of all isotherms on the $`RB`$ graph. The corresponding state should be considered as the critical one with the temperature-independent resistance at the critical field $`B_c3.5`$ T. But instead of seeking scaling parameters we shall compare experimental data with the microscopic theory of the superconducting transition in dirty limit formulated in terms of quantum corrections to the classical Drude conductivity $`\sigma _0=e^2/h(k_Fl)`$ where $`k_F`$ is the Fermi wavevector and $`l`$ is the elastic mean free path. This comparison became executable due to recent progress in calculation of the corrections due to superconducting fluctuations Galitsky .

All quantum corrections fall into two categories — a one-particle correction, usually called weak localization, and those due to e-e interactions. The latter are divided into a diffusion channel correction (also known as Aronov–Altshuler term) and Cooper channel corrections (also known as superconductive fluctuations corrections which include Aslamazov–Larkin, Maki–Thompson and DOS terms). Weak localization and Aronov–Altshuler corrections diverge at $`T0`$, Cooper channel corrections diverge at $`TT_c(B)`$ with $`T_c(B)`$ being mean field transition temperature. When the superconductivity is suppressed by the magnetic field, $`T_c(B)0`$ and all corrections are important.

Recently, Galitski and Larkin Galitsky  succeeded in extending calculations in the Cooper channel for two-dimensional superconductors to the low temperature $`TT_c(0)`$ and high magnetic field $`BB_{c2}(0)`$. The correction to the conductivity in the dirty limit $`\delta \sigma `$ is obtained as the sum of contributions of ten Feynman diagrams in the first (one-loop) approximation and can be written in the form

$$\delta \sigma =\frac{4e^2}{3\pi h}\left[\mathrm{ln}\frac{r}{b}\frac{3}{2r}+\psi (r)+4(r\psi ^{}(r)1)\right],$$
(2)
where $`r=(1/2\gamma ^{})(b/t)`$, $`\gamma ^{}=e^\gamma =1.781`$ is the exponential of Euler’s constant, and $`t=T/T_{c0}1`$ and $`b=(BB_{c2}(T))/B_{c2}(0)1`$ are reduced temperature and magnetic field.

To compare these calculations with the experiment, we added to the correction (2) an additional term to account for Aronov-Altshuler contribution, which is assumed to be field independent. Weak localization was omitted because we are interested in the region of rather strong magnetic fields where this correction was expected to vanish. Finally, we arrived at the formula

$$R_{_{\mathrm{}}}^1(B,T)=\sigma _0+\delta \sigma (B,T)\alpha \frac{e^2}{h}\mathrm{ln}(T/T^{}).$$
(3)
Inserting $`T_{c0}=11.8`$ K and experimental value of the classical conductivity $`\sigma _0=1/R_{_{\mathrm{}}}`$(7 T, 20 K), and choosing $`T^{}=20`$ K to make the last term zero at 20 K and $`\alpha =1/2`$ to match the temperature dependence of the experimental curve at 7 T, we get the plot of Fig. 2 which can be compared with the experimental one, Fig. 1b. (Note, that on Fig. 2 curves are labelled by reduced field values, those in units of $`B_{c2}(0)`$. The same can’t be done on Fig. 1, because experimental value of $`B_{c2}(0)`$ is a bit uncertain.)

As one can see, the picture bears clear resemblance to the experiment — there is separation between low-field curves which “bend down”, and high-field which “bend up”; there is also high field negative magnetoresistance at low temperature. There are two remarkable points: i) the scales of variation of resistance both with temperature and magnetic field are correct; ii) the region and the magnitude of the negative magnetoresistance are in reasonable agreement with the experiment as well.

However, the similarity is qualitative. It is difficult to make it quantitative and both the experiment and the theory are responsible for this.

Disadvantage of the experiment is hidden in macro-inhomogeneity of the film. It follows from Fig. 2 that small 2%–3% changes of $`B_{c2}(0)`$ lead to drastic shift in the shape of curves $`R(T)`$, especially near the critical value of $`B`$. Inevitable dispersion of the values of $`B_{c2}(0)`$ along the film smoothes the curves and clears away the extremum. Hence, one should scarcely expect to find in the experimental assortment of curves one similar to the theoretical curve labelled 1.02 (plotted by the dotted line on Fig. 2).

The expression (2) is apparently very sensitive to the function $`B_{c2}(T)`$. Basically, this function is an implicit parameter of the theory. In Ref. Galitsky, , authors used for $`B_{c2}(T)`$ the mean-field function from the Werthamer–Helfand–Hohenberg theory. It is doubtful, that this theory is applicable to high-resistive 2D objects, especially, as the shape of transition in 2D case should be affected by the vortex motion (Berezinsky–Kosterlitz–Thouless theory).

As a side note, a comment about the finite-size scaling equation (1) related to SIT. Certainly, expression (3) does not have form of equation (1) and no genuine scaling exists. However, in a restricted region of values of $`T`$ and $`B`$ representation of the theoretical curves in the form (1) can be done. This is illustrated by Fig. 3 where calculated data from the region $`0.98<B/B_{c2}(0)<1.2`$ and $`0.1<T/T_c<0.15`$ are used for the tracing. As the “critical” magnetic field $`B^{}=1.016B_{c2}(0)`$, the crossing point of several isotherms $`R(B)`$ was taken; $`B^{}`$ is the field where the minimum of the isomagnetic curve $`R(T)`$ is located in the middle of the chosen temperature region. (Actually, in the limited range of parameters $`B`$ and $`T`$ scaling always exists provided that several curves $`R(B)`$ have a common crossing point.) It follows that the scaling tracing is necessary but not sufficient element of the analysis of the SIT, especially taking into account that we always deal with the limited temperature range in the experiment.

Appearance of the negative correction to conductance in the microscopic theory of the superconductive fluctuations Galitsky  is very remarkable. It confirms that the superconducting correlations may lead at fields above the critical one not to the drop but to upsurge of the resistance. This can be regarded as the tendency toward the Bose-insulator, which could be distinguished from the Aronov–Altshuler term because it leads to the negative magnetoresistance. All the materials mentioned above can be lined up demonstrating continuous crossover from the Bose-insulator and gigantic negative magnetoresistance in InO to faint low-temperature uprise of the resistance and its tiny drop in strong magnetic fields in MoSi and NdCeCuO. In essence, these films are similar to each other: they are uniform, highly disordered films, with the resistance close to quantum value $`h/4e^2`$. Nevertheless, experimental observations on InOx and, for example, on NdCeCuO are quite different and there is a reason for it.

There is little doubt that at low enough temperature the growth of $`R(T)`$ we observe in the high magnetic field, i.e. in the normal state, will turn exponential. According to phenomenological estimate suggested by Larkin and Khmel’nitskii LKh , the crossover happens when the corrections to the conductivity reach the level of the conductivity itself. The condition $`\sigma _0(e^2/h)\mathrm{ln}T`$ gives crossover temperature

$$T_{LKh}\frac{\epsilon _F}{k_Fl}e^{2(k_Fl)},$$
(4)
where $`\epsilon _F`$ and $`k_F`$ are the Fermi energy and the Fermi wavevector and $`l`$ is the elastic mean free path LKh . Below this temperature there will definitely be a superconductive state at low field and pronounce insulating behavior at high field and there would be clear reason to apply SIT framework. So, the quantum corrections to the conductivity and the quantum phase transition phenomena are manifested at different temperature regions.

Though $`T_{LKh}`$ may be very low for normal metal ($`T_{LKh}1`$ mK already for $`k_Fl5`$), there are clear experimental indications that crossover to bosonic insulator behavior (that is, to the SIT framework) in the intermediate field range, where pair correlations are still important, occurs at higher temperature ourInO . This is consistent with theoretical observation Shep2  that the attractive interaction stimulates localization by combining single particles into pairs.

By equating two last terms in the relation (3) to the $`\sigma _0`$ and solving ensued equation one gets crossover temperature to bosonic insulator $`T_0`$ as the function of the magnetic field. These curves for $`\sigma _0`$ equal $`5e^2/h`$ (or $`k_Fl=5`$), $`7e^2/h`$ and $`9e^2/h`$ are presented on Fig.4 by solid lines. Thin solid lines present levels of $`T_{LKh}`$ determined by using only last term in the relation (3) and corresponding $`\sigma _0`$. As the equation (2) is valid only in the fields close to $`B_{c2}(0)`$, the parts of curves in the higher fields, where $`T_0(B)`$ approaches $`T_{LKh}`$, are indicated qualitatively by dotted lines. In agreement with Refs. ourInO, ; Shep2, , the crossover to activation behavior in the medium-range fields occurs at the temperatures more than order of magnitude higher than $`T_{LKh}`$. At the same time, the crossover temperature falls off exponentially with increasing classical conductivity so that for the actual value of our experiment it becomes infinitesimal. That’s why field-induced SIT is so manifest in the InOx, whereas it is not observed in MoGe or NdCeCuO, and there is no slightest sign of it in the Al film (note that according to the scaling hypothesis 4band  any metal film should become insulating at $`T=0`$ if the superconductivity is destroyed by the magnetic field).

To summarize, we compared experimental data obtained on two-dimensional NdCeCuO superconductor in magnetic field at low temperature with the calculations of quantum corrections to the conductivity and found reasonable agreement. Lack of the activation behavior at high fields (on the ”insulating side of transition”) was the main reason which made inferior comparison of the same data with the model of field-induced SIT. Apparently, this happened because the temperature range turned out to be too high for this specific material. The type of the resistance dependence on the temperature is the guide in choosing the theoretical approach. To employ framework of the SIT in its full for NdCeCuO, further essential lowering of the temperature is necessary.

This work was supported by grants RFBR-02-02-16782, RFBR-02-02-08004, RFBR-03-02-16368 and by grant of Ministry of Science.