text:
# Relics of subluminous supernovae in metal-poor stars

## 1 INTRODUCTION

Recently, Aoki et al. (2002) revealed a unique elemental abundance pattern in a carbon-rich extremely metal-poor star CS29498-043. The chemical composition of CS29498-043 has a large excess of $`\alpha `$-elements such as Mg and Si as well as C and N compared to Fe-group elements. A similar characteristic was previously identified in another C-rich star CS22949-037 (McWilliam et al., 1995; Norris, Ryan, & Beers, 2001; Depagne et al., 2002). In CS22949-037, oxygen is also found to be significantly overabundant (Depagne et al., 2002). Since a mixing process at the asymptotic giant branch of stellar evolution enhances exclusively C and N, the lack of Fe is likely to give rise to the enhancements of C, N, Mg and Si (Norris et al., 2001; Aoki et al., 2002). The lack of Fe, of course, must be initially imprinted in stars at formation because the low-mass stars considered here are unable to destroy Fe nuclei in their envelope.

Recent studies on the chemical compositions of metal-poor halo stars have claimed that these stars might have inherited the abundance pattern for the ejecta of the preceding few supernovae (SNe) (Audouze & Silk, 1995) or a single SN (Shigeyama & Tsujimoto, 1998, ST98 in the following). If this is the case, these stars would have formed from the interstellar matter (ISM) comprising the ejecta of a single SN. Depagne et al. (2002) compared the relative elemental abundances observed in CS22949-037 with those of the theoretical SN models proposed by Woosley & Weaver (1995) and obtained fairly good agreement with a particular model. Furthermore, since the mass of the ISM eventually swept up by the preceding SN is theoretically derived, the observed abundance ratios with respect to H on the surfaces of stars would allow the mass of each element contained in the SN to be estimated (ST98; Tsujimoto & Shigeyama, 1998). A comparison of the derived mass of each element with that of theoretical SN models also leads to the conclusion that the apparent enhancement of $`\alpha `$-elements in the above two stars is a result of the deficiency of Fe in the ejecta of the preceding SN (see §2 for quantitative discussion). This indicates that there existed subluminous SNe with a small amount of 56Ni as a heat source, which decays to 56Fe through 56Co with a half-life of $``$77 d.

The least luminous SNe ever observed, SN 1997D (Turatto et al., 1998) and SN 1999br (Hamuy, 2002, and references therein), are in fact Fe-poor. From analyses of the light curve of SN 1997D (Turatto et al., 1998; Benetti et al., 2001), it has been found that this SN ejects only $`0.002M_{}`$ of 56Ni. The explosion energy $`E`$ is estimated to be $`E4\times 10^{50}`$ erg and the progenitor mass to be $`25M_{}`$. Hamuy (2002) reported that SN 1999br ejects a similar amount of 56Ni ($`0.0016M_{}`$) with an explosion energy of $`6\times 10^{50}`$ erg and that the ejected mass, which gives the lower bound of the progenitor mass, has been estimated to be 14 $`M_{}`$ from the light curve analyses. It is well known that there is a distinct class of SNe with similar progenitor mass that yield Fe comparable to SN 1987A or slightly higher ($`0.070.15M_{}`$, e.g., ST98). Such SNe will be referred to as normal SNe in the following. As mentioned in Turatto et al. (1998), the low Fe mass of SN 1997D may be due to fallback of material onto the collapsed remnant. In this case, Fe as well as other Fe-group elements such as Cr, Mn, Co and Ni produced in the deep interior of the core would fall back. In contrast, lighter elements synthesized in the outer layers of the star would be ejected and be excluded from any fallback. As a result, lighter elements are significantly enhanced with respect to Fe-group elements. The observed elemental abundance patterns of CS29498-043 and CS22949-037 exhibit characteristics that fit this explanation.

In the next section, the features of SNe that induced the formation of the stars CS29498-043 and CS22949-037 will be deduced from the observed elemental abundance pattern in light of the SN-induced star formation scenario proposed by ST98. Subsequently in §3, nucleosynthesis for Fe-group elements is discussed based on the implications of the \[Cr, Mn, Co/Fe\] ratios observed in these two stars, which is surprisingly at variance with the present theoretical expectations from spherically symmetric SN models.

## 2 Stars with relics of subluminous supernovae

The subluminous SNe discovered so far are less energetic varieties of SNe, although theoretically some energetic SNe with high masses such as 35 $`M_{}`$ could also be subluminous (Woosley & Weaver, 1995). Such low energy reduces the mass $`M_{\mathrm{SW}}`$ of ISM (mostly hydrogen) eventually swept up by a subluminous SN, obeying the relation $`M_{\mathrm{SW}}E^{0.97}`$ (ST98). Since the metallicity \[X/H\] of extremely metal-poor stars born from the ISM swept up by a SN remnant (SNR) is inversely proportional to $`M_{\mathrm{SW}}`$, stars that inherit the ejecta of a subluminous SN should have larger \[X/H\] as long as the other conditions are the same.

As discussed in ST98, the abundance ratios \[Mg/H\] in metal-poor stars can be used as an indicator of the progenitor mass of the SN inducing the formation of the stars, since theoretical models for massive stars predict that the mass of Mg inside a star increases with increasing stellar mass (e.g., Woosley & Weaver, 1995; Umeda & Nomoto, 2002). In addition, the amount of ejected lighter elements such as Mg can be predicted more precisely than that of heavier elements such as Fe because the ejection of the outer layer is less affected by theoretical uncertainty concerning the SN explosion mechanism (see, however, model Z20A of Woosley & Weaver, 1995, for example).

The relationships between the metallicity \[Mg/H\] of stars and the mass $`M_{\mathrm{ms}}`$ of the SN progenitor thus obtained are shown in Figure 1 for $`E=4\times 10^{50}`$ erg (solid line) and $`E=1\times 10^{51}`$ erg (dotted line), assuming the metallicity \[Mg/H\] of stars born from the shell formed by an SNR is approximated well by the average \[Mg/H\] inside the ISM swept up by the SNR. If CS29498-043 inherited the ejecta of a subluminous SN with explosion energy equal to that of SN 1997D, the SN progenitor mass becomes 25 $`M_{}`$, which intriguingly agrees with that inferred from the light curve analyses of SN 1997D. The SN progenitor mass for CS22949-037 estimated in the same manner is $`M_{\mathrm{ms}}20`$$`M_{}`$. The masses of ejected Fe are estimated to be 0.006 $`M_{}`$ (CS29498-043) and 0.003 $`M_{}`$ (CS22949-037) using the observed \[Fe/H\] ratios for the two stars and assuming $`E=4\times 10^{50}`$ erg. Both of the estimated Fe masses indicate subluminous SNe such as SNe 1997D and 1999br.

Figure 2 shows the enhancement factors of each element defined by the abundance ratios of the element with respect to Fe in the metal-poor stars discussed in this letter relative to those ratios in the ejecta of normal SNe with the corresponding progenitor masses. The abundances of normal SNe have been inferred from the elemental abundance pattern observed for metal-poor stars in the manner described in Tsujimoto & Shigeyama (1998). Thus derived abundance ratios correspond to the average values of metal-poor stars observed by McWilliam et al. (1995) at the metallicity of \[Mg/H\]$`2.3`$ and $`2.8`$, respectively (see Fig. 1). That is, an SN with $`E=1\times 10^{51}`$ erg and $`M_{\mathrm{ms}}=25`$ $`M_{}`$ will imprint a metallicity \[Mg/H\]$`=2.3`$ on the descendant stars, while an SN with the same energy but with $`M_{\mathrm{ms}}=20`$ $`M_{}`$ corresponds to \[Mg/H\]$`=2.8`$. As seen in Figure 2, the abundance of elements heavier than Si in these stars is at least a factor of 5 smaller than the average value of metal-poor stars, which is attributed to fallback onto the remnant in subluminous SNe. A moderate enhancement factor for Si in comparison with Fe-group elements might imply that fallback occurred inside the Si layer. The enhancement factors for Ca through Co have similar values in CS22949-037, while Cr through Fe are simlar in CS29498-043. A good coincidence in the enhancement factors is seen especially in CS22949-037, probably due to more accurate abundance determinations. This result suggests that there is little difference in the ratios of Ca, Cr, Mn, or Co to Fe between subluminous SNe and normal SNe with the same progenitor mass, with intriguing implications for nucleosynthesis, explosion, and mixing in SN models.

## 3 Implications for explosive nucleosynthesis and supernova explosion

Figure 3 shows the observed correlations between \[Cr/Fe\], \[Mn/Fe\], \[Co/Fe\] and \[Mg/H\] for metal-poor halo stars (McWilliam et al., 1995), CS29498-043, and CS22949-037. As shown in these figures, the observed \[Cr, Mn, Co/Fe\] trends are quite clear; the \[Cr, Mn/Fe\] ratios increase with increasing \[Mg/H\], whereas the \[Co/Fe\] ratio decreases. These clear correlations suggest that the \[Cr/Fe\], \[Mn/Fe\], \[Co/Fe\] ratios are closely correlated with the SN progenitor mass (Tsujimoto & Shigeyama, 1998; Nakamura et al., 1999), while the abundance ratios of $`\alpha `$-elements such as Mg, Si, and Ca to Fe do not have clear correlations. It is also clear from the same figure that the observed \[Cr/Fe\], \[Mn/Fe\], \[Co/Fe\] ratios of CS29498-043 and CS22949-037 do not follow the trend presented by the other stars. Five horizontal arrows in this figure indicate how these ratios would shift when these stars inherited the ejecta of an SN with $`E=1\times 10^{51}`$ erg instead of $`E=4\times 10^{50}`$ erg. The lengths of the horizontal shifts correspond to differences in the swept-up hydrogen mass between SNe with $`E=4\times 10^{50}`$ and $`1\times 10^{51}`$ erg. All the arrows pointing to the lines of the best fit for the observational data by $`\chi ^2`$ test, suggesting that the \[Mn/Fe\], \[Cr/Fe\], and \[Co/Fe\] ratios of these two stars do not deviate from those of normal SNe with the same progenitor mass, although for CS29498-043 an uncertainty remains due to the large error in the observed \[Mn/Fe\] ratio and a lack of Co data. As a consequence, there is little difference in the ratios of Cr, Mn, or Co to Fe between subluminous SNe and normal SNe given the same progenitor mass. This turns out to be incompatible with spherically symmetric SNe, as follows.

The elements Mn and Cr are the products of incomplete Si burning, whereas Co is produced by complete Si burning. Fe is produced in both burning processes. The region of complete Si burning is located at the bottom of the ejecta. If an SN explosion were spherically symmetric, different amounts of Fe could only be attributable to different combinations of the mass cut separating the ejecta from the remnant and the explosion energy. Since the mass cut determines the ejected fraction of the region in which complete Si burning takes place, smaller Co/Fe ratios would be expected from SNe with smaller amounts of Fe, corresponding to subluminous SNe. Larger Mn/Fe and Cr/Fe ratios would be expected from these SNe, as shown by Umeda & Nomoto (2002). Lower explosion energy would suppress complete Si burning and thus reduce both the amount of Fe and the Co/Fe ratio. Therefore, any combination of explosion energy and mass cut in spherically symmetric SNe is incompatible with the observed invariance of the abundance ratios of Fe-group elements to Fe. Since Rayleigh-Taylor instabilities take place after the nuclear reactions freeze out as was shown by a multi-dimensional calculation (Kifonidis et al., 2000), mixing does not alter the relative abundances of Fe-group elements determined in the way mentioned above. The same calculation has also shown that the mixing process is effective. Thus it is expected that fallback that may occur afterwards would hardly change these relative abundances, either. Mixing is unlikely to explain the observed invariance of the relative abundances among SNe with different explosion energies.

We will discuss a possibility that the observed ratios among Fe-group elements for the two stars reflect those in the ISM swept up by an SNR rather than those in individual SN ejecta. Indeed, the gas swept up by the SNR from which stars are formed contains not only heavy elements ejected from the SN itself but also those already present in the ISM. In general, for metal-poor regime such as \[Mg/H\] $`\stackrel{<}{}`$–2, the stellar abundance pattern is essentially determined by the SN yield due to a negligible contribution from the ISM (ST98; Tsujimoto, Shigeyama, & Yoshii, 1999). A lack of iron-group elements in the ejecta of a sub-luminous SN could enhance contribution from the ISM to the stellar abundance of iron-group elements. Since the \[Cr/Fe\], \[Mn/Fe\], \[Co/Fe\] ratios in the ISM should be equal to nearly zero, which corresponds to the plateau value seen in halo stars with –2$`<`$\[Fe/H\]$`<`$–1, the two stars discussed in this letter should show these abundance ratios nearly equal to the solar values if this is the case. However, the observed ratios shown in Figure 3 do not support this possibility. Therefore it is likely that these stars contain little heavy elements from the ISM. Furthermore, combinations of these ratios predicted from spherically symmetric SN models with those in the ISM can not make the observed ones.

Stellar roation might be a factor to produce a variety in energies of supernova explosions with the same progenitor mass. It is sometimes suggested that a rapid rotaion leads to a jet-like explosion (e.g., Woosley, 1993). In fact, it has been a preferred explanation for matter mixing suggested by the emission line profiles observed for SN 1987A (Nagataki, 2000) and SN 1998bw (Maeda et al., 2002), both of which are core-collapse SNe. Furthermore, the optical spectropolarimetry of several core-collapse SNe suggests highly asymmetric explosion mechanisms such as the formation of bipolar jets (Wang et al., 2001). Detailed numerical studies on nucleosynthesis and ejection mechanism in jet-like explosions might be required to explain the observed trend of Fe-group elements in halo field stars.

An alternative solution may be in a hot bubble formed by neutrino emission. If neutrinos from the nascent neutron star deposit substantial energy to make a hot bubble behind the blast wave, the hot bubble may carry exclusively the product of complete Si burning to the outer layer by buoyancy. A numerical simulation by Wilson & Mayle (1988) has shown that this mechanism leads to a less energetic SN. This may suppress fallback of Co and Fe in subluminous SNe, though other numerical studies have shown that a hot bubble does not result in a SN explosion (e.g., Bruenn, Mezzacappa, & Dineva, 1995). The current understanding of the involved physics is, however, still too premature to extract a definite conclusion on this issue.

## 4 CONCLUSIONS

The recently discovered new class of carbon-rich stars including CS29498-043 and CS22949-037, which exhibit a significant enhancement of $`\alpha `$-elements with respect to Fe-group elements, has been shown in this study to inherit the elemental abundance pattern of a precursor SNe that ejected a small but adequate amount of Fe with low explosion energy. These SNe are found to resemble the two recently discovered subluminous SNe 1997D and 1999br in terms of the ejected Fe mass, explosion energy, and progenitor mass. Subluminous supernovae may therefore be studied more easily by looking at the surface of nearby metal-poor stars than by searching for similar occurrences elsewhere.

The present study revealed that there is little difference in the \[Cr/Fe\], \[Mn/Fe\], \[Co/Fe\] ratios between subluminous SNe and normal SNe given the same progenitor mass, which has intriguing implications in nucleosynthesis and explosion in SN models. The observed facts are apparently at odds with the prediction by spherically symmetric SN explosion models.

We are grateful to the anonymous referee for making useful comments. This work was supported in part by a Grant-in-Aid for Scientific Research (11640229, 12640242) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan.


format: text

text:
# Stellar Rotation from GAIA Spectra

## 1. Introduction

At present there are approximately 20.000 stars with measured $`v_{rot}\mathrm{sin}i`$ (Glebocki, Gnacinski and Stawikowski 2000). Sorting them by stellar type and rotational velocity (Munari 2002, Soderblom 2001) shows that early type stars (O, B, A, early F) have high rotational velocities (50-400 km/s) and late type stars (late F, G, K, M) are slow rotators (in majority with $`v_{rot}\mathrm{sin}i<`$50 km/s), which is mainly attributed to the fact that late type stars have convective envelopes, while early type stars do not. In view of measuring stellar rotational velocities with GAIA, it is necessary to estimate the accuracy of measured $`v_{rot}\mathrm{sin}i`$, which helps us in answering the question, whether slow rotators can be discerned from non-rotators with GAIA.

## 2. Scientific questions related to stellar rotation

Data obtained by GAIA with millions and millions of stellar spectra and measured rotational velocities will certainly crucially contribute to the current knowledge on stellar rotation and its influence on various aspects of stellar physics. To shortly name a few (Maeder and Meynet 2000, Soderblom 2001, Patten and Simon 1996):

* The effect of stellar rotation on stellar structure and evolution: rotation of a star influences its structure, luminosity, position on the HR diagram, life time etc. Rotation can via rotationally induced mixing lead to He and N enrichment, chemically peculiar stars, it can cause turbulence, influence stellar winds and spots.
* The impact of differential rotation on stellar structure, surface phenomenology, models for generating magnetic fields etc.
* Stellar rotation as the indicator of age (for low mass ZAMS stars) through the magnetic breaking and the indicator of photospheric activity.
* In binary systems: stellar rotation is important in understanding the number of questions, such as the question of synchronization between rotational and orbital periods and the effectivness of tidal energy dissipation.
* Rotational velocity as an indicator of angular momentum of solar type stars and how it is affected by the presence of massive distant planets.
* Stellar rotation as an indicator of age, total mass and binding energy in open clusters and as a clue to test various theories of fragmentation of the parent cloud.
* The orientation of rotational axis - are they really randomly oriented and therefore not connected with the rotation of the Galaxy, which is the question that only wide GAIA statistics can answer.

## 3. Simulating GAIA accuracy of $`v_{rot}\mathrm{sin}i`$

The important issue of GAIA contribution to stellar rotation physics is the accuracy of obtained $`v_{rot}\mathrm{sin}i`$. Here are presented results of simulating the determination of rotational velocities on synthetic stellar spectra. We used Kurucz synthetic stellar spectra database in the GAIA wavelength region.

To start the simulation we choose the spectrum with known T, Z, log g, vradand original $`v_{rot}\mathrm{sin}i`$, artificially add Poison distributed noise and afterwards fit it with noise-free spectra with various $`v_{rot}\mathrm{sin}i`$. As the best fitting spectrum we take the spectrum with minimum $`\chi ^2`$ and repeat the test N times. The accuracy is determined as the standard deviation error between the original and recovered $`v_{rot}\mathrm{sin}i`$:

$$\sigma ^2=\frac{1}{N}\underset{i=1}{\overset{N}{}}\left[v_i(orig)v_i(rec)\right]^2.$$

We used synthetic spectra of four star types: late G type giant (T=4750, Z=-0.5, log g=1.0), K dwarf (T=4750, Z=-0.5, log g=4.5), G type giant (T=5500, Z=-0.5, log g=2.0) and early F type star (T=7250, Z=0.0, log g=4.5), all with $`v_{rot}\mathrm{sin}i`$=10 km/s at different resolutions: R=5000, 8615, 17230 and 20000. Noise added corresponds to signal to noise ratio in the range of S/N=10-300 and the number of trials is N=5000.

As the fitting spectra we first use the spectra with stellar parameters (T, Z, log g, vrad) exactly the same as in the original spectrum, so that they differ only in $`v_{rot}\mathrm{sin}i`$. The obtained accuracy in such simulations is shown on Figure 1: as expected the accuracy is much better for high resolution and improves with better signal to noise. It should be stressed though, that this is the most ideal case, since other stellar parameters are precisely known.

Due to possible uncertainties in these, the accuracy of rotational velocity will generally become worse. To estimate the influence of those uncertainties on the accuracy of obtained rotational velocity, we performed the same simulations by fitting the original spectrum with spectral templates that had one of the parameters (T, Z, log g or vrad) offset from its true value.

The least crucial factor turns out to be the uncertainty in log g. The dotted curves in Figure 2 show the accuracy of rotational velocity if fitted by spectra having log g offset for 0.5 from its original value.

More crucial is an accurate temperature. The error of 250 K in some cases keeps the rotational velocity error at 5-10 km/s, not improving with S/N, while the error of 125 K (Figure 3) is small enough that at least at high resolutions the accuracy is better than 1 km/s at high S/N.

One of crucial factors is an accurate metallicity. The error of 0.5 in metallicity leads to inaccurate rotational velocities (error $`>`$10 km/s), the smaller error of 0.25 is about a factor of 2 better and 0.1 is quite good for high resolutions (Figure 4). Note that the accuracy of the rotational velocity does not seem to suffer so much if the metallicity is underestimated than if it is overestimated by the same amount.

Crucial is also the accuracy of the obtained radial velocity. As shown in Figure 5, the error in radial velocity of 10 km/s leads to an error in rotational velocity of more than 10 km/s. If the radial velocity error is 5 km/s the error of rotation velocity is about 2 times smaller.

## 4. Conclusion

Table 1 gives the summary of stellar rotation accuracy obtained by simulations based on synthetic stellar spectra for different resolution and signal to noise ratio in the continuum. Columns 3-7 report the rotational velocity errors if the template used to recover the spectrum is ideal or mismatched for a given amount in temperature, metallicity, surface gravity or radial velocity.

The results we obtained can be used only as a rough indicator of GAIA capabilities. We used synthetic stellar spectra without allowing for modeling uncertainties, chemical composition peculiarities etc. The results of simulations show that effects of combined errors in two of the template parameters (e.g. T and Z, or Z and log g) on the accuracy of rotational velocity are not always easily predictable. One should also be aware that the crowding of stellar spectra in GAIA focal plane may smear out some line profile details, which may be crucial in determining the rotational broadening of spectral lines. Nevertheless, we believe that measuring even small rotational velocities with GAIA is not out of reach, provided that the resolution is better than 10000, or preferably better than 15000.

## References

Glebocki R., Gnacinski P., Stawikowski A. 2000, AcA, 50, 509G

Maeder, A., Meynet, G. 2000, ARA&A38, 143

Munari, U. 2002, 2nd GAIA RVS Workshop, Asiago

Patten, B. M., Simon, T. 1996, ApJS106, 489

Soderblom, D. R. 2001 in Encyclopedia of Astronomy and Astrophysics 3137, ed. A. Baker, D. Emerson, N. Hankins & J. Matthias, IOP Publishing and Nature Publishing Group.


text:
# Entanglement versus relaxation and decoherence in a quantum algorithm for quantum chaos

## Abstract

We study analytically and numerically the behavior of the concurrence (a measure of the entanglement of formation) of a pair of qubits in a quantum computer operating an efficient algorithm for quantum chaos. Our results show that in an ideal algorithm the entanglement decays exponentially with the diffusive relaxation rate induced by classical chaos. This decay reaches a residual level which drops exponentially with the number of qubits $`n_q`$. Decoherence destroys the residual entanglement with a rate exponential in $`n_q`$.

Enormous interest into quantum information and computation Nielsen and Chuang (2000) has generated serious efforts in characterizing and understanding quantum entanglement which is considered as the ultimate origin of quantum power (see a recent review Bruß (2002)). A quantitative measure of the entanglement of formation, namely the concurrence $`C`$, was introduced and shown to be able to characterize an arbitrary state of two qubits Hill and Wootters (1997); Wootters (1998). Being closely related to the von Neumann entropy $`S`$ of the reduced density matrix $`\rho `$ of two qubits111Given a density matrix $`\rho `$ for a pair of two-level systems, the concurrence $`C`$ is defined as $`min\{0,\lambda _1\lambda _2\lambda _3\lambda _4\}`$, where the $`\lambda _i`$ are the eigenvalues, in decreasing value order, of the hermitian matrix $`\sqrt{\sqrt{\rho }\stackrel{~}{\rho }\sqrt{\rho }}`$, with $`\stackrel{~}{\rho }=(\sigma _y\sigma _y)\rho ^{}(\sigma _y\sigma _y)`$. The von Neumann entropy $`S`$ of $`\rho `$ monotonically increases from 0 to 1 when $`C`$ goes from 0 to 1 Hill and Wootters (1997); Wootters (1998)., this quantity was recently found to have interesting applications to quantum phase transitions in interacting spin systems Osterloh et al. (2002). In parallel, the properties of entanglement were investigated in a quantum model of coupled tops where it was shown that there exists a typical value of entanglement which is determined by the chaotic behavior of the dynamics of the model Miller and Sarkar (1999) and that the growth rate of entanglement of initially decoupled tops is increased by the underlying classical chaos Bandyopadhyay and Lakshminarayan (2002). On the same line, it was recently shown that contrary to intuition even a heat bath may create entanglement between two qubits Braun (2002).

All these studies Osterloh et al. (2002); Miller and Sarkar (1999); Bandyopadhyay and Lakshminarayan (2002); Braun (2002) clearly demonstrated how rich entanglement properties can be in interacting quantum systems. However, in the context of quantum computation it is much more crucial to analyze the evolution of entanglement in a specific algorithm performing an operational task. Indeed, it is expected that the entanglement is very sensitive to noise and decoherence Paz and Zurek (2001); Schack and Caves (1993, 1996) and the understanding of its behavior in an operating algorithm can lead to better strategies in the control of decoherence and imperfection effects. As for our knowledge such direct investigations have not been performed until now. Therefore, in this paper we study the behavior of the concurrence in an efficient algorithm for the quantum sawtooth map which has been proposed recently in Benenti et al. (2001). The algorithm for this model has a number of important advantages: all $`n_q`$ qubits are used in an optimal way and no ancillae are required, one map iteration in the Hilbert space of size $`N=2^{n_q}`$ is performed in $`O(n_q^2)`$ quantum gates, and the algorithm is based on the quantum Fourier transform (QFT) which is one of the main elements of various quantum algorithms Nielsen and Chuang (2000). This allows to simulate a complex dynamics in the regime of quantum chaos with a small number of qubits. Since the entanglement can be efficiently measured experimentally (see e.g. Horodecki and Ekert (2002)) the experimental observation of the concurrence behavior discussed here can be realized on NMR Weinstein et al. (2001); Vandersypen et al. (2001) or ion-trap Gulde et al. (2003) based quantum computers with about 6 - 10 qubits and a few hundreds of gates.

Contrary to the situation discussed in Miller and Sarkar (1999); Bandyopadhyay and Lakshminarayan (2002), our results show that in the exact quantum algorithm the underlying classical chaos leads to an exponential decrease of the concurrence $`C`$ down to some residual level $`\overline{C}`$ which characterizes the global system coherence. On the other hand, the presence of noise in the quantum gates leads to a destruction of this coherence with a rate $`\mathrm{\Gamma }`$ growing exponentially with the number of qubits $`n_q`$. This shows that entanglement can be very sensitive to decoherence.

The dynamics of the classical sawtooth map Dana et al. (1989); Chen et al. (1990) is given by:

$$\overline{n}=nk\frac{dV(\theta )}{d\theta },\overline{\theta }=\theta +T\overline{n}mod2\pi $$
(1)
where $`V(\theta )=\theta ^2/2`$, $`pi\theta <\pi `$ and the bars denote the variables after one iteration. After rescaling $`y=Tn`$ and $`x=\theta `$, it is clear that the dynamics depends only on the parameter $`K=kT`$. Due to the discontinuity in the derivative of $`V(\theta )`$, the Kolmogorov-Arnold-Moser (KAM) theorem cannot be applied to the map (1) and its dynamics becomes chaotic and diffusive for arbitrarily small values of the chaos parameter $`K>0`$ Dana et al. (1989). For $`K1`$ the diffusion is governed by a nontrivial cantori regime which has been worked out in Dana et al. (1989). In this case the rescaled diffusion rate $`D_0(K)=(\mathrm{\Delta }y)^2/t1.2\pi ^2K^{2.5}/3`$ is much smaller than the quasi-linear diffusion rate corresponding to the random phase approximation $`D_{ql}=\pi ^2K^2/3`$ (the latter becomes valid only at $`K1`$). The diffusion rate in $`n`$ is $`D=D_0(K)/T^2`$.

The quantum sawtooth map Benenti et al. (2001); .Borgonovi et al. (1996); Borgonovi (1998); Casati and Prosen (1999); Prange et al. (1999) is the quantized version of the classical map, to which it corresponds in the limit $`k\mathrm{}`$, $`T0`$, and $`K=kT=\text{const}`$. One step in the quantum map is given by the unitary operator $`\widehat{U}`$ acting on the wave function $`\psi (\theta )`$:

$$\overline{\psi (\theta )}=\widehat{U}\psi (\theta )=e^{ik\theta ^2/2}e^{iT\widehat{n}^2/2}\psi (\theta ),$$
(2)
where $`\widehat{n}=i/\theta `$ (we set $`\mathrm{}=1`$); when $`n_q`$ qubits are used for the discretization, so that the total number of levels is $`N=2^{n_q}`$, we assume periodic boundary conditions on a torus in the phase representation ($`\psi (\theta +2\pi )=\psi (\theta )`$) and in the momentum representation ($`\psi (n+N)=\psi (n)`$), as discussed in Benenti et al. (2001); Izrailev (1990). As a result, the phase $`\theta `$ takes only $`N=2^{n_q}`$ discrete equidistant values in the interval $`\pi \theta <\pi `$, and so does the momentum ($`n=0,1,\mathrm{},N1`$, i.e. $`n+N`$ is identified with $`n`$). The sawtooth map approximately describes the dynamics in the stadium billiard .Borgonovi et al. (1996), the phenomenon of dynamical localization which is similar to the Anderson localization in disordered potentials and cantori induced localization Borgonovi (1998); Casati and Prosen (1999); Prange et al. (1999). Thus this simple map describes a rich and complex dynamics and represents an interesting testing ground for efficient quantum computation. It is especially important to understand how this complex quantum dynamics will be affected by imperfections in realistic quantum computers.

The numerical simulation of the map (2) is based on the quantum algorithm described in Benenti et al. (2001) and is implemented on the basis of a quantum computer language developed in Bettelli et al. (2003) which is well adapted for the experimental operation of quantum gates via classical computer software. In this way the dynamics of up to 20 qubits can be easily simulated on a laptop.

To investigate the behavior of the concurrence in the quantum map (2) we compute $`C`$ for the two most significant qubits which determine the first two binary digits $`a_{1,2}`$ in the expansion of momentum $`n`$: the reduced density matrix $`\rho `$ for this qubit pair is obtained by tracing out all other $`n_q2`$ less significant qubits (the digits $`a_i`$ with $`3in_q`$ in the expansion of $`n=(a_1a_2a_3\mathrm{}a_i\mathrm{}a_{n_q})`$). After that $`C`$ is computed from $`\rho `$ as described in Hill and Wootters (1997); Wootters (1998). In this way we obtain the concurrence value $`C`$ on a global scale of the whole system which is decomposed into 4 equal parts with $`N/4`$ quantum states in each of them. In addition we fix $`T=2\pi L/N`$ in the regime of quantum resonance so that $`L`$ gives the integer number of classical phase space cells embedded in the quantum torus of size $`N`$ Benenti et al. (2001); Izrailev (1990) (the classical dynamics is periodic in $`n`$ with period $`2\pi /T`$). In the following we also take $`L`$ to be a multiple of 4 to have an integer number of classical cells in the 4 parts of the partition in the momentum $`n`$.

Typical examples of the dependence of $`C`$ on the number of map iterations $`t`$ are shown in Fig.1. According to these data $`C(t)`$ decays exponentially down to a residual value $`\overline{C}`$ and, in the limit of large $`N`$, the decay rate $`\gamma `$ becomes independent of $`N`$. It is natural to compare this rate with the rate of classical relaxation. Indeed, due to underlying classical chaos, the probability distribution $`f_n=|\psi _n|^2`$ over $`n`$ is described by the Fokker-Planck equation $`f_n/t=D^2f_n/^2n/2`$ which gives the relaxation to equipartition with the rate

$$\gamma _c=2\pi ^2D/N^2=D_0(K)/2L^2.$$
(3)
The comparison between this classical value $`\gamma _c`$ and the rate of concurrence decay $`\gamma `$ is given in Fig.2. It clearly shows that the decay rate of $`C(t)`$ is given by the classical rate: $`\gamma =\gamma _c`$. It is important to stress that this relation remains valid also in the nontrivial cantori regime ($`K1`$) and that the quantum decay reproduces all oscillations of the classical diffusion (see inset in Fig.2).

The properties of the residual value of the concurrence $`\overline{C}`$ are analyzed in Fig.3. We will argue in the following that they can be understood in terms of the system conductance $`g`$ (Edwards and Thouless (1972); Thouless (1974)); in view of this we express $`\overline{C}`$ vs. $`g=2\gamma _c/\mathrm{\Delta }=ND_0(K)/L^2`$ where, up to a constant factor, $`\mathrm{\Delta }=1/N`$ is the level spacing and $`2\gamma _c=D_0(K)/L^2`$ is the Thouless energy (see Mirlin (2000) for a recent review). In spite of strong fluctuations the data presented in Fig.3 can be described by the global average dependence $`\overline{C}1/\sqrt{g}1/\sqrt{N}`$. Indeed, for $`K=0.5,L=4`$ the system size varies by 3 orders of magnitude and the fit gives an algebraic decay with power $`\alpha =0.56\pm 0.02`$ being close to $`1/2`$. We attribute the presence of strong fluctuations to the fact that the value $`\overline{C}`$ is averaged only over time but there is no averaging over parameters. Thus, from the point of view of disordered systems Mirlin (2000) $`\overline{C}`$ represents only one value for one realization of disorder.

We propose the following explanation of the results presented in Figs.1-3. For a state $`|\psi `$ like in Fig.1, we can write $`|\psi =_{a_1a_2}|a_1a_2|\varphi _{a_1a_2}`$ where $`a_{1,2}=0`$ or $`1`$. Then the value of the concurrence $`C`$ is proportional to the difference of two scalar products, $`C|Q_{14}Q_{23}|`$, where $`Q_{14}=2\sqrt{|\varphi _{00}|\varphi _{11}|^2}`$ and $`Q_{23}=2\sqrt{|\varphi _{01}|\varphi _{10}|^2}`$. From this relation and the fact that the initial state is symmetrically distributed with respect to the transformation $`nNn`$, it follows that $`C`$ is proportional to the difference $`|W_{11}+W_{00}W_{01}W_{10}|`$ where $`W_{a_1a_2}`$ is the total probability inside the part $`(a_1a_2)`$. In the classical limit this probability difference relaxes to zero with the classical relaxation rate $`\gamma _c`$, and that’s why $`\gamma =\gamma _c`$ in agreement with the data of Fig.2.

The residual value $`\overline{C}`$ is determined by the quantum fluctuations of the previous difference of scalar products. In fact, due to the discretization of map (2), the symmetry $`nNn`$ is broken and $`|\varphi _{00}`$ becomes different from $`|\varphi _{11}`$. Therefore in the scalar product $`Q_{14}`$ (and $`Q_{23}`$) the $`N/4`$ terms have random signs and thus $`Q_{14}1/\sqrt{N}`$ (each term is of the order of $`1/N`$). In this estimate we assumed a summation over all $`N`$ wave function components. However for finite values of the conductance $`g`$ only the states inside the Thouless energy interval $`2\gamma _c`$ have a significant scalar overlap Mirlin (2000); Blanter (1996), and thus we can make a conjecture that $`N`$ should be replaced by the effective number of components which is of the order of $`N_{\text{eff}}\gamma _c/\mathrm{\Delta }g`$. According to this, $`\overline{C}1/\sqrt{g}`$ in agreement with the data of Fig.3.

The existence of a residual level of concurrence for an ideal quantum algorithm reflects the fact that the global behavior of the whole system remains coherent. In fact, the Poincaré theorem guarantees that for very large times the concurrence will have a revival close to the initial value (however, this will happen on an exponentially large time scale). The situation becomes qualitatively different in presence of external decoherence represented by noisy gates. In our numerical simulations, noisy gates are modeled by unitary rotations by an angle randomly fluctuating in the interval $`(ϵ/2,ϵ/2)`$ around the perfect rotation angle. The presence of this external decoherence leads to a decrease of the residual value of $`C`$ as illustrated in Fig.4: the constant level is replaced by an exponential decay which gives $`\overline{C}\mathrm{exp}(\mathrm{\Gamma }t)`$.

In order to obtain the dependence of $`\mathrm{\Gamma }`$ on the parameters we extracted it from the fit of the averaged ratio of $`C`$ under a noisy evolution to its value in the ideal one. An example of such ratio and the corresponding fit is shown in Fig.4. To suppress fluctuations we averaged over $`20`$ realizations of the noisy evolution. Moreover, the fit was restricted to the plateau regime, where the exact concurrence is fluctuating around its residual value (the initial diffusive relaxation was excluded from the fit).

The results for $`\mathrm{\Gamma }`$ obtained in this way are presented in Fig.5. Quite naturally we find that $`\mathrm{\Gamma }ϵ^2`$, as it was also seen in other simulations of quantum algorithms with noisy gates (e.g. Chepelianskii and Shepelyansky (2002)).This scaling becomes better and better for large $`ϵ`$ values where $`\mathrm{\Gamma }`$ is larger. However, more surprisingly there is an exponential growth of $`\mathrm{\Gamma }`$ with the number of qubits $`n_q`$ ($`\mathrm{\Gamma }\sqrt{N}`$). This result is very different from those obtained in Benenti et al. (2001); Chepelianskii and Shepelyansky (2002), where the time scale for fidelity and the decoherence rate for tunneling oscillations varied polynomially with $`n_q`$. We see two possible reasons for the exponential sensitivity of the residual concurrence to decoherence. At first, in our case $`\mathrm{\Gamma }`$ is computed over a very large time interval, for which the quantum dynamics already reached its asymptotic behavior (plateau for the residual concurrence); it is known that on very large times the eigenstates are exponentially sensitive to imperfections due to the chaotic structure of the wave functions (see results and discussions in Benenti et al. (2002)). Another possible reason can be related to the fact that the residual value of the concurrence on the plateau is on its own exponentially small and maybe this is the reason why it becomes so sensitive to decoherence. Further investigation on the decoherence effects for the concurrence are required to understand in a better way this exponential sensitivity of $`C`$.

In our studies we restricted ourselves to the investigation of entanglement only between two qubits. The problem of the characterizing the entanglement of a larger number of qubits represents an interesting challenge but at the same time it is much more complicated Bruß (2002). However even the relatively simple case of two qubits shows nontrivial links between concurrence and such interesting physical phenomena as quantum phase transitions Osterloh et al. (2002) and statistical relaxation.

In summary, our studies show that the decay of the concurrence in an operating quantum computer is determined by the underlying relaxation rate of the classical dynamics. We show that the residual level of entanglement in an ideal algorithm scales as the inverse square root of the conductance of the system. This residual entanglement is destroyed by decoherence, whose effective rate grows exponentially with the number of qubits.

This work was supported in part by the EC contracts RTN QTRANS and IST-FET EDIQIP and the NSA and ARDA under ARO contract No. DAAD19-01-1-0553.


text:
# Large-Scale Simulations of Diffusion-Limited 𝑛-Species Annihilation

## I introduction

Diffusion-limited reactions have attracted much interest in recent years reactions ; kang85 . The kinetics of such systems is dominated by local fluctuations in the concentration of the reactants, thus posing a formidable problem which has not yet been solved: there exists no comprehensive theoretical approach for the analysis of diffusion-limited processes.

Few select models are amenable to exact analysis. These include one-species annihilation, $`A+A0`$, and two-species annihilation $`A+B0`$ (see dbook  and references therein). In one dimension, the particle density for one-species annihilation decays as $`c(t)t^{1/2}`$, while for two-species annihilation (with equal initial concentrations and same diffusion constants for the two species) $`c(t)t^{1/4}`$. In either case the result is strikingly different from the mean-field kinetics of the corresponding reaction-limited process, $`c(t)1/t`$. To bridge the gap between these disparate behaviors, ben-Avraham and Redner dba  proposed the n-species annihilation model, where particles belonging to the $`n`$ species $`A_1,A_2,\mathrm{},A_n`$ diffuse on the line and react immediately uppon encounter, according to the scheme:

$$A_i+A_j0,i,j=1,2,\mathrm{},n,ij.$$
(1)
For $`n=2`$ we recover two-species annihilation, while in the limit $`n\mathrm{}`$ encounters between like-particles are improbable and the model is equivalent to one-species annihilation. For intermediate values of $`n`$, one expects $`c(t)t^{\alpha (n)}`$.

In dba  it was proposed, following a heuristic scaling argument and treating fluctuations via the Van Kampen $`\mathrm{\Omega }`$-expansion vankampen , that $`\alpha (n)=(2n3)/(4n4)`$. This was supported by numerical simulations of lattices of typically $`10^6`$ sites, and up to $`10^6`$ time steps. (In one time step, all of the particles in the system move one lattice spacing each, on average.) Recently, we have conducted extensive numerical simulations dexin , following the method of Renormalized Reaction-Cells (RRC) benavraham87 ; dawkins00 ; shafrir . The systems involved are up to $`2^{28}2.7\times 10^8`$ sites long, and the processes were simulated to completion (until no further reactions are possible), for up to $`10^{15}`$ time steps. The new data leads us to the conjecture that $`\alpha (n)=(n1)/2n`$ dexin . We also find a correction to the main decay mode, of the form $`c(t)t^{\alpha (n)}(A+Bt^{\mathrm{\Delta }(n)})`$, $`\mathrm{\Delta }(n)=1/2n`$. The same results were found, independently (and unbeknownst to us), by Deloubrière et al., tauber . In their theoretical derivation, they consider a simplified version of $`n`$-species annihilation, where domains of alternating species loose particles to reactions at one and the same rate, in a synchronous fashion. The approximation is more than reasonable, yet it does not rigorously apply to the original model, and analysis of corrections is certainly beyond its scope. Moreover, the simulations in tauber  are comparable in size to those in dba . In what follows, we report the results of our large-scale simulations, which strongly support the conclusions of tauber . We also propose a scaling relation for the correction exponent $`\mathrm{\Delta }`$ for $`n`$-species annihilation, and possibly for other reaction models where particles segregate into distinct domains.

## II scaling

As is well known, local fluctuations in the concentrations of the different species drive the kinetics of $`n`$-species annihilation kang85 ; dbook ; dba . An initially random homogeneous distribution of the particles evolves into a continuously growing mosaic of domains of alternating surviving species. Two lengthscales characterize the emerging distribution and dominate the system: the inter-domain distance — the distance between the last particle in a domain and the first particle in the domain next to it — $`\mathrm{}_{\mathrm{AB}}(t)`$, and the domain length, $`\mathrm{}(t)`$ redner . These quantities grow with time as

$$\mathrm{}(t)t^\beta ,\mathrm{}_{\mathrm{AB}}(t)t^\gamma .$$
(2)
Once domains form, reactions might take place only at the domain boundaries, and particles have to diffuse across the domain gap $`\mathrm{}_{\mathrm{AB}}`$ to react with other species. This takes a typical time of $`\mathrm{\Delta }t\mathrm{}_{\mathrm{AB}}^2/D`$, where $`D`$ is the diffusion constant. The change in particle concentration during time $`\mathrm{\Delta }t`$ equals the total number of domain boundaries divided by the lattice size $`L`$; $`\mathrm{\Delta }c(L/\mathrm{})/L=1/\mathrm{}`$. Thus,

$$\frac{\mathrm{\Delta }c}{\mathrm{\Delta }t}\frac{D}{\mathrm{}\mathrm{}_{\mathrm{AB}}^2}.$$
(3)
On substituting the relations (2) and $`c(t)t^\alpha `$, we derive the scaling rule

$$2\gamma +\beta \alpha =1.$$
(4)
Due to the underlying transport mechanism, we expect that domains grow diffusively, as $`\mathrm{}t^{1/2}`$, so $`\beta =1/2`$, and in effect there is only one independent exponent: $`2\alpha \gamma =1/2`$. The general scaling form holds also for two-species annihilation in the presence of drift (and with ard-core repulsion between like species), where $`\alpha =1/3`$, $`\beta =7/12`$, and $`\gamma =3/8`$ shafrir .

## III simulation results

The $`n`$-species annihilation process is simulated as follows. The sites of a one-dimensional lattice are either empty or occupied by a particle (of one of the $`n`$ species). Periodic boundary conditions are imposed, so the lattice is effectively a ring. At each Monte Carlo step a particle is chosen randomly and is moved to the nearest site to its right or left, with equal probabilities. If the target site is occupied by a particle of a different species, then both particles are removed from the system, mimicking the reaction (1). If the target site is occupied by a particle of the same species, then the move is disallowed and it does not take place. Regardless of the outcome, time is incremented by $`1/N(t)`$, where $`N(t)`$ is the total number of extant particles.

As the simulation proceeds, the particle concentration declines and the typical distance between particles increases. The time spent on simulating the diffusive motion of the particles until they interact grows even faster, as the square of the distance between them. Because of that, computer simulations are limited to relatively short times. This problem is overcome by the RRC method benavraham87 ; dawkins00 ; shafrir .

In the RRC method the particles occupy cells, rather than sites. Each time the concentration halves, the cells are renormalized: every two cells are merged into one, and time is renormalized accordingly. The typical time required to diffuse out of a renormalized cell twice as large as that of the previous generation is four times longer. Thus, physical time is simulated faster with each renormalization step and the process can be simulated to completion. Other details for the implementation of the RRC method are discussed in shafrir .

Simulations were performed on DEC Alpha processors running Linux. Since each lattice site requires 6 bytes (for species, number of particles, and a pointer to a list of populated sites that is used for fast selection at each Monte Carlo step), with 2 Gigabytes ($`2^{31}`$ bytes) memory we were able to simulate lattices of up to $`2^{28}`$ sites. The compiler was given special $`\mathrm{\#}`$pragma pack(1) instructions to circumvent word alignment (which would allocate $`32`$ bytes for our $`6`$-byte site).

To test the technique, we have simulated the cases of $`n=2`$ and $`n=3`$ on lattices of $`2^{16}=65,536`$ sites, in both the RRC and the traditional simulation method. These lattices are small enough to enable the simulation of the process by the traditional method to completion. On the other hand, the system is large enough to let us examine the effect of the renormalizations: with $`2^{16}`$ sites and $`c(0)=1/16`$ the RRC method requires $`12`$ renormalizations. In Fig. 1a we compare the particle concentration as obtained by the two methods. In Fig. 1b we plot the local slope of the curves, the exponent $`\alpha (t)`$. The renormalizations are discernable only in this second, more stringent test, but the overall agreement is excellent. Similar results were obtained for the domain size and the distance between domains.

Having gained some confidence in the RRC method, we proceed to larger simulations. In Fig. 2 we show the surviving number of particles, $`N(t)`$, at time $`t`$, for $`n=3`$, $`4`$, and $`5`$, and several lattice sizes. In Fig. 3, we plot the local decay exponent $`\alpha (t)`$ for our largest simulations of $`n=3`$. The maximum of the curve at $`t10^4`$ agrees with the earlier prediction that $`\alpha =3/8`$ dba . (Indeed, simulations in dba  yielded a somewhat smaller value than the theoretical $`3/8`$, in perfect agreement with current results.) However, $`\alpha (t)`$ is seen to diminish with time, suggesting a long-time asymptotic limit of $`\alpha 1/3`$. This limiting value is confirmed in the data collapse (especially at long times) of Fig. 4, where we plot $`t^\alpha c(t)`$ vs. $`t^\beta /L`$ for various system sizes, and $`\alpha =1/3`$, $`\beta =1/2`$. Independent measurements show that $`\beta =1/2`$, as assumed, to within $`2\%`$, and the data collapse of Fig. 4 deteriorates with other choices for the values of $`\alpha `$ and $`\beta `$.

We have analyzed in this fashion $`n=3`$, $`4`$, and 5-species annihilation, and measured the exponents $`\alpha `$, $`\beta `$, and $`\gamma `$. Our results are summarized in Table 1. In all cases, the scaling relation (4) seems to hold, and $`\beta =1/2`$ to within numerical errors. Looking for a simple expression of these results, that would have the appropriate limits for the known cases of $`n=2`$ (two-species annihilation) and $`n\mathrm{}`$ (one-species annihilation), we were led to the conjecture dexin

$$\alpha =\frac{n1}{2n},\beta =\frac{1}{2},\gamma =\frac{2n1}{4n},$$
(5)
a result derived independently by Deloubrière et al. tauber

Finally, let us address the issue of corrections to scaling of the concentration decay. We look for corrections of the form

$$c(t)t^a(A+Bt^\mathrm{\Delta }),$$
(6)
where $`A`$ and $`B`$ are constants. Our strategy consists of performing a least squares linear fit of $`A+Bt^\mathrm{\Delta }`$ to $`t^\alpha c(t)`$, for different powers $`\mathrm{\Delta }`$, and searching for the value of $`\mathrm{\Delta }`$ which minimizes the error. The scaling form (6) is expected to work only after the asymptotic regime sets in, and before finite-size effects begin, and the sticky part of our procedure is deciding which times demarcate this region. Experimenting with different choices gives us a feel for the errors involved. In Fig. 5 we show best fits for the region $`t=10^6`$$`10^{12}`$, for $`n=3`$ on a $`L=2^{28}`$ lattice, where our data is most reliable. The results are most compatible with $`\mathrm{\Delta }=1/6`$ (for $`n=3`$). Similar tests for other values of $`n`$ lead us to the conjecture that $`\mathrm{\Delta }(n)=1/2n`$.

The correction exponent can be understood by a simple-minded argument. In deriving (3) we have assumed that the typical distance between reacting particles, at the edges of adjacent domains, is $`\mathrm{}_{\mathrm{AB}}`$. While this is correct, we note that, had the distribution of particles been homogeneous, the distance between reacting pairs would be typically $`\mathrm{}_{AA}L/ct^\alpha `$, quite different from the assumed $`\mathrm{}_{\mathrm{AB}}t^\gamma `$. Using $`\mathrm{\Delta }t\mathrm{}_{AA}^2/D`$ in (3), instead of $`\mathrm{}_{\mathrm{AB}}^2/D`$, yields a faster decay; $`ct^{(1\beta )}`$. Diffusion provides a natural drive toward a homogeneous distribution, and so it is conceivable that this faster mode of decay is manifested as a correction to the main behavior, $`ct^a`$. It follows from (6) that the correction exponent is

$$\mathrm{\Delta }=1\beta \alpha =\frac{1}{2n},$$
(7)
where the last equality applies to $`n`$-species annihilation, provided that the conjecture (5) holds. The more general relation works well for two-species annihilation with drift, where $`\alpha =1/3`$, $`\beta =7/12`$, and $`\mathrm{\Delta }=1/12`$ shafrir .

## IV Summary and discussion

We have presented large scale simulation results for diffusion-limited $`n`$-species annihilation, in one dimension, using the RRC method. Our simulations contradict previous work dba  and are in favor of new theoretical arguments advanced by Deloubrière et al. tauber . We have also provided a new scaling relation for the correction-to-scaling exponent $`\mathrm{\Delta }`$, valid for diffusion-limited reactions in one dimension, where the particles segregate into alternating domains. The corrections to the main decay mode are large, and explain the failure of dba  to obtain the correct asymptotic behavior with the size of simulations employed at that time. An important conclusion to be drawn is that predicting asymptotic behavior from the typical size of simulations used commonly in the field is dangerous. More advanced techniques and larger simulations seem to be imperative.

###### Acknowledgements.

We thank S. Redner for numerous illuminating discussions. We gratefully acknowledge the NSF (PHY-0140094) for partial support of this work.


text:
# Contents

## Chapter 0 Preface

This article reviews recent progress in magnetic light scattering in one- and two-dimensional quantum spin systems. These systems received considerable interest from both theoretical and experimental points of view. Following the investigations of the two-dimensional superconducting cuprates and the search for related transition metal oxides a fascinating field of copper oxide compounds, vanadates, manganites and nickelates opened up. These compounds show effects of strong electronic correlations and in particular magnetism in low dimensions.

The theory of magnetism in one dimension, on the other hand, has a history reaching back to the origin of quantum mechanics. This is due to the fact that a spin chain allows more easily analytical or numerical solutions. It was found that the suppression of “trivial” long-range order sets the stage for an enormous complexity of possible ground states, exotic quasiparticles and many-body states. Understanding these effects is the most intriguing challenge at present.

A central concept in describing these low-dimensional quantum spin systems is that of a spin liquid. This ground state is dominated by strong quantum fluctuations, pronounced spin-spin correlations and a suppression of long-range magnetic order. The Heisenberg chain with isotropic-antiferromagnetically coupled spins (s=1/2) represents such a state in the sense that the spin-spin correlations decay algebraically. It is therefore often denoted as a critical spin liquid. An interesting situation occurs when transitions lead to a sudden change of the excitation spectrum, e.g., the opening of an excitation gap or the formation of long-range magnetic order. These quantum phase transitions are driven or controlled by the exchange coupling parameters, the exchange topology or by spin vacancies. The excitation gap may be realized with or without a spontaneously broken translational symmetry.

The spin-Peierls transition and the related charge ordering instability discovered in the inorganic compounds $`\mathrm{CuGeO}_3`$ and $`\mathrm{NaV}_2\mathrm{O}_5`$, respectively, represent the case of broken translational symmetry. These compounds allow to investigate the excitation spectrum going from a homogeneous gapless to a dimerized state just as a function of temperature. In the two-leg spin ladder system $`\mathrm{SrCu}_2\mathrm{O}_3`$ and the chain/ladder system $`\mathrm{Sr}_{14}\mathrm{Cu}_{24}\mathrm{O}_{41}`$ an excitation gap is realized without breaking translational symmetry. These compounds are discussed as model systems for an electronic mechanism of high temperature superconductivity. The steady improvement of understanding also leads to surprising reinterpretations of compounds that have been investigated for years. The formerly canonical example of a spin ladder, the vanadium compound $`(\mathrm{VO})_2\mathrm{P}_2\mathrm{O}_7,`$ is now recognized as a spin chain with strongly alternating coupling constants. This result has profound consequences for the interpretation of its low energy excitations. Very important compounds that bridge one and two dimensions and still do not show long-range magnetic order are the spin frustrated system $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$ and the 1/5-depleted square lattice system $`\mathrm{CaV}_4\mathrm{O}_9`$.

Light scattering experiments or other spectroscopic methods like inelastic neutron scattering have been used to investigate both the above cited and many more compounds. One of the most significant aspects of light scattering experiments is the observation of magnetic singlet bound states. These states originate from strong triplet-triplet interaction and characterize the excitation spectrum of the spin system. Recent theoretical progress has enabled a more detailed understanding of these effects. Parameters like dimerization, frustration, interchain coupling, and spin-phonon coupling have an important impact on the ground state and the excitations of a quantum spin system.

This review is organized as follows: After a brief description of the excitations and the phase diagram of quantum spin systems given in Chapter 1, important low-dimensional spin systems, recent experimental results and their interpretations are discussed in Chapter 2. Up to now no comprehensive review on this rapidly growing field exists that also considers material aspects. Therefore we try to balance between well established results and very recent developments. In Chapter 3 magnetic light scattering in low-dimensional spin systems is reviewed. The following Chapters 4 and 5 discuss magnetic bound states and quasielastic scattering. Chapter 6 finally sums up some aspects of the present knowledge in this field and gives an outlook to future developments.

This review benefited from discussions and collaborations with many colleagues, too numerous to mention here completely.

We would like to thank especially E.Ya. Sherman and R. Valentí for carefully reading the manuscript. PL and GG would like to express their gratitude towards the recent and previous members of the light scattering group in Aachen, e.g., K.Y. Choi, J. Pommer, A. Ionescu, M. Fischer, M. Grove, G. Els, and P.H.M. van Loosdrecht.

Further valuable discussions with F. Mila, W. Brenig, K. Ueda, G.S. Uhrig, C. Pinettes, M. Udagawa, V. Gnezdilov and Yu. G. Pashkevich are gratefully acknowledged. We are obliged to C. Geibel, F. Steglich, H. Kageyama, M. Isobe, Yu. Ueda, J. Akimitsu, H. Tanaka, M. Johnsson, P. Millet, A. Revcolevschi, B.C. Sales, W. Assmus, S. Barilo and their collaborators for providing numerous samples. Furthermore, PL acknowledges gratefully a stipend of the Venture Business Laboratory for a visit at the Institute of Integrated Arts and Sciences (IIAS) Hiroshima University.

## Chapter 1 Excitations in low-dimensional spin systems

In strongly correlated electron systems with integer number of electrons per site the low energy excitations are usually given by the spin degrees of freedom. This situation is properly described by the Heisenberg exchange spin Hamiltonian. If, in addition, the exchange is restricted to low dimensions, then spin chains, spin ladders, and respective systems with a more complex exchange geometry are realized. These systems exhibit a number of unusual properties which are related to strong quantum fluctuations. These properties will be addressed in the following.

One-dimensional s=1/2 spin systems (spin chains) with uniform nearest neighbor exchange coupling show according to the Lieb-Schultz-Mattis theorem a degeneracy of the singlet ground state with triplet excitations. Assuming negligible spin anisotropies even for T=0 the ground state is gapless and not magnetically ordered. It is described by the Bethe Ansatz. The spin-spin correlations are algebraically decaying typical for a quantum critical state. Triplet excitations in this system are not described as magnons (bosons) but as massless domain wall-like s=1/2 spinons (fermions). These spinons are created as pairs, e.g., by an exchange process. Their dynamical structure factor is therefore given by a gapless two-particle continuum restricted by a lower and an upper dispersing boundary. In Fig. 1.1 a sketch of the spinon creation and the spinon continuum is given. The spectral weight of the continuum is dominant close to its lower boundary.

A quantum phase transition from a gapless critical state into a gapped state (disordered spin liquid) is induced by dimerization, i.e. an alternation $`\delta `$ of the coupling constants to nearest neighbors $`\mathrm{J}_{\mathrm{nn}}^\pm `$=(1$`\pm \delta `$)$`\mathrm{J}_{\mathrm{nn}}`$ along the chain or by a sufficient frustration $`\alpha `$=$`\mathrm{J}_{\mathrm{nnn}}`$/$`\mathrm{J}_{\mathrm{nn}}`$ due to next nearest neighbor antiferromagnetic exchange $`\mathrm{J}_{\mathrm{nnn}}`$. With dimerization the spinons are confined into massive triplet excitations. This confinement of spinon and antispinon composite objects (triplets) is discussed similar to the quark confinement in particle physics. The resulting quantum disordered ground state is characterized by short-ranged exponential decaying spin-spin correlations. In many cases the system is allowed to be described as an arrangement of spin dimers. The resulting lifted degeneracy of triplet and singlet excitations leads to an energy gain of the system. Fig. 1.2 shows a sketch of the excitation processes in a dimerized chain with the respective energy dispersion. The dimerization or alternation of the coupling constants connects k=0 and k=$`\pi `$ therefore a small part of the continuum spectral weight ($``$$`\delta `$2) is also expected for k$``$$`\pi `$/2. In Fig. 1.2 this contribution is neglected. A (T=0) phase diagram of dimerized and frustrated spin chains is given in Fig. 1.3. The points ($`\delta `$=0,$`\alpha `$=0) and ($`\delta `$=0,$`\alpha `$=0.5) correspond to the Bethe Ansatz and the Majumdar-Ghosh point. For $`\delta `$=0 and $`\alpha _c`$$`<`$0.3 the gap remains numerically small.

The simplest representative of the quantum disordered state, however, is the two-leg spin ladder with an approximately equal or larger exchange coupling along the rungs with respect to the coupling along the legs of the ladder. The singlet ground state is composed of spin dimers on the rungs. Here, the term spin liquid is even more appropriate as it is not based on a broken translational symmetry. An excitation in this picture of strong dimerization corresponds to the breaking of one dimer. The energy related to this process is the singlet-triplet gap $`\mathrm{\Delta }_{01}`$, see Fig. 1.2. A coupling of more than two chains to three-, four- or five-leg ladders leads to the experimentally proven conjecture that ladders with an even number of legs have a spin gap while odd-leg ladders are gapless . 222It should be mentioned that the combined effect of dimerization and interchain interaction may also lead to a vanishing spin gap in a quantum spin system. In a certain parameter space of a two-leg ladder with additional frustration, 4-spin cyclic exchange or interchain interactions, quantum phase transitions to gapless phases have been observed. On the other hand, there is theoretical evidence for a spin gap in doped three-leg ladders for a certain set of exchange coupling constants.

In the limit of a large number of coupled chains a two-dimensional Heisenberg system is obtained and the magnitude of the respective spin gap is going to zero. This limit may also be used to understand the two-dimensional high temperature superconductors (HTSC). Weakly doped two- and three-leg ladders have been theoretically investigated in this context. However, also in two dimensions spin dimer ground states with a gapped excitation spectrum are realized. This happens either if the exchange topology is modified to favor a dimer ground state, e.g., in removing 1/5 of the spins from a square lattice, or due to strong frustration (next nearest neighbor interaction).

The triplet-triplet interactions that are responsible for the opening of the gap also lead to magnetic bound states, i.e. triplet excitations that are bound to singlet, triplet or quintuplet states. The former two states are characterized by a well-defined excitation with an energy reduced with respect to the energy of the two-particle continuum of “free” triplets. If interchain or magnetoelastic interactions are dominant bound states consist of soliton-antisoliton pairs. Neglecting these effects the binding energy of a bound state in a dimerized spin chain originates from frustration. The maximum number of bound states of a spin chain is restricted to one singlet and one triplet state. In spin ladders with an additional diagonal frustration the number of bound states and their binding energy is less limited and increasing with frustration. In Fig. 1.4 the excitation spectrum of a homogeneous two-leg s=1/2 ladder is shown with a singlet and a triplet bound state at the lower boundary of the two-particle continuum. A quantum phase transition at a critical frustration into a gapless phase is understood as a condensation process of large many-particle bound states together with a general softening of the excitation spectrum. More generally, magnetic bound states may therefore be used to study the triplet-triplet interaction, determine the coupling parameters and the phase diagram of the system.

If defects, either as localized non-magnetic vacancies or as mobile carriers, are introduced into a quantum spin system its excitation spectrum may change drastically. In a 2D square lattice doping with mobile carriers destroys long-range Néel order and leads to the opening of a pseudo gap in the spin and charge excitations. The effect of localized spin vacancies in dimerized spin chains or in spin ladders is different. Here, a transition from a gapped into an ordered and gapless Néel-type state is induced by a seemingly negligible amount of vacancies. This effect is based on the doping of only weakly bound spinons by every vacancy. Thereby additional excitations are introduced in the gap corresponding to staggered moments for sites far from the vacancy and increasing spin-spin correlations. Interchain interaction leads to the occurrence of magnetic order at finite temperature.

The coexistence of true long-range magnetic order and dimerization is possible if spatial variations of the competing order parameters are taken into account. This means that the excitation spectrum of such a system has two features, the gapped triplet mode due to dimerization and the gapless “spin wave mode”. The emergence of antiferromagnetism keeps the structure and the energy scales of these modes essentially unchanged as the transfer of spectral weight from the gapped to the spin wave mode is realized with only a small reduction of the gap and an increase of the spin wave velocity. The latter mode is damped with a broadening proportional to the square of the wave vector. Similar arguments have recently been used to describe the interplay or competition of disorder-induced antiferromagnetism and superconductivity in Heavy Fermion compounds. Comprehensive experimental studies concerning the effect of spin vacancies exist for spin ladder and dimerized spin chain systems. Some results including light scattering data will be presented in the next chapter.

The effect of mobile carriers on the gapped excitation spectrum of a spin liquid is directly related to the problem of an electronic mechanism for HTSC and not yet understood completely. In the 2D CuO2 square lattice the doped holes are believed to form self-organized slowly fluctuating arrays of metallic stripes in which the motion of holes shows a locally quasi-one-dimensional character. A spin gap or pseudo gap is then the result of the spatially confined Mott-insulating regions of the material in the proximity of the metallic stripes. This effect has been described by the term “topological doping”. Corresponding theoretical studies of weakly doped two- and three-leg ladders confirm these ideas in the sense that a tendency toward a binding effect of holes either into a superconducting condensate or charge ordered ground state exists. The excitation spectrum of the latter system is of special interest as it is separated into a gapless Luttinger-liquid (odd channel) and an insulating gapped spin liquid phase (even channel). In some sense this spectrum represents or mimics the scenario of spin and charge separation discussed for HTSC. A quantum phase transition into a superconducting state with d-wave character has been predicted for the three-leg ladder at higher doping levels.

In a very simplified picture the two channels of the three-leg ladder may be understood as a “plain” ladder coupled to a chain. Holes hop back and forth from the chain to the ladder system. In the ladder they prefer to form pairs minimizing the number of broken dimers. Hopping back into the chain system this correlation is “partially transferred” into the conducting channel. The experimental part of this problem is far from being completely settled and therefore under progress. Although only few spin ladder systems are available unusual experimental results pointing to k-dependent relaxation rates different from the undoped material exist. In the following chapter a thorough review of the presently known inorganic low-dimensional compounds will be given focusing on the questions discussed above.

## Chapter 2 Important inorganic quantum spin systems

### 2.1 Structural Considerations and Important Parameters

### 2.2 Quasi-Zero-Dimensional Compounds

#### 2.2.1 Dimer compounds

#### 2.2.2 The Cu-Tetrahedra Systems $`\mathrm{Cu}_2\mathrm{Te}_2\mathrm{O}_5\mathrm{X}_2`$

### 2.3 Quasi-Two-Dimensional Compounds

#### 2.3.1 The 1/5-Depleted System $`\mathrm{𝐂𝐚𝐕}_\mathrm{𝟒}𝐎_\mathrm{𝟗}`$

#### 2.3.2 The Shastry-Sutherland System $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$

#### 2.3.3 High Temperature Superconductors

#### 2.3.4 The diluted square lattice system $`𝐊_\mathrm{𝟐}𝐕_\mathrm{𝟑}𝐎_\mathrm{𝟖}`$

### 2.4 Low-Dimensional Cuprates: new compounds related to high-temperature superconductors

#### 2.4.1 The Two-Leg Ladder $`\mathrm{SrCu}_2\mathrm{O}_3`$

#### 2.4.2 The Chain/Ladder System $`(\mathrm{Sr},\mathrm{Ca})_{14}\mathrm{Cu}_{24}\mathrm{O}_{41}`$

### 2.5 Low-Dimensional Vanadates

#### 2.5.1 The $`\mathrm{AV}_2\mathrm{O}_5`$ Family of Compounds

#### 2.5.2 The Alternating Chain System $`(\mathrm{VO})_2\mathrm{P}_2\mathrm{O}_7`$

### 2.6 Low-Dimensional Halides and Pnictides

#### 2.6.1 The Chain System $`\mathrm{𝐊𝐂𝐮𝐅}_\mathrm{𝟑}`$

#### 2.6.2 The $`\mathrm{ACuCl}_3`$ Family of Compounds

#### 2.6.3 The Haldane System $`\mathrm{𝐂𝐬𝐍𝐢𝐂𝐥}_\mathrm{𝟑}`$

#### 2.6.4 $`\mathrm{𝐘𝐛}_\mathrm{𝟒}\mathrm{𝐀𝐬}_\mathrm{𝟑}`$

### 2.7 Inorganic Systems with Spin-Peierls and related Instabilities

#### 2.7.1 $`\mathrm{𝐂𝐮𝐆𝐞𝐎}_\mathrm{𝟑}`$

#### 2.7.2 $`\mathrm{𝐍𝐚𝐕}_\mathrm{𝟐}𝐎_\mathrm{𝟓}`$

### 2.8 Magnetic Parameters of Selected Spin Systems

## Chapter 3 Magnetic light scattering

### 3.1 Two-Magnon Scattering

### 3.2 Light scattering in high temperature superconductors

### 3.3 Light Scattering in quasi-one-dimensional systems

#### 3.3.1 The Limit of Large Dimerization

#### 3.3.2 Nonzero Temperature

#### 3.3.3 Defect-Induced Light Scattering

### 3.4 Spinon Light Scattering in $`\mathrm{CuGeO}_3`$

### 3.5 Raman scattering from spin-1/2 ladders

## Chapter 4 Magnetic bound states

### 4.1 Bound states in $`\mathrm{CuGeO}_3`$

#### 4.1.1 Three-magnon scattering in $`\mathrm{CuGeO}_3`$

#### 4.1.2 Defect-induced bound states in $`\mathrm{CuGeO}_3`$

### 4.2 Bound States in $`\mathrm{NaV}_2\mathrm{O}_5`$

#### 4.2.1 Experimental observations in $`\mathrm{NaV}_2\mathrm{O}_5`$

#### 4.2.2 Deficiency and substitutions on the Na Site

#### 4.2.3 Theoretical considerations

### 4.3 Bound states in $`(\mathrm{VO})_2\mathrm{P}_2\mathrm{O}_7`$

### 4.4 Bound states in $`\mathrm{SrCu}_2(\mathrm{BO}_3)_2`$

#### 4.4.1 Effect of substitutions

#### 4.4.2 Polarization dependence

## Chapter 5 Quasielastic scattering in low-dimensional spin systems

## Chapter 6 Conclusions and outlook

## Acknowledgement


text:
# Three-Fluid Description of the Sympathetic Cooling of a Boson-Fermion Mixture

## I Introduction

Recent attempts to cool down a trapped gas of fermionic atoms, such as 40K  or 6Li , have reached temperatures below the Fermi temperature. Degenerate atomic Fermi gases offer the intriguing possibility to create a paired fermion state such as a Bardeen-Cooper-Schrieffer phase in a novel and highly controllable system. However, until now the temperatures obtained in the experiments with ultracold fermionic atoms are not low enough to observe such phenomena. In this paper, we propose a straightforward description of the cooling mechanism, discuss its inherent limitations, and we investigate strategies to alter the cooling procedure in order to achieve lower temperatures.

The cooling process used to create Bose-Einstein condensates in trapped bosonic gases is evaporative cooling, whereby the most energetic atoms are removed from the gas, leaving the remnant colder after rethermalization. However, due to the antisymmetrization requirement of the fermionic wave function, identical fermionic atoms cannot undergo the $`s`$-wave collisions necessary for rethermalization, and thus they cannot be directly cooled using evaporative cooling. This problem was circumvented by simultaneously trapping two atom species or two different hyperfine spin states of a given atom species. In , 6Li is trapped together with its bosonic isotope 7Li. Whereas the $`s`$-wave scattering length of 6Li is zero by symmetry, the $`s`$-wave scattering length for a collision between a 6Li and a 7Li atom is 2.2 nm, which is enough to thermalize both species together. This version of evaporative cooling is denominated ‘sympathetic cooling’.

We describe the sympathetic cooling process, using the three-fluid model presented in Sec. II, as a sequence of evaporation and rethermalization steps so that in the evaporation steps the energetic atoms are removed from the mixture and in the rethermalization steps the new (lower) equilibrium temperature is reached, as explained in Sec. III. The results obtained with this model are discussed in Sec. IV, and an improvement on the current sympathetic cooling process is proposed to reach lower temperatures ($`<0.1`$ $`T/T_F`$).

## II Three-fluid model

The trapped system which will be subjected to sympathetic cooling consists of $`N_b`$ bosons with mass $`m_b`$ and scattering length $`a_b`$; and $`N_f`$ fermions with mass $`m_f`$, all at initial temperature $`T`$ in a parabolic trap with trapping frequency $`\omega `$ and characteristic length $`a_{\text{HO}_{b,f}}=\sqrt{\mathrm{}/\left(m_{b,f}\omega \right)}`$. The scattering length for a collision between a boson and a fermion will be denoted by $`a_f`$. The three-fluid model  distinguishes the fermion gas, the Bose-Einstein condensed (BEC) bosons and the non-condensed or ‘thermal’ bosons. The number of BEC bosons is $`N_c=N_b(1(T/T_C)^3)`$ for $`T<T_C`$ and zero otherwise, where

$$T_C=\frac{\mathrm{}\omega }{k_B}\sqrt[3]{\frac{N_b}{\zeta \left(3\right)}}\left(10.73N_b^{1/3}1.33\frac{a_b}{a_{\text{HO}_b}}N_b^{1/6}\right).$$
(1)
The number of thermal bosons is $`N_t=N_bN_c`$. We use the notation

$$g=\frac{4\pi \mathrm{}^2a_b}{m_b}\text{ and }f=\frac{4\pi \mathrm{}^2a_f}{m_r}$$
(2)
where $`m_r=2m_bm_f/\left(m_b+m_f\right)`$ is twice the reduced mass of a fermion and a boson. $`n_c(𝐫)`$ is the density of BEC bosons, $`n_t(𝐫)`$ is the density of thermal bosons, and $`n_f(𝐫)`$ is the density of fermions. In the three-fluid model, the BEC bosons, the thermal bosons and the fermions are described as ideal gases in effective one-body potentials. The effective one-body potential for the BEC bosons ($`V_c`$), for the thermal bosons ($`V_t`$) and for the fermions ($`V_f`$) are given by:

$`V_c(𝐫)`$ $`=`$ $`{\displaystyle \frac{m_b\omega ^2}{2}}r^2+gn_c\left(𝐫\right)+2gn_t(𝐫)+fn_f(𝐫)`$ (3)
$`V_t(𝐫)`$ $`=`$ $`{\displaystyle \frac{m_b\omega ^2}{2}}r^2+2gn_c(𝐫)+2gn_t(𝐫)+fn_f(𝐫)`$ (4)
$`V_f(𝐫)`$ $`=`$ $`{\displaystyle \frac{m_f\omega ^2}{2}}r^2+fn_c(𝐫)+fn_t(𝐫)`$ (5)
These potentials determine the densities

$`n_c(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{g}}\left[\mu _b{\displaystyle \frac{m_b\omega ^2}{2}}r^22gn_t(𝐫)fn_f(𝐫)\right]`$ (6)
$`n_t(𝐫)`$ $`=`$ $`{\displaystyle \frac{d^3𝐤}{(2\pi )^3}\frac{1}{\mathrm{exp}\left[\beta \left({\displaystyle \frac{\mathrm{}^2k^2}{2m_b}}+V_t(𝐫)\mu _b\right)\right]1}}`$ (7)
$`n_f(𝐫)`$ $`=`$ $`{\displaystyle \frac{d^3𝐤}{(2\pi )^3}\frac{1}{1+\mathrm{exp}\left[\beta \left({\displaystyle \frac{\mathrm{}^2k^2}{2m_f}}+V_f(𝐫)\mu _f\right)\right]}}`$ (8)
In these expressions, $`\beta =1/(k_BT)`$ and $`\mu _f,\mu _b`$ is the chemical potential of the fermion gas and boson gas respectively. The set of equations (3)-(5) and (6)-(8) form a self-consistent description of the three fluids, in which the chemical potentials are set by fixing the total number of fermions and the total number of bosons; $`N_b=d^3𝐫(n_c(𝐫)+n_t(𝐫))`$ and $`N_f=d^3𝐫`$ $`n_f(𝐫).`$ Note that the interaction between the bosonic atoms, and the interaction between the bosonic and fermionic atoms, is taken into account on the level of the mean-field approximation.

The total energy of the BEC bosons ($`E_c`$), of the thermal bosons ($`E_t`$) and the fermions ($`E_f`$) is given by

$`E_c`$ $`=`$ $`{\displaystyle d^3𝐫\left[\frac{m_b\omega ^2}{2}r^2n_c\left(𝐫\right)+\frac{1}{2}g\left[n_c\left(𝐫\right)\right]^2\right]}`$ (9)
$`E_t`$ $`=`$ $`{\displaystyle d^3𝐫\left[\frac{d^3𝐤}{(2\pi )^3}\frac{\frac{m_b\omega ^2}{2}r^2+\frac{m_b}{2}k^2}{\mathrm{exp}\left[\beta \left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}+V_t(𝐫)\mu _b\right)\right]1}+g\left[n_t\left(𝐫\right)\right]^2+2gn_t\left(𝐫\right)n_c\left(𝐫\right)\right]}`$ (10)
$`E_f`$ $`=`$ $`{\displaystyle d^3𝐫\left[\frac{d^3𝐤}{(2\pi )^3}\frac{\frac{m_f\omega ^2}{2}r^2+\frac{m_f}{2}k^2}{\mathrm{exp}\left[\beta \left({\displaystyle \frac{\mathrm{}^2k^2}{2m_f}}+V_f(𝐫)\mu _f\right)\right]+1}+fn_f\left(𝐫\right)\left(n_c\left(𝐫\right)+n_t\left(𝐫\right)\right)\right]}.`$ (11)
From the above set of equations (3)-(11), the temperature, and the chemical potential of the bosons and the fermions can be calculated if the total numbers of bosons and fermions, and their total energy are given. Note that the Bose-Einstein condensate is treated in the Thomas-Fermi approximation as if it were at zero temperature. For the thermal bosons we have experimented using the Maxwell-Boltzman distribution instead of a Bose-Einstein distribution in the expressions for the density (7) and the energy (10). In the experiment, this approximation is routinely performed in order to determine the temperature of a trapped atomic gas by fitting the Maxwell-Boltzman distribution to the tail of the measured velocity distribution or density profile. Also in our calculations, we found that the Maxwell-Boltzman distribution provides an accurate description of the thermal bosons.

## III Sympathetic cooling

We describe the cooling process as a two-step process. The first step is the evaporative stage, whereby either just the bosons can be ‘evaporated’, or just the fermions, or both bosons and fermions. The effect of the evaporative stage is truncating the radial distribution by removing the outer atoms. The energy of these atoms is higher than the average energy of the mixture and so this removal allows us to obtain a lower temperature. The second step is the rethermalization stage, where the collisions between the bosons and between bosons and fermions bring the truncated (non-equilibrium) distribution obtained in the evaporative stage into a new equilibrium distribution. This new equilibrium distribution will be at a lower temperature if the atoms removed have more energy than the average atom energy. We will describe the sympathetic cooling as a discrete succession of evaporation and rethermalization stages. In reality, the evaporation and rethermalization goes on simultaneously. However, previous studies have shown that the discrete treatment of the cooling process provides an accurate description .

The evaporation is realized by a radiofrequent (or microwave) field which induces a transition between the trapped hyperfine state of a selected atom species and an untrapped hyperfine state. The frequency suited for this transition depends on the distance from the atom to the center of the trap, due to the Zeeman effect of the spatially varying magnetic trapping field. By tuning the frequency, atoms that are beyond a cut-off distance $`R`$ from the center of the trap are removed. If the bosons are evaporated in this way, then the number of bosons which are removed from the mixture is

$$\mathrm{\Delta }N_t=\underset{r>R}{}d^3𝐫n_t\left(𝐫\right),$$
(12)
and the total energy removed in the evaporation is

$$\mathrm{\Delta }E_t=\underset{r>R}{}d^3𝐫\left[\frac{d^3𝐤}{(2\pi )^3}\frac{\frac{m_b\omega ^2}{2}r^2+\frac{m_b}{2}k^2}{\mathrm{exp}\left[\beta \left({\displaystyle \frac{\mathrm{}^2k^2}{2m}}+V_t(𝐫)\mu _b\right)\right]1}+g\left[n_t\left(𝐫\right)\right]^2+fn_f\left(r\right)n_t\left(r\right)\right].$$
(13)
If fermionic atoms are removed from the mixture, then the corresponding number of removed fermions $`\mathrm{\Delta }N_f`$ and the total evaporated energy of the fermion gas $`\mathrm{\Delta }E_f`$ are found from expressions analogous to those above. Note that evaporating the Bose-Einstein condensed bosons will not lead to cooling. The Bose-Einstein condensed bosons are always the lowest-energetic bosons in the mixtures - removing them would only lead to heating of the system. Thus it is important that the cut-off distance $`R_C`$ is larger than the Thomas-Fermi radius of the condensate.

The evaporative stage results in a new total number of bosons $`N_b^{}=N_b\mathrm{\Delta }N_b`$, and a new total energy $`E_{\text{total}}^{}=E_c+E_t+E_f\mathrm{\Delta }E_t`$, if only bosons are removed from the trap. If also fermions are evaporated, $`\mathrm{\Delta }E_f`$ has to be subtracted too, and $`N_f^{}=N_f\mathrm{\Delta }N_f`$ is the new number of fermions. Using the equations (3)-(11), the new system parameters $`N_f^{},N_b^{},E_{\text{total}}^{}`$ are expressed as a function of the new temperature $`T^{}`$ and the new chemical potentials $`\mu _b^{}`$ and $`\mu _f^{}`$. The new temperature and chemical potentials are then solved from these equations. After this rethermalization stage, a new evaporative stage can be initiated.

Losses due to three-body collisions lead to heating, since the three-body collisions will be prominent in the condensate which harbors the coldest atoms. Such losses introduce another depletion of the number of atoms. For the bosons,

$`\mathrm{\Delta }N_b^{\text{heat}}`$ $`=`$ $`K_3t{\displaystyle 𝑑𝐫\text{ }[n_c(𝐫)]^3},`$ (14)
$`\mathrm{\Delta }E_b^{\text{heat}}`$ $`=`$ $`K_3t{\displaystyle 𝑑𝐫\text{ }[n_c(𝐫)]^3\left(\frac{m_b\omega ^2}{2}r^2+gn_c\left(𝐫\right)+2gn_t\left(𝐫\right)+fn_f\left(r\right)n_c\left(𝐫\right)\right)}.`$ (15)
where $`K_3`$ is the three-body collisional cross section and $`t`$ is the time duration of one step. For $`{}_{}{}^{7}Li`$, $`K_3`$ is about $`2.6\times 10^{28}cm^6/s`$ . This additional loss process can be taken into account during the evaporative stage.

## IV Results and discussion

### A Evaporating the bosonic atoms

In the experiments on sympathetic cooling by Schreck et al.  a mixture of bosons and fermions is cooled from 2 mK to around 9 $`\mu `$K, starting with $`2\times 10^8`$ 7Li bosons and about $`10^6`$ 6Li fermions, in an anisotropic trap ($`\omega _r=2\pi \times 4050`$ Hz, $`\omega _a=2\pi \times 71.17`$ Hz). We have simulated this process and compared with the experimental data . The suited number of steps can be determined by noting that the time needed for the system to rethermalize is 200 ms and that the reported cooling process takes 35 s . The number of 200 ms steps must be of the order of 100-200. For the evolution of the cut-off radius in time, we have used the same as in . The agreement between theory and experiment, shown in figure 1, is reasonably good.

In the next simulation we have started from $`9`$ $`\mu `$K, trying to cool further by evaporating only bosons, as is done in the experiment reported in Ref. . We simulate the evaporative cooling of a system, starting with $`N_f=4\times 10^3`$ fermions and $`N_b=10^6`$ bosonic atoms at an initial temperature of $`T=9`$ $`\mu K`$ which is $`6.12`$ times the Fermi temperature for this amount of fermions. Figure 2 shows the cooling process as a function of the time (number of evaporation-rethermalization steps), when the evaporation is performed at decreasing cut-off radius starting from $`20`$ $`a_{HO}`$ and decreasing with $`0.32`$ $`a_{HO}`$ each time step. The lowest temperature which is obtained is $`0.22`$ $`T_F`$, in agreement with the temperature reported by Schreck et al. . We have performed the simulation for different values of the three body boson decay constant ($`K_3`$) and one can see that in our model there is only an effect on the number of bosons, not on the final temperature.

Figure 3 illustrates the influence of the initial cut-off radius on the cooling process, showing the temperature at step 10,20,30,40 and 50 of the cooling process, as a function of the initial cut-off radius. During each step the cut-off radius was decreased by $`0.32`$ $`a_{HO}`$. An optimal radius exists around 18 $`a_{HO}`$, which is small enough to remove enough bosons, but sufficiently large in order not to remove too many bosons in the initial steps and not being energy selective. The simulations have a maximum of 50 steps or stop when the cut-off radius becomes smaller than the condensate. In figure 4 we investigate how the temperature obtained by cooling is affected by reducing the cut-off radius with different amounts during the cooling process. At each time step, the cut-off radius is decreased by a size $`a`$. The different curves in the figure show the temperature as a function of time (step in the cooling process) at different values of $`a`$. We find that a faster rate of decrease of the radius leads to a faster cooling, but that no lower temperatures can be reached. If the radius is decreased too fast, it reaches the radius of the condensate too soon, and the cooling process is stopped. For a fixed time of the cooling process, which is determined by the loss rate of the trap, an optimal value of $`a`$ can be determined. For the 50 steps that we took in figure 4, this is $`a`$ between 0.32 and 0.6 trap oscillator lengths. Figure 5 shows the dependence of the temperature on the initial number of bosons, keeping the number of fermions fixed. We find the counterintuitive result that the initial ratio between bosons and fermions has only a minor effect on the final temperature that can be reached: it is always around 0.2 $`T_F`$. This can be understood by noticing that the number of thermal bosons (which are responsible for the cooling) is almost independent of the total number of bosons for temperatures below the condensation temperature. ($`N_t=N_b\left(\frac{T}{T_c}\right)^3T^3\zeta \left(3\right)`$).

Truscott et al.  note that the cooling process seems to stall as soon as the specific heat of the bosons is smaller than the specific heat of the fermions, because the energy which can be carried away by evaporating a bosonic atom then becomes smaller than the energy of a fermionic atom. This result is confirmed by the present simulation, which explains why the final temperature is insensitive to the initial ratio of bosons to fermions. Basically, when the number of thermal bosons is small compared to $`N_F`$, the energy which can be carried away by evaporating the bosons is small compared to the energy of the Fermi gas. The temperature does not decrease significantly and the bosons are quickly lost, because the spatial extent of the cloud remains big.

### B An improved cooling process: sequential evaporation

This leads us to suggest another cooling process. When evaporative cooling by removing bosons has led to a final temperature, and continued evaporation of the bosons no longer leads to substantial further cooling, we suggest to switch to evaporating the fermions, or fermions and bosons. Removing the bosons from the trap does not cool the system any more because the bosons can no longer carry away the necessary energy - but the fermions still can.

Assume that by the evaporation of bosons, the temperature has already decreased to $`0.3`$ $`T_F,`$ and there are $`N_b=5\times 10^4`$ and $`N_f=5\times 10^4`$ atoms left in the trap. This system is experimentally achievable by evaporating bosons alone. Figure 6 shows the resulting temperature as a function of time, after switching from evaporating the bosons to evaporating only the fermions. The cut-off radius was chosen at $`R_C=15`$ $`a_{HO}`$ and decreases each time step. In this way the temperature can be lowered to $`0.06`$ $`T_F.`$ As also illustrated in the figure, the number of bosons has almost no influence on the cooling process.

Losses due to three particle collisions involving two bosons and one fermion is expected to limit the achievable temperature, as pointed out by Timmermans . The effect of this type of loss can be investigated in the present formalism in the same way as the effect of the three-boson collisions. The results are shown in the inset of figure 6, for a cooling process starting with $`5\times 10^4`$ bosons and $`5\times 10^4`$ fermions. The different curves correspond to different values of the three particle recombination constant $`K_{3BF}`$. If the product of the time to achieve equilibrium ($`t`$) and the recombination constant is $`K_{3BF}t10^7`$ we find that the lowest achievable temperature is raised by less than $`0.05`$ $`T_F`$. For $`K_{3BF}t10^8`$ there is no appreciable effect on the final temperature of the cooling process.

## V Conclusions

We have proposed a versatile, and straightforwardly implementable three-fluid model to describe the sympathetic cooling of bosons and fermions. Our model reproduces well the cooling observed in an experiment with 6Li and 7Li, and we find that the final temperature $`T/T_F0.2`$ which is reached is not sensitive to the initial ratio of the number of bosons to the number of fermions. Lower temperatures can be reached if the stage where the only the bosonic atoms are evaporated is followed by a stage where also the fermionic atoms are evaporated. When this procedure is applied, a temperature of $`T/T_F=0.06`$ could be reached, which in combination with the tunability of the interaction strength through the Feshbach resonance, brings the Fermi gas closer to temperatures where pairing may be observed.

## VI Acknowledgments

Two of the authors (M. W. and J. T.) are supported financially by the Fund for Scientific Research - Flanders (Fonds voor Wetenschappelijk Onderzoek – Vlaanderen). This research has been supported by the GOA BOF UA 2000, IUAP, the FWO-V projects Nos. G.0071.98, G.0306.00, G.0274.01, WOG WO.025.99N (Belgium).