text:
# Little-Parks effect in a superconducting loop with magnetic dot

## Abstract

We have studied the nucleation of superconductivity in a mesoscopic Al loop, enclosing magnetic dot with the perpendicular magnetization. The superconducting phase boundary $`T_c(B)`$, determined from transport measurements, is asymmetric with respect to the polarity of an applied magnetic field. The maximum critical temperature has been found for a finite applied magnetic field, which is antiparallel to the magnetization of the dot. Theoretical phase boundary shows a good agreement with the experimental data.

A great deal of experimental and theoretical work has been devoted to the investigation of mesoscopic superconductors (e.g. ref. victor ; vict1  and references therein). Over the past few years hybrid superconductor/ferromagnet systems have attracted a lot of attention ja ; pokro ; rusi ; misko ; miskoa ; miskob ; danac ; erdin ; mjvb ; martin . It has been revealed that a magnetic dot, embedded into a superconductor, attracts a vortex, resulting in an enhanced pinning when the magnetic moment $`m`$ of magnetic dot and an applied magnetic field $`B`$ have parallel orientation miskoa ; mjvb ; mjvbb . Moreover, it was recently demonstrated that perpendicularly magnetized magnetic dots could efficiently be used not only to enhance the flux pinning, but also to tune the nucleation of superconductivity martin .

Theoretical studies of individual mesoscopic superconductor/ferromagnet structures in the disk geometry were carried out in the framework of the Ginzburg-Landau formalism, with the emphasis on the vortex phases and transitions deep in the superconducting state misko ; miskob . A mesoscopic superconducting Al disk with a perpendicularly magnetized magnetic dot was recently fabricated and its superconducting phase boundary $`T_c(B)`$ was measured ja . These measurements demonstrated that the superconducting critical temperature is systematically higher when the applied magnetic field is parallel to the magnetization of the dot. However, in this sample the superconducting phase boundary seems to have been affected by the proximity effects between the superconducting disk and magnetic dot.

In this paper the nucleation of superconductivity in a mesoscopic superconducting Al loop, with a magnetic dot at the centre of the opening, is investigated. Given the spatial separation of the superconducting loop and magnetic dot, it is clear that the dot interacts with the loop exclusively through the stray field and proximity effects are absent. Therefore, any peculiar features in the $`T_c(B)`$ phase boundary can undoubtedly be attributed to the magnetic interaction between the superconducting condensate and the stray field of the dot.

The sample was prepared by electron beam lithography in three steps on a SiO2 substrate using positive PMMA950K and the corresponding co-PMMA electron beam resists. The lithography was carried out in a JEOL 5600 scanning electron micrograph modified for electron beam writing. The contacts pads, with the size of $`500\mu \mathrm{m}`$, leads and alignment markers were patterned first and $`5\mathrm{nm}`$ of Cr and $`30\mathrm{nm}`$ of Au were thermally evaporated. Afterwards, the resist was removed with a lift-off procedure, using warm acetone and ultrasonic agitation. The sample was covered again with the e-beam resist. After the alignment, the mesoscopic loop and contacts were patterned. The designed width of the loop and contacts is $`400\mathrm{nm}`$, whereas the inner radius of the loop is $`600\mathrm{nm}`$. In order to avoid proximity effects, the mesoscopic contacts next to the loop are $`22\mu \mathrm{m}`$ apart from the Cr/Au leads. $`35\mathrm{nm}`$ of Al was thermally evaporated and a lift-off was carried out. In the final step, the sample was yet again covered with the e-beam resists. After a quite precise alignment procedure, magnetic dot with the radius of $`200\mathrm{nm}`$ was patterned. A magnetic dot consisting of $`2.5\mathrm{nm}`$ Pd buffer layer and 10 bilayers of $`0.4\mathrm{nm}`$ Co and $`1\mathrm{nm}`$ Pd was subsequently evaporated. The resist was removed in a lift-off.

Fig. 1 shows an electron micrograph of the structure. The bright area at the centre is magnetic dot. The actual dimensions of the loop and magnetic dot were obtained from atomic force microscopy. The inner radius of the loop equals $`R_i=0.55\mu \mathrm{m}`$, outer radius is equal to $`R_o=1.05\mu \mathrm{m}`$, whereas the Co/Pd magnetic dot has the radius of $`r_d=0.27\mu \mathrm{m}`$. The mesoscopic contacts next to the loop, used for measuring voltage, are closer to the loop than initially designed. This effectively increases the outer radius of the loop, since the supercurrent can spread into the contacts, thereby affecting both the periodicity and background of the phase boundary victor .

The superconducting phase boundary $`T_c(B)`$ was obtained by four-point transport measurements in a cryogenic setup at temperatures down to $`1.11`$K, applying the magnetic field perpendicularly to the sample surface. A transport current with the rms of $`150`$nA and frequency $`27.7`$Hz was used. The response was measured by a PAR124A lock-in amplifier. The phase boundary was determined resistively by sweeping the magnetic field at a very low rate, whilst keeping the temperature constant. $`T_c(B)`$ curve was afterwards extracted from the magnetoresistance using the resistive criterion $`R_n/2`$, where $`R_n`$ is the resistance in the normal state. The temperature and field steps were $`500\mu \mathrm{K}`$ and $`5\mu \mathrm{T}`$, respectively, with the temperature stability $`100\mu \mathrm{K}`$. In order to minimize the influence of possible temperature gradients, the temperature was being monitored by two independent temperature sensors placed at two different positions on the sample holder.

The resistance of the sample at room temperature is $`6.6\mathrm{\Omega }`$, the resistance in the normal state at low temperatures is $`R_n=2.5\mathrm{\Omega }`$, whereas the maximum critical temperature $`T_{cm}`$ is $`1.3422`$K.

The multilayer Co/Pd magnetic dot provides a perpendicular magnetization martin ; gar .

The hysteresis loop in the perpendicular magnetic field of the co-evaporated plane film with the same composition, obtained by the magneto-optical Kerr measurements at room temperature, is presented in Fig. 2. The shaded area indicates the range in which the applied magnetic field was varied in the experiment. The remanence is complete, whereas the coercive field is approximately $`70\mathrm{mT}`$. When saturated, Co/Pd structures patterned at the sub-micrometre scale are in the single-domain state martin .

Fig. 3 shows the calculated stray field produced by a Co/Pd magnetic dot, averaged over the thickness of the superconducting loop. The stray field is normalized to the magnetization of the magnetic dot. The darker shaded area shows the position of the magnetic dot, whereas the lighter shaded areas indicate the position of the superconducting loop. The arrow shows the direction of the magnetization $`m`$ of the dot. The stray field was calculated assuming a single-domain state of the magnetic dot. Prior to the measurements the dot was magnetized perpendicularly in $`300`$mT. As the applied magnetic field in our experiment never exceeded $`30`$mT (Fig. 2), it has been assumed that the magnetization of the magnetic dot remains unaltered during the measurements martin .

The fluxes of the stray field (see Fig. 3) through the opening of the loop ($`r<0.55\mu \mathrm{m}`$) and loop itself ($`0.55<r<1.05[\mu \mathrm{m}]`$) have different orientations and, in qualitative terms, there should be a competition between the effects of these two fluxes on the superconducting condensate. Assuming that the loop is sufficiently big compared to the coherence length at zero temperature $`\xi (0)`$ ( $`R_i,R_o>>\xi (0)`$, where $`R_i`$ and $`R_o`$ are the inner and outer radius of the loop, respectively) the loop is similar to a stripe at temperatures much lower than the maximum critical temperature the system can attain $`T_{cm}`$ (we will be referring to these temperatures as low temperatures). Clearly, the superconducting state of the stripe depends only upon the flux value within the stripe. In other words, at low temperatures the behaviour of the loop is governed by the flux through the loop and not by the flux in the opening of the loop. Therefore, at low temperatures the shift of the $`T_c(B)`$ phase boundary along the $`B`$-axis is expected to be in the direction of the dot magnetization. At temperatures close to the maximum critical temperature (we will be referring to these temperatures as high temperatures), the direction of the shift of the $`T_c(B)`$ phase boundary is a result of an interplay between the opposite fluxes. Note that the superconducting phase boundary does not depend on the local profile of the stray field in the opening of the loop, but on the total flux through the opening generated by magnetic dot.

To describe quantitatively the $`T_c(B)`$ phase boundary the Ginzburg-Landau theory has been used. The free energy of the superconducting structure can be expressed as

$`F=2\pi {\displaystyle _{R_i}^{R_0}}r𝑑r\left[\alpha \psi +\beta \psi ^3+{\displaystyle \frac{\mathrm{}^2}{4m}}\left(L{\displaystyle \frac{2\pi }{\mathrm{\Phi }_0}}A\right)\psi ^2+{\displaystyle \frac{\mathrm{}^2}{4m}}\left({\displaystyle \frac{d\psi }{dr}}\right)^2\right]`$ (1)
where $`\psi `$ is the modulus of the order parameter, $`A`$ is the vector potential, an integer number $`L`$ stands the winding number (vorticity) of the order parameter, $`\mathrm{\Phi }_0`$ is the superconducting flux quantum, whereas $`\alpha `$ and $`\beta `$ represent the Ginzburg-Landau parameters. The integration is performed over the area of the loop. In Eq. (1) cylindrical symmetry of the order parameter near the phase transition to the normal state has been taken into account. As the thickness of sample is much smaller than $`\xi (T)`$ within the temperature range of interest, there is no modulation of the order parameter in the direction of the applied field. For this reason, the stray field can be averaged out over the thickness of the sample and the problem is reduced to the 1D case. Near the superconducting phase boundary the total magnetic field is equal to the sum of the stray field and an applied magnetic field. Therefore, the vector potential may be expressed as $`A(r)=A_s(r)+Br/2`$, where $`A_s(r)`$ is the vector potential of the stray field and $`B`$ is an applied magnetic field. We represent $`\psi (r)`$ as a Fourier expansion within the range $`R_irR_o`$ and consider the coefficients of the expansion the variational parameters. The values of these parameters can be calculated by a minimization of the free energy, yielding the solution of the full Ginzburg-Landau system of equations. The transition to the normal state corresponds to $`\psi =0`$. We have found that for our sample and experimental conditions, only the first term in the Fourier expansion is important, that is, $`\psi `$ is constant within the sample, whereas all other terms are negligibly small and practically do not affect the position of the phase boundary.

Fig. 4 presents the superconducting phase boundary, given as $`1T_c(B)/T_{cm}`$ versus the normalized flux $`\mathrm{\Phi }/\mathrm{\Phi }_0`$. $`T_{cm}`$ is the maximum critical temperature and $`\mathrm{\Phi }_0`$ is the superconducting flux quantum. Applied fields parallel to the magnetization of the dot are taken as positive and vice versa. The open symbols are experimental data, whereas the solid line corresponds to the theoretical curve for $`R_i=0.55\mu m`$, and $`r_d=0.27\mu m`$. The coherence length at zero temperature $`\xi (T=0)`$ has been used as a fitting parameter as it depends upon the preparation procedure and can show moderate variations even amongst samples prepared simultaneously. Due to the influence of the mesoscopic contacts, the outer radius has to be used as a fitting parameter, since the supercurrent is not confined only by the ring, but can well flow into the mesoscopic contacts, thus effectively increasing the outer radius. The best fit has been obtained for the coherence length at zero temperature $`\xi (0)=96`$nm and the outer radius $`R_o=1130`$nm. The latter is quite close to the actual radius $`R_o=1050`$nm.

The theoretical $`T_c(B)`$ curve is in a good agreement with the experimental data. The asymmetry, shift and periodicity of the experimental data have been reproduced well. The asymmetry and the shift can be explained by considering Eq. (1). From this equation it follows that the maximum critical temperature, for a fixed $`L`$, is attained when the applied field is equal to

$`B={\displaystyle \frac{4\mathrm{\Phi }_0}{\pi (R_o^4R_i^4)}}{\displaystyle _{R_i}^{R_o}}r^2𝑑r\left({\displaystyle \frac{2\pi A}{\mathrm{\Phi }_0}}{\displaystyle \frac{L}{r}}\right)`$ (2)
The corresponding critical temperature is

$$1\frac{T_c}{T_c\mathrm{}}=\frac{2\xi (0)^2}{R_o^2R_i^2}\left[_{R_i}^{R_o}r𝑑r\left(\frac{2\pi }{\mathrm{\Phi }_0}A\frac{L}{r}\right)^2\frac{4}{R_o^4R_i^4}\left(_{R_i}^{R_o}r^2𝑑r\left(\frac{2\pi }{\mathrm{\Phi }_0}A\frac{L}{r}\right)\right)^2\right]$$
(3)
where $`T_c\mathrm{}`$ stands for the bulk critical temperature of Al. For the sake of simplicity, we assume that the stray field does not change considerably over the width of the loop. Denoting the average field in the opening of the loop and within the loop as $`B_{op}`$ and $`B_l`$, respectively Eqs. (2) and (3) reduce to

$$B=B_l\frac{2\mathrm{\Phi }_0}{\pi (R_0^2+R_i^2)}\left(\frac{\pi R_i^2(B_lB_{op})}{\mathrm{\Phi }_0}L\right)$$
(4)
and

$$1\frac{T_c}{T_c\mathrm{}}=\frac{2\xi (0)^2}{R_o^2R_i^2}\left(\frac{\pi R_i^2(B_lB_{op})}{\mathrm{\Phi }_0}L\right)^2\left[\mathrm{ln}(\frac{R_0}{R_i})\frac{R_o^2R_i^2}{R_o^2+R_i^2}\right]$$
(5)
The maximum critical temperature $`T_{cm}`$ can be found by the minimization of the right-hand side of the Eq. (5) with respect to the integer $`L`$. Eq. (4) gives the shift of the $`T_c(B)`$ curve along the $`B`$-axis close to the $`T_{cm}`$, caused by the presence of an additional magnetic field. This shift is a sum of two contributions. The first term is just a compensation of the stray field $`B_l`$ within the loop. The second term describes the effect of the inhomogeneity of the stay field. This contribution behaves rather nontrivially due to the discreteness of $`L`$. For instance, when $`0<\pi R_i^2(B_{op}B_l)<0.5\mathrm{\Phi }_0`$, it follows from Eq. (5) that $`L=0`$, and the second contribution to the shift is negative. If the size of the dot and/or the intensity of the stray field are increased $`0.5\mathrm{\Phi }_0<\pi R_i^2(B_{op}B_l)<\mathrm{\Phi }_0`$, $`L`$ switches to $`1`$, and the second contribution becomes positive. Physically, this means that it is energetically more favorable for the system to change the vorticity by $`1`$ in order to decrease the modulus of the current in the loop. As a result, the current and the total shift, defined by Eq. (4), can change the sign. Thus, we may conclude, that the contribution to the shift, which is caused by the stray field inhomogeneity, is a periodic function of $`B_{op}B_l`$ and vanishes in the case of a homogeneous field, $`B_{op}=B_l`$. The total shift of the $`T_c(B)`$ curve along the $`B`$-axis close to the $`T_{cm}`$ can be of both signs for the same orientation of magnetization of the dot.

Our calculations have shown that at low temperatures the shift of the phase boundary is controlled only by the average stray field within the loop and is equal to -$`B_l`$, which is in accordance with the preliminary qualitative analysis. If the additional magnetic field were homogeneous, the shifts of the $`T_c(B)`$ phase boundary for low and high temperatures would have the same direction. The inhomogeneity of the stay field gives rise to two different and mutually independent shifts of the phase boundary in the low and high temperature ranges, which, in turn, twist the whole superconducting phase boundary, as presented in Fig. 4.

We would like to note that the discrepancy between the theoretical and the experimental data for low temperatures can be attributed to the effect of contacts, since they change the geometry of the sample, which cannot be entirely accounted for only by using the outer radius of the loop as a fitting parameter.

To summarize, we have fabricated a mesoscopic superconducting loop with a magnetic dot at the centre of the opening and investigated the onset of superconductivity in this structure by measuring the superconducting $`T_c(B)`$ phase boundary. The phase boundary has been found to be asymmetric with respect to the polarity of an applied magnetic field. The structure exhibits the maximum critical temperature for a finite value of the applied magnetic field, which is antiparallel to the magnetization of the magnetic dot. The theoretical superconducting phase boundary, obtained in the framework of the Ginzburg-Landau theory, is in a good agreement with the experimental data. It has been demonstrated that the inhomogeneous stray field shifts the $`T_c(B)`$ phase boundary along the $`B`$-axis independently in the low and high temperature ranges. In the limit of a thin loop with a large radius, these shifts are determined by the average value of the stray field in the opening of the loop and within the loop.

The authors would like to thank G. Rens for AFM measurements. This work has been supported by the Belgian IUAP, the Flemish FWO and the Research Fund K. U. Leuven GOA/2004/02 programmes, as well as by the ESF programme ”VORTEX”. W. V. P. acknowledges the support from the Research Council of the K.U. Leuven and DWTC.


format: text

text:
# Lifetime Measurement of the ³𝑃₂ Metastable State of Strontium Atoms

## Abstract

We have measured the lifetime of the $`5s5p{}_{}{}^{3}P_{2}^{}`$ metastable state of strontium atoms by magneto-optically trapping the decayed atoms to the ground state, which allowed sensitive detection of the rare decay events. We found that the blackbody radiation-induced decay was the dominant decay channel for the state at $`T=300`$ K. The lifetime was determined to be $`500_{130}^{+280}`$ s in the limit of zero temperature.

preprint: APS/123-QED

Precision measurements require the isolation of the physical system under study from environmental perturbations. The long coherence time thus obtained can be exploited to carry out ultra-high resolution spectroscopy Bergquist2000  or to create and control the macroscopic quantum coherence in atomic systems, such as Bose-Einstein Condensate (BEC) BECNature  and Cooper pairing BCS . Furthermore, the origin of environmental decoherence has recently attracted much attention in the context of quantum computation/communication, in which an entangled quantum system should maintain its internal and/or motional state coherence QC ; Bloch . The lower-lying $`{}_{}{}^{3}P`$ metastable states KatoriICAP ; Nagel ; Hemmerich ; Loftus ; ErtmerMg ; Udem ; Wilpers  of alkaline earth species are intriguing candidates for these studies. A long metastable lifetime may allow an optical spectroscopy at the 1 mHz level Katori\_Scottland , enabling one to realize an ultra-precise atomic clock. In addition, the possibilities of evaporatively cooling them to reach BEC are discussed theoretically SrBECDerevianko ; SrBECKokoouline .

Laser-cooled and trapped atoms in ultra-high vacuum condition offer an ideal sample for these studies, as the collisional interactions with container walls or residual gases can be substantially removed. However, in some cases, the radiation from the surrounding walls dramatically affects the evolution of the internal state coherences.The room-temperature blackbody radiation (BBR) has its intensity peak around a wavelength of $`10\mu `$m. Such BBR manifests itself most in Rydberg states Gallagher ; Farley\_Wing ; Hollberg , where infrared transitions with large electric dipole moments can be found. However, because of large energy differences, it hardly excites atoms in the lower-lying states. Its influence, therefore, has been rarely discussed in laser cooling and trapping. On the other hand, it is known that the BBR causes small but non-negligible ac Stark shifts in ultra-precise spectroscopy Itano ; Bauch .

As for the $`{}_{}{}^{3}P`$ metastable state of heavier alkaline-earth atoms, the upper-lying $`{}_{}{}^{3}D`$ states are connected to the metastable state by electric dipole transitions with mid-infrared wavelength. Therefore, the excitation of the metastable state by the BBR may significantly alter its effective lifetime or introduce detectable blackbody shifts in precision spectroscopy. In this Letter, we study the influence of the room-temperature BBR on the metastable state lifetime of Sr, which is recently predicted to be 1050 s Derevianko . Because the lifetime is significantly longer than the collision-limited lifetime of tens of seconds that is realized in neutral atom traps at a vacuum pressure of $`10^{10}\mathrm{torr}`$ Katori\_Shimizu , observing the survival of metastable atoms in such traps will not be practical Katori\_Shimizu ; Ertmer . Instead of directly observing the decay of the metastable atoms by emitted photons Walhout , we monitored the occurrence of the rare decay events to the ground state by magneto-optically trapping the atoms on the $`{}_{}{}^{1}S_{0}^{}{}_{}{}^{1}P_{1}^{}`$ transition. In this way, we used the magneto-optical trap (MOT) as a photon amplifier with a gain of $`10^{10}`$ to detect the decay events with unit quantum efficiency Nagourney ; Bergquist .

Figure 1 shows the relevant energy levels for 88Sr. Three radiative decay channels from the metastable $`5s5p{}_{}{}^{3}P_{2}^{}`$ state have been identified by Derevianko Derevianko : 1) The M2 transition to the ground state, 2) Decay to the $`{}_{}{}^{3}P_{1}^{}`$ state through the M1, E2, and M3 transition, 3) Decay to the $`{}_{}{}^{3}P_{0}^{}`$ state through the E2 transition, where E$`n`$ or M$`n`$ stands for an electric or a magnetic $`2^n`$-pole transition. The theory suggests that the contribution of the channel 3) is only 0.1% Derevianko . Therefore, 99.9% of the population in the $`{}_{}{}^{3}P_{2}^{}`$ state finally relaxes to the $`{}_{}{}^{1}S_{0}^{}`$ ground state, as the atoms that decayed to the $`{}_{}{}^{3}P_{1}^{}`$ state further decay to the $`{}_{}{}^{1}S_{0}^{}`$ ground state in $`22\mu \mathrm{s}`$. On this basis, we use the $`{}_{}{}^{1}S_{0}^{}`$ state population to detect the decay of the $`{}_{}{}^{3}P_{2}^{}`$ state.

We measured the number of atoms $`N_S(t)`$ in the $`{}_{}{}^{1}S_{0}^{}`$ ground state by applying the MOT on the $`{}_{}{}^{1}S_{0}^{}{}_{}{}^{1}P_{1}^{}`$ transition. The change of the atom number in the MOT is given by the rate equation,

$$\frac{d}{dt}N_S(t)=\left(\gamma _p+\mathrm{\Gamma }_c\right)N_S(t)+\left(\gamma _r+\mathrm{\Gamma }_q\right)N_P(t).$$
(1)
The atoms in the ground state are collisionally lost at a rate $`\mathrm{\Gamma }_c(10\mathrm{s})^1`$ or leaked out of the MOT transition via the $`5s4d{}_{}{}^{1}D_{2}^{}`$ state with the typical rate of $`\gamma _p(20\mathrm{ms})^1`$. The atoms are supplied by the $`{}_{}{}^{3}P_{2}^{}`$ metastable state with its atom number of $`N_P(t)`$, which either radiatively decay ($`\gamma _r`$) or collisionally quench ($`\mathrm{\Gamma }_q`$) to the ground state, where $`\mathrm{\Gamma }_q\gamma _r`$ holds as discussed later. Ignoring both of the collisional loss terms of $`\mathrm{\Gamma }_c`$ and $`\mathrm{\Gamma }_q`$, the steady-state solution of Eq. (1) for $`t\gamma _p^1`$ is given by,

$$\gamma _r=\frac{N_S}{N_P}\gamma _p,$$
(2)
where we assumed $`N_P(t)`$ to be constant as the decay of the population in the $`{}_{}{}^{3}P_{2}^{}`$ state with the collisional decay rate of $`\mathrm{\Gamma }_m(10\mathrm{s})^1`$ is small enough in the time scale of interest. The metastable decay rate $`\gamma _r`$, therefore, can be determined by the ratio $`N_S/N_P`$ of population distributed in both of the states and the MOT decay rate $`\gamma _p`$. The number of atoms in the $`{}_{}{}^{1}S_{0}^{}`$ ground state can be derived by observing the MOT fluorescence intensity of $`I_S=\eta N_S`$, where $`\eta `$ is the photon counting rate per atom trapped in the MOT. Similarly, $`N_P`$ can be obtained by transferring the metastable state population into the ground state and measuring the MOT fluorescence intensity of $`I_P=\eta N_P`$. By the ratio of these fluorescence intensities, $`\gamma _r/\gamma _p(=I_S/I_P)`$ is accurately determined regardless of the coefficient $`\eta `$.

The apparatus for magneto-optically trapping strontium atoms is similar to that described in Ref. Katori . The MOT fluorescence was collected by a lens with a solid angle of $`10^5`$ and then sent to a photomultiplier tube (PMT). An interference filter was placed in front of the PMT to block except the 461 nm light. The output signal was then sent to a multi-channel scaler. The two transitions, $`{}_{}{}^{3}P_{2}^{}{}_{}{}^{3}S_{1}^{}`$ and $`{}_{}{}^{3}P_{0}^{}{}_{}{}^{3}S_{1}^{}`$ (Fig. 1), were used to pump the $`{}_{}{}^{3}P_{2}^{}`$ metastable state population into the $`{}_{}{}^{1}S_{0}^{}`$ ground state.

We first loaded the MOT from an atomic beam on the $`{}_{}{}^{1}S_{0}^{}{}_{}{}^{1}P_{1}^{}`$ transition for 0.3 s. In the mean time, the $`{}_{}{}^{3}P_{2}^{}`$ metastable state was populated via the weak branching decay to the $`{}_{}{}^{1}D_{2}^{}`$ state from the $`{}_{}{}^{1}P_{1}^{}`$ state, which is estimated to be $`10^5`$ branch . About $`10^7`$ atoms in the low-field-seeking state were thus accumulated in the magnetic trap, which is formed by the quadrupole magnetic field used for the MOT. The field gradient was 100 G/cm along its axis of symmetry. At the typical peak density of $`10^9\mathrm{cm}^3`$, two-body collisional loss rate is estimated to be much smaller than the collisional loss term $`\mathrm{\Gamma }_c`$ in Eq. (1), where we assumed an inelastic collisional loss rate of $`10^{11}\mathrm{cm}^3/\mathrm{s}`$ SrBECDerevianko ; SrBECKokoouline . We then turned off the MOT lasers and closed the shutter that blocked both the atomic beam and the thermal radiation from the oven. We waited for 0.3 s so that all atoms except magnetically trapped metastable atoms diffused out of the trap region. After that, we turned on the MOT lasers again to capture atoms that were radiatively decayed from the $`{}_{}{}^{3}P_{2}^{}`$ state and recorded the fluorescence intensity. At $`t=1\mathrm{s}`$, we irradiated both of the pumping lasers to transfer the metastable state population into the $`{}_{}{}^{1}S_{0}^{}`$ ground state and determined the number of atoms trapped in the $`{}_{}{}^{3}P_{2}^{}`$ state by the MOT fluorescence intensity.

Figure 2 shows a change of fluorescence intensity averaged over $`10^2`$ measurements, where the background level was subtracted by alternating the same procedure with and without loading atoms into the magnetic trap. The number of atoms in the $`{}_{}{}^{1}S_{0}^{}`$ state at $`t=1`$ s was determined by exponentially extrapolate the whole decay curve $`I_S(t)`$ of the MOT fluorescence to obtain $`I_S(1)`$ with better statistics. The fluorescence decay in $`0<t<1`$ s was mainly caused by the collisional atom loss in the magnetic trap with the decay rate $`\mathrm{\Gamma }_m=(6.6\mathrm{s})^1`$ at the background gas pressure of $`6.6\times 10^{10}\mathrm{torr}`$. By transferring the metastable state population into the ground state at $`t=1\mathrm{s}`$, the fluorescence intensity sharply rose up to its maximum in about 5 ms. We approximated $`I_P(1)`$ by the peak fluorescence intensity shown in the inset of Fig. 2. The signal then decayed double-exponentially, consisting of the MOT decay with $`\gamma _p=(14.1\mathrm{ms})^1`$ due to the branching loss and the much slower collisional decay $`\mathrm{\Gamma }_m`$ of the metastable atoms recaptured in the magnetic trap. The metastable state lifetime is calculated by applying the ratio $`N_S/N_P=I_S(1)/I_P(1)`$ and the measured MOT decay rate $`\gamma _p`$ in Eq. (2). The measurement shown in Fig. 2 gave an effective radiative lifetime of $`\gamma _r^1=104_7^{+8}`$ s, which is only one tenth of the theoretical lifetime Derevianko . This shortening can be attributed to the BBR-induced decay via the $`5s4d{}_{}{}^{3}D`$ state as discussed later. However, before discussing the BBR-induced decay, we checked the other decay channels.

The reduction of the lifetime may be caused by the fine structure mixing of the metastable $`{}_{}{}^{3}P_{2}^{}`$ state with the $`{}_{}{}^{3}P_{1}^{}`$ or the $`{}_{}{}^{1}P_{1}^{}`$ state in the presence of the trapping magnetic field. Assuming a magnetic field of 10 G, the magnetically induced decay rates are estimated to be 0.1% of the natural decay rate $`\gamma _0`$ of the $`{}_{}{}^{3}P_{2}^{}`$ state. Actually, we measured the metastable lifetime under various magnetic field gradient, however, the change of the lifetime was within the statistical errors of 5 % as expected. Second, due to collisions with background gases, atoms in the metastable state may be 1) kicked out of the magnetic trap with the rate $`\mathrm{\Gamma }_m`$ or 2) quenched into the ground state with $`\mathrm{\Gamma }_q`$. For the former issue, since we compared the fluorescence intensity $`I_S`$ and $`I_P`$ just before and after the population transfer, the collisional atom loss in the transferring period of $`\tau _p=5`$ ms may cause an error. However its fraction is estimated to be as small as $`1e^{\mathrm{\Gamma }_m\tau _p}10^3`$. For the latter issue, the collisional quench to the ground state may cause pressure-dependent shortening of the metastable state lifetime. To check this influence, we increased the background gas pressure up to $`1.6\times 10^9\mathrm{torr}`$ and measured the lifetime. However, the change in the decay rate was well within statistical uncertainties.

The rapid decay can be attributed to the metastable state quenching due to the BBR field that transfers atom population in the $`{}_{}{}^{3}P_{2}^{}`$ state to the short-lived $`{}_{}{}^{3}P_{1}^{}`$ state via the $`{}_{}{}^{3}D`$ states. By solving the coupled rate equations, the steady-state value of the BBR-induced decay rate $`\gamma _B(T)`$ is expressed as,

$$\gamma _B(T)=\gamma _D\frac{7}{36}\overline{n}(T),$$
(3)
where $`\gamma _D`$ is the radiative decay rate of the $`5s4d{}_{}{}^{3}D`$ state, and $`\overline{n}(T)=\{\mathrm{exp}(hc/k_BT\lambda )1\}^1`$ is the BBR photon occupation number at temperature $`T`$ with the wavelength $`\lambda 3.01\mu \mathrm{m}`$ for the $`{}_{}{}^{3}P_{2}^{}{}_{}{}^{3}D`$ transition. Using the radiative lifetime $`\gamma _D^1=2.9\pm 0.2`$ $`\mu `$s of the $`{}_{}{}^{3}D`$ states Miller , the BBR-induced decay rate is calculated to be $`\gamma _B(T_0)=8.03\times 10^3\mathrm{s}^1`$ for $`T_0=299.5`$ K. Therefore the intrinsic decay rate $`\gamma _0`$ of the $`{}_{}{}^{3}P_{2}^{}`$ state at $`T0`$K is determined as $`\gamma _0=\gamma _r(T_0)\gamma _B(T_0)=(2.3\pm 0.7)\times 10^3\mathrm{s}^1`$, which corresponds to a lifetime of $`440_{100}^{+200}`$ s. To confirm that this lifetime shortening originates in the BBR excitation, we measured the decay rate $`\gamma _r(T)`$ as a function of the ambient temperature $`T`$ by heating up the vacuum chamber that enclosed the MOT. We note that this temperature change altered the vacuum pressure in the range of $`(6.69.7)\times 10^{10}\mathrm{torr}`$. However, the resultant collisional losses did not affect the measured lifetime as mentioned previously. The measured decay rate $`\gamma _r(T)`$ is plotted in Fig. 3 by filled squares, where error bars indicate one standard deviation. The monotonic increase of the decay rate as the temperature clearly supports that this lifetime shortening is caused by the BBR.

The BBR-induced decay rate is estimated as follows. Initially, the chamber that surrounds the atoms is assumed to be in thermal equilibrium. However, in heating up the chamber, care should be taken if it is in thermal equilibrium. The vacuum chamber is made of stainless steel with 150 mm diameter and has six viewing ports made of BK7 glass with 40 mm diameter. It was heated by a ribbon heater wound around the chamber body. We monitored the temperature of the viewing ports (at their rim and center) and the body of the chamber one hour after changing the heating power: Whereas the temperature of the rim was the same as that of the chamber body $`T`$, that of the center was somewhat lower by $`\delta T(>0)`$. Assuming the linear change of the temperature between the rim and the center, the average temperature for the viewing ports was approximated as $`T_p=T\delta T/3`$. Owing to this temperature inhomogeneity, atoms are not in the BBR field of thermal equilibrium. Assuming a spherical radiant cavity, the effective photon occupation number for the body temperature of $`T`$ is calculated as Chandos ,

$$\overline{n}_{\mathrm{eff}}(T)=\frac{\mathrm{\Omega }_b\epsilon _b\overline{n}(T)+\mathrm{\Omega }_p\epsilon _p\overline{n}(T_p)}{\mathrm{\Omega }_b\epsilon _b+\mathrm{\Omega }_p\epsilon _p}.$$
(4)
Here $`\mathrm{\Omega }_b`$ is the solid angle covered by the chamber body and $`\epsilon _b`$ is its spectral emissivity at $`\lambda 3\mu `$m. Likewise we defined $`\mathrm{\Omega }_p`$ and $`\epsilon _p`$ for the viewing ports, which fraction in the solid angle was $`\mathrm{\Omega }_p/\mathrm{\Omega }_b10^1`$. We used this effective photon occupation number $`\overline{n}_{\mathrm{eff}}(T)`$ to calculate the BBR-induced decay rate given by Eq. (3). We took the spectral emissivity of both objects as fitting parameters so that the measured decay rate $`\gamma _r(T)`$ should have a constant offset $`\gamma _0`$ to the calculated BBR-induced decay rate $`\gamma _B(T)`$ at each temperature. After the least squares fitting, we obtain the metastable decay rate of $`\gamma _0=(2.1\pm 0.4)\times 10^3\mathrm{s}^1`$ and the spectral emissivities of $`\epsilon _b=1.6\%`$ and $`\epsilon _p=88\%`$ for stainless steel and BK7 glass, respectively. The obtained emissivity for BK7 glass showed fair agreement with $`\epsilon =96\%`$ derived from the index of refraction $`n=1.48`$ for BK7 glass at $`\lambda 3\mu `$m. Reference data for polished stainless steel, however, was not available.

In the above discussion, we have not included the uncertainty of the $`{}_{}{}^{3}D`$ state radiative lifetime, which gives coupling strength between the $`{}_{}{}^{3}P`$ and $`{}_{}{}^{3}D`$ states and thus significantly affects the BBR-induced decay rate. In Ref. Miller , the observed radiative lifetime of the $`{}_{}{}^{3}D`$ state is given as $`\tau _D=2.9\pm 0.2\mu \mathrm{s}`$. This uncertainty brings another statistical error of $`\pm 0.6\times 10^3\mathrm{s}^1`$ to the metastable decay rate. Thus, we finally obtain the metastable decay rate $`\gamma _0=(2.0\pm 0.7)\times 10^3\mathrm{s}^1`$, or the metastable lifetime of $`500_{130}^{+280}`$ s.

In summary, we have determined the lifetime of the $`5s5p{}_{}{}^{3}P_{2}^{}`$ state of $`{}_{}{}^{88}\mathrm{Sr}`$ to be $`500_{130}^{+280}`$ s in the limit of zero temperature. Because the room-temperature BBR considerably shortens the metastable lifetime, care should be taken when dealing with the state to form a BEC SrBECDerevianko ; SrBECKokoouline  or any other applications that require long coherence time. We, therefore, need to prepare cold environment to suppress thermal photons. For example, by just lowering the ambient temperature down to 275 K, the BBR-induced decay rate becomes comparable to its intrinsic decay rate. A straightforward comparison of the measured lifetime with the theory will be possible by further lowering the temperature down to 217 K, where the BBR quenching rate is expected to be 1% of the natural decay rate.

The authors thank K. Okamura and M. Takamoto for their technical support. This work was supported by the Grant-in-Aid for Scientific Research (B) (12440110) from the Japan Society for the Promotion of Science.


text:
# Limits of sympathetic cooling of fermions: The role of the heat capacity of the coolant

## I Introduction

Sympathetic cooling of fermionic atoms has proven to be an efficient tool to produce strongly degenerate Fermi gases truscott1 ; schreck1 ; roati2002 ; hadzibabic2002 ; hadzibabic2003 . It is often thought that when the specific heat of the coolant becomes less than that of the fermions the cooling process slows drastically and even comes to a halt truscott1 ; onofrio1 ; presilla2003 . This belief derives from the following argument based on standard thermodynamics huang1 . A coolant is not only described by its temperature but also by its heat capacity. The equilibrium temperature between two bodies initially at different temperatures $`T_f`$ and $`T_b`$ is $`T=(C_fT_f+C_bT_b)/(C_f+C_b)`$, where $`C_f`$ and $`C_b`$ are the heat capacities of the two bodies constance . Therefore, regardless of having a very cold object $`b`$ ($`T_bT_f`$), it will not act as a good coolant unless it is also true that $`C_bC_f`$.

However, this argument does not include the evaporation of atoms from the coolant hess1986 ; luiten1996 , a key ingredient of sympathetic cooling. In the following, we shall consider the ideal case of sympathetic cooling, where all losses of particles due to inelastic collisions are absent. One can then evaporate sufficiently slowly that the Fermi gas and the coolant remain close to thermal equilibrium. This is in the spirit of the above argument, based on heat capacities, as heat capacity is a property of an equilibrium system. In this case, the only limitation of resources is the number of particles in the coolant. We then introduce an idealized model of sympathetic cooling: the atomic motion is treated in the semi-classical approximation and cooling proceeds in discrete steps. In each step, the atoms of the coolant go through an evaporation phase, and then the fermions and the coolant thermalize (see Fig. 1 below). This is a generalization of the discrete scheme used for evaporative cooling in Ref. Mewes , to the case of sympathetic cooling. A similar semi-classical and discrete model for sympathetic cooling was put forward in Ref. wouters2002 , with the addition of features like collisional losses and space dependent effects in the evaporation process, but without a determination of the limits of the cooling process. Using our idealized model, we shall then show that the argument based on the heat capacities is approximately correct for fermions cooled by an ideal Boltzmann gas prudence , but incorrect when cooled by an ideal Bose gas below the Bose condensation temperature $`T_c`$. In fact, in the latter case, it is *preferential* to have $`C_bC_f`$ to minimize the number of particles of the coolant whose evaporation is required to reach zero temperature proviso .

## II The Model

We assume that the Fermi gas, denoted by the subscript $`f`$, and the coolant, denoted by the subscript $`b`$, are initially in thermal equilibrium at a common temperature $`T`$. One evaporates all particles of the coolant with an energy greater than $`\eta k_BT`$ Mewes ; hess1986 ; luiten1996 ; cct1997 , where $`\eta `$ is a constant, typically 5 to 7 in current experiments. This has the effect of removing a number of particles $`\delta N_b`$ and an energy $`\delta E_b`$ from the coolant:

$`\delta N_b={\displaystyle _{\eta k_BT}^{\mathrm{}}}𝑑ϵ\rho _b(ϵ)N_b(ϵ),`$ (1)
$`\delta E_b={\displaystyle _{\eta k_BT}^{\mathrm{}}}𝑑ϵ\rho _b(ϵ)N_b(ϵ)ϵ,`$ (2)
where $`\rho _b(ϵ)`$ is the density of states of a particle in $`b`$ and $`N_b(ϵ)`$ is the mean occupation number of a state of energy $`ϵ`$. In what follows we restrict ourselves to the semi-classical regime in which the number of trap levels with an energy less than $`k_BT`$ is much larger than unity, both for species $`f`$ and $`b`$. Then the case of a density of states of the form

$$\rho _{f,b}(ϵ)=A_{f,b}ϵ^{\delta +1/2}$$
(3)
may be considered: this encompasses both a harmonic potential $`(\delta =3/2)`$ and a box $`(\delta =0)`$ dos\_comment . Here the subscript $`f`$ ($`b`$) stands for the fermions (coolant). In the two limiting cases of a non-degenerate gas $`N_b(ϵ)=\mathrm{exp}[\beta (ϵ\mu )]`$ and a Bose gas below $`T_c`$, $`N_b(ϵ)=1/[\mathrm{exp}(\beta ϵ)1]`$, the mean energy of an evaporated particle cct1997  is independent of the chemical potential $`\mu `$ and proportional to $`T`$:

$$\frac{\delta E_b}{\delta N_b}=f(\eta )k_BT.$$
(4)
The function $`f(\eta )`$ is different for a Bose gas below $`T_c`$ and for a Boltzmann gas, but has the same large $`\eta `$ behaviour in both cases formules :

$$f(\eta )=\eta +1+\frac{\delta +\frac{1}{2}}{\eta }+𝒪\left(\frac{1}{\eta ^2}\right).$$
(5)
Once the particles are removed, the system is left to rethermalize. The decrease in temperature $`\delta T`$ is obtained by conservation of energy,

$`E_f(N_f,T)+E_b(N_b,T)\delta E_b`$
$`=E_f(N_f,T\delta T)+E_b(N_b\delta N_b,T\delta T).`$ (6)
Note that $`E_f`$ is the total energy of the fermions, not the Fermi energy. The parameter $`\eta `$ is chosen sufficiently large so that $`\delta N_bN_b`$ and $`\delta TT`$. Equation (6) may be expanded to leading order in $`\delta N_b`$ and $`\delta T`$:

$$\frac{\delta T}{\delta N_b}=\frac{f(\eta )k_BT{\displaystyle \frac{E_b}{N_b}}(N_b,T)}{C_f(N_f,T)+C_b(N_b,T)},$$
(7)
where the heat capacities are defined by $`C_f_TE_f`$ and $`C_b=_TE_b`$.

Equation (7), valid for a Boltzmann gas or a degenerate Bose gas, will be studied for both a three dimensional box and a harmonic potential. In the former case, the mean field energy per particle of the coolant, due to the interaction with the Fermi gas, is constant and therefore plays no role. In the latter case, it is not exactly constant, but we shall assume that it is negligible as compared to $`k_BT`$. Additionally, for simplicity it will be assumed that the initial state of the fermions is degenerate. In this case the total energy is

$`E_f(N_f,T)=E_f(N_f,0)`$
$`+\alpha N_fk_B{\displaystyle \frac{T^2}{T_F}}\left\{1+\kappa \left({\displaystyle \frac{T}{T_F}}\right)^2+𝒪\left[\left({\displaystyle \frac{T}{T_F}}\right)^4\right]\right\},`$ (8)
where $`\alpha \pi ^2(\delta +3/2)/6`$, $`\kappa \pi ^2(\delta +1/2)(\delta 3)/20`$, $`\delta `$ is the exponent in the density of states and the Fermi temperature $`T_F`$ is given by

$$A_f(k_BT_F)^{\delta +3/2}=(\delta +3/2)N_f.$$
(9)
In Eq. (8) we shall keep terms only to second order in $`T`$ in the following.

## III Cooling by an ideal Bose gas below $`T_c`$

A non-interacting Bose gas with a temperature below $`T_c`$ has a total energy

$$E_b(N_b,T)=\frac{\nu A_b}{\delta +5/2}(k_BT)^{\delta +5/2},$$
(10)
where $`\nu (\delta +5/2)\mathrm{\Gamma }(\delta +5/2)\zeta (\delta +5/2)`$ is a constant and $`\mathrm{\Gamma }`$ and $`\zeta `$ are gamma and zeta functions, respectively abramowitz1 . Note that $`E_b`$ does not depend on the total number of bosons, since the occupation of the excited states is saturated to its maximal value in the presence of a condensate. Therefore $`E_b/N_b=0`$ and Eq. (7) reduces to

$$\frac{\delta T}{\delta N_b}=\frac{f(\eta )T_F}{2\alpha N_f\left[1+{\displaystyle \frac{3\nu }{\pi ^2}}{\displaystyle \frac{A_b}{A_f}}\left({\displaystyle \frac{T}{T_F}}\right)^{\delta +1/2}\right]}$$
(11)
where we used the fact that $`A_f`$ is related to the Fermi temperature $`T_F`$ by Eq. (9). When species $`f`$ and $`b`$ are in a box, $`A_b/A_f`$ is simply the ratio $`V_b/V_f`$ of the spatial volumes occupied by the two species. When species $`f`$ and $`b`$ are trapped harmonically, $`A_b/A_f`$ is given by the ratio $`(\omega _f/\omega _b)^3`$ where $`\omega _{f,b}`$ are the geometric means of the three oscillation frequencies for species $`f`$ and $`b`$, respectively.

The term in $`T^{\delta +1/2}`$ in the denominator of Eq. (11) represents the ratio of the two heat capacities:

$$\frac{C_b}{C_f}=\frac{3\nu }{\pi ^2}\frac{A_b}{A_f}\left(\frac{T}{T_F}\right)^{\delta +1/2}.$$
(12)
If one measures the efficiency of the cooling by the temperature decrease per number of bosons evaporated, i.e., $`\delta T/\delta N_b`$, it is most efficient to have $`C_bC_f`$. In this case the first term in the denominator dominates and the derivative is maximized. In Fig. 2 is shown the temperature as a function of the number of bosons evaporated, as can be obtained from the solution of the differential equation Eq. (11) for the case of harmonic trapping and $`\eta =6`$. It is apparent that it is advantageous to evaporative cooling to choose a trap frequency ratio $`\omega _f/\omega _b1`$, in order to maintain $`C_b/C_f`$ small at typical experimental temperatures $`T/T_F0.2`$ or $`0.1`$.

A simple condition may be obtained from Eq. (11) for the total initial number of bosons necessary to cool to $`T=0`$:

$$N_{b,0}>\mathrm{\Delta }N_b=\frac{C_{f,0}}{k_Bf(\eta )}\left[1+\frac{1}{\delta +3/2}\frac{C_{b,0}}{C_{f,0}}\right],$$
(13)
where $`T_0`$ is the initial temperature, $`C_0`$ are the initial heat capacities and $`\mathrm{\Delta }N_b`$ is the number of bosons evaporated. As the coolant is assumed to be initially a Bose gas with a condensate present, the initial number of bosons must exceed the saturated value of the population of the excited states:

$$N_{b,0}>A_b\mathrm{\Gamma }(\delta +3/2)\zeta (\delta +3/2)(k_BT_0)^{\delta +3/2}.$$
(14)
Note that, as an arbitrarily large number of particles can be stored in the condensate, no connection can be established between the number of particles and the heat capacity of the coolant. Thus condition (13) cannot be considered as a constraint on $`C_{b,0}/C_{f,0}`$ implic .

Recall that $`\delta T/T1`$ was required to replace the finite differences in (6) by derivatives. When the heat capacity of the Bose gas is initially much smaller that that of the Fermi gas, one finds that this condition is satisfied for any value of $`\eta `$, including $`\eta =0^+`$ where $`f(0)`$ is finite and strictly positive. In this case, Eq. (11) and its consequences apply also when one evaporates all the non-condensed atoms at each cooling cycle.

## IV Cooling by a Boltzmann gas

A Boltzmann gas has a total energy

$$E_b(N_b,T)=(\delta +3/2)N_bk_BT.$$
(15)
Equation (7) then reduces to

$$\frac{\delta T}{\delta N_b}=\frac{[f(\eta )(\delta +3/2)]T}{2\alpha N_f{\displaystyle \frac{T}{T_F}}+(\delta +3/2)N_b}.$$
(16)
The explicit dependence of the energy on $`N_b`$ introduces a negative contribution in the numerator. However, the total numerator remains positive, since one can show, by integration by parts, that the condition

$$f(\eta )>\delta +3/2$$
(17)
is satisfied for all $`\eta >0`$.

One finds that, by making a change of variables

$$x\frac{2\alpha }{\delta +3/2}\frac{T}{T_F}\frac{N_f}{N_b}=\frac{C_f}{C_b},$$
(18)
the equation becomes separable. This leads to the constant of motion

$$K=\frac{\lambda }{\lambda 1}\mathrm{ln}|\lambda 1x|+\frac{1}{\lambda 1}\mathrm{ln}x\mathrm{ln}(N_b/N_f),$$
(19)
where

$$\lambda \frac{f(\eta )(\delta +3/2)}{\delta +3/2}>0.$$
(20)
The asymptotic behavior of the temperature in the limit of the complete evaporation of the Boltzmann particles, i.e., $`N_b0`$, falls into two regimes. The first is given by $`x_0>\lambda 1`$, where $`x_0`$ is the initial value of $`x`$. In this case, the temperature tends to a finite non-zero value due to an insufficient number of Boltzmann particles. The ratio of the final to initial temperature is

$$\frac{T_{\mathrm{final}}}{T_0}=\left[1\left(\lambda 1\right)\frac{C_{b,0}}{C_{f,0}}\right]^{\lambda /(\lambda 1)}.$$
(21)
The second regime is given by $`x_0<\lambda 1`$, which implies $`\lambda >1`$. The inequality $`x_0<\lambda 1`$ is equivalent to the requirement that there be a minimal initial number of Boltzmann particles:

$$N_{b,0}>\frac{C_{f,0}}{k_B[f(\eta )2(\delta +3/2)]}.$$
(22)
Note that $`\lambda >1`$ implies that the denominator is positive. Equation (22) may be equally expressed as a requirement on the initial heat capacities of the fermions and Boltzmann particles:

$$C_{f,0}<C_{b,0}(\lambda 1).$$
(23)
In this case the temperature tends to zero as

$$\frac{T}{T_F}\left(\frac{N_b}{N_f}\right)^\lambda .$$
(24)
In Fig. 3 are shown three sample trajectories of Eq. (16). For $`x_0>\lambda 1`$ the final temperature is non-zero when all of the Boltzmann particles have been evaporated. For $`x_0<\lambda 1`$ the final temperature is zero and all of the Boltzmann particles are used in the evaporation. A third curve shows the case $`x_0=\lambda 1`$, which divides the two regimes: this corresponds to constant $`x`$, which implies $`T/T_FN_b/N_f`$.

Note that, for a harmonic potential, $`f(\eta )`$ can be calculated exactly for a Boltzmann gas formules . The critical value of $`\eta `$ which corresponds to $`\lambda 1=0`$ is then $`\eta _c=4.59\mathrm{}`$. For $`\eta <\eta _c`$ the final temperature of the fermions after evaporation of all the Boltzmann particles is always non-zero.

## V Conclusion

We have demonstrated that, in an idealized experiment with no particle losses besides those due to evaporation of the coolant, sympathetic cooling of degenerate Fermi gases is not limited by the heat capacity of the coolant when the latter is a Bose gas with a condensate present. In fact, it is *advantageous* to obtain $`C_f/C_b1`$ in order to minimize the number of bosons that must be evaporated, where $`C_{f,b}`$ are the heat capacities of the fermions and bosons. In this case, $`T=0`$ may be obtained provided that the initial number of bosons is sufficiently large:

$$N_{b,0}>2.7N_fT_0/T_F,$$
(25)
where $`N_f`$ is the number of fermions, $`T_0`$ is the initial temperature which is much less than the Fermi temperature $`T_F`$, and $`\eta =6`$, a typical value of the energy cut-off for evaporative cooling, was chosen.

In contrast, in the case of a Boltzmann coolant, the final temperature of the fermions depends on the ratio of the initial values of the heat capacities $`C_{f,0}/C_{b,0}`$. For typical values of the evaporation cut-off parameter $`\eta `$, $`C_{f,0}/C_{b,0}`$ must be smaller than a critical value on the order of unity to reach $`T=0`$ after evaporation of all the Boltzmann particles.

Acknowledgments: We would like to thank Thomas Bourdel and Jean Dalibard for useful discussions, and Zoran Hadzibabic and Roland Combescot for useful comments on the manuscript. This work was supported by NSF grant no. MPS-DRF 0104447. Laboratoire Kastler Brossel is a research unit of l’Ecole Normale Supérieure and of l’Université Pierre et Marie Curie, associated with CNRS. We acknowledge financial support from Région Ile de France.


text:
# Giant pulses in Millisecond Pulsars

## 1. Introduction

Giant pulses (GPs) have been reported in four pulsars (B0531+21, B1937+21, B1821$``$24 and B0540$``$69) todate (Staelin & Reifenstein 1968; Lundgren et al. 1995; Cognard et al. 1996; Romani & Johnston 2001; Johnston & Romani 2003). Three of these pulsars are millisecond pulsars (MSPs) and also show strongly pulsed hard X-ray profiles (Takahashi et al. 2001). The radio GPs occur in a narrow phase window close to the high energy non-thermal pulse indicating a common magnetospheric origin. All these pulsars have a high value of $`B_{lc}`$. We have used GMRT to search for GPs in candidate MSPs with a non-thermal high energy emission and a range of $`B_{lc}`$ and report detection of such pulses in two more pulsars, PSR J0218+4232 and B1957+20.

We obtained about 3600 s of data on each pulsar using 20 to 22 GMRT antennae in an incoherent mode at 610 MHz with 16 MHz bandwidth. The expected RMS noise in the above configuration of GMRT for a sampling time of 258 $`\mu s`$ is about 1 Jy. The data in two subbands (8 MHz each) were dedispersed to a common sky frequency (610 MHz). The periods with a peak greater than 3.5 times RMS in both bands at the same sample were identified as Large amplitude pulses (LAPs). This procedure disperses any narrow interference spike discriminating against interference. A number of marginal LAPs, i.e. pulses with a peak between 3.0 to 3.5 times RMS, were also identified.

## 2. Results

We searched 2.2 million periods for PSR J0218+4232 ($`B_{lc}`$ 3.2 $`\times 10^5`$ G) and found three significant LAPs. Figure 1a shows the integrated profile for this pulsar and the detected LAPs are marked with filled circles. The largest of these had an intensity 51 times the mean intensity (Intensity of LAP, $`I_{lap}1341`$ Jy-$`\mu `$s). Our data also consisted of 9 marginal LAPs, 7 of which occur between phase 0.89 - 1.2, which is the phase interval corresponding to one of the high energy peaks. We searched about 1 million periods for PSR B1957+20 ($`B_{lc}`$ 3.8 $`\times 10^5`$ G) and found one significant LAP, shown in Figure 1b, with an intensity 129 times the mean intensity ($`I_{lap}925`$ Jy-$`\mu `$s). In addition, we also detected 5 marginal LAPs.

PSRs B1957+20 and J0218+4232 have the fourth and sixth highest values of $`B_{lc}`$ respectively of known radio pulsars. Hence, these new detections support a connection between the magnitude of $`B_{lc}`$ and the existence of GPs and LAPs. The data for pulsars with $`B_{lc}`$ marginally below $`10^5`$ G are being analyzed currently.

## References

Cognard, I. et al. 1996, ApJ, 457, 81

Johnston, S., & Romani, R. W. 2003, private communication

Lundgren, S. C. et al. 1995, ApJ, 453, 433

Romani, R. W., & Johnston, S. 2001, ApJ, 557, L93

Staelin, D. H., & Reifenstein, E. C. 1968, Science, 162, 1481

Takahashi, M. et al. 2001, ApJ, 554, 316


text:
# Gamma-ray lines from cosmic-ray interactions with interstellar dust grains

## 1 Introduction

Observations of nuclear interaction $`\gamma `$-ray lines from the interstellar medium (ISM) would provide a unique tool to study Galactic cosmic-ray ions at non-relativistic energies, as well as the physical conditions of the emitting regions. If the lines produced in the gaseous phase are expected to be significantly Doppler-broadened, some lines produced in interstellar dust grains can be very narrow, because some of the excited nuclei can stop in solid materials before emitting $`\gamma `$-rays (Lingenfelter and Ramaty,, 1977; Ramaty et al.,, 1979). The latter are prime candidates for detection with $`\gamma `$-ray telescopes having high spectral resolution, such as the INTEGRAL spectrometer (SPI).

An illustrative $`\gamma `$-ray spectrum is shown in figure (1). It was obtained in a recent experiment performed at Orsay (Kiener et al., in preparation), in which a thick sample of the Allende meteorite ($``$5 mm in diameter) was bombarded with 10-MeV protons. The Allende meteorite belongs to the class of carbonaceous chondrites and has a composition close to cosmic, except for the volatile elements. The strongest lines (apart from the 511-keV emission) originate from reactions with the most abundant isotopes: 16O, 24Mg, 28Si and 56Fe. One point of particular interest is that the 16O\* line at 6129 keV is much narrower than the two other 16O\* lines at 6916 and 7115 keV. This is because the 6.13-MeV state is relatively long-lived ($`T_{1/2}`$=18.4 ps), such that the recoiling 16O\* nuclei can come to rest in the target before the 6.129-keV $`\gamma `$-ray is emitted, whereas the 6.92- and 7.12-MeV states ($`T_{1/2}`$=4.7 and 8.3 fs, respectively) mostly de-excite in flight. The 6129-keV line, as well as the intense lines at 847 keV from 56Fe\* ($`T_{1/2}`$=6.1 ps), 1369 keV from 24Mg\* ($`T_{1/2}`$=1.35 ps) and 1779 keV from 28Si\* ($`T_{1/2}`$=475 fs) are the prime candidates to scrutinize the interstellar dust grains through their $`\gamma `$-ray emission.

We have calculated in detail the shapes of these four lines as they are produced by cosmic-ray interactions with the ISM. The interaction model and the results are presented in sections (2) and (3), respectively. In section (4), we evaluate the intensities of the predicted emissions from the inner Galaxy.

## 2 Interaction model

To calculate the line profiles, we used a Monte-Carlo method similar to that described in Ramaty et al., (1979). Each event in the simulation corresponds to a $`\gamma `$-ray producing nuclear reaction of a cosmic-ray proton or $`\alpha `$-particle with a target nucleus either in the gas or locked up in a spherical and homogeneous dust grain of the ISM. The $`\gamma `$-rays can be produced by inelastic scattering reactions, e.g. 16O($`p`$,$`p^{}`$$`\gamma _{6129}`$)16O, by $`(p,n)`$ reactions followed by electron capture or positron emission, e.g. 56Fe$`(p,n)`$56Co($`ϵ`$-$`\beta ^+`$)56Fe$`{}_{847}{}^{}{}_{}{}^{}`$ ($`T_{1/2}`$(56Co)=77.2 d), or by spallation reactions, e.g. 20Ne($`p`$,$`p`$$`\alpha `$$`\gamma _{6129}`$)16O. The total cross sections are from Kozlovsky et al., (2002), except for the reactions 24Mg$`(p,n)`$24Al($`\beta ^+`$)24Mg$`{}_{1369}{}^{}{}_{}{}^{}`$ ($`T_{1/2}`$(24Al)=2.053 s) and 28Si$`(p,n)`$28P($`\beta ^+`$)28Si$`{}_{1779}{}^{}{}_{}{}^{}`$ ($`T_{1/2}`$(28P)=270.3 ms). We estimated the contribution of the 24Mg$`(p,n)`$24Al channel from the data of Kiang et al., (1989) and assumed the same relative contribution for the reaction 28Si$`(p,n)`$28P.

The energy distribution of the recoiling excited nuclei were calculated from the differential cross sections. The latter were obtained from a large number of experimental data, completed with extensive coupled-channel calculations with the code ECIS94 (Raynal,, 1994). Two examples of differential inelastic scattering cross sections are shown in figure (2).

The stopping powers of the nuclei recoiling in the grain material were calculated with the code TRIM (Ziegler,, 1985). For these calculations, the grain composition was assumed to be (MgSiFe)O4, which is characteristic of interstellar silicates (e.g. Dwek et al.,, 1997).

To estimate the cosmic-ray proton interstellar spectrum, we used the disk-halo propagation model of Jones et al., (2001). Coulomb and ionization energy losses were taken from Mannheim and Schlickeiser, (1994). For the proton inelastic cross sections, we used the empirical formula given in Moskalenko et al., (2002). Solar-modulated spectra, calculated from the force-field approximation, were fitted to measured proton fluxes (fig. 3). We found the best fit to be provided by the source spectrum $`\dot{Q}(p)dpp^{2.4}dp`$, where $`p`$ is the proton momentum. For simplicity, we used the same form for the $`\alpha `$-particle LIS spectrum as for the protons, with an abundance ratio $`\alpha /p`$=0.1.

## 3 Gamma-ray line profiles

Calculated $`\gamma `$-ray line profiles are shown in figure (4) for different interstellar grain size distributions. We assumed that all the available refractory elements Mg, Si and Fe, and half of the O are locked up in silicate grains, whereas the volatile elements Ne and S (which contribute through spallation reactions to the 6129- and 1779-keV lines, respectively) are in the gaseous phase (see Savage and Sembach,, 1996; Dwek et al.,, 1997). The dotted curves were obtained for the commonly used MRN (Mathis et al.,, 1977) size distribution: an $`a^{3.5}`$ power law in grain radii $`a`$ from $`a_{min}`$=5 nm (see Dwek et al.,, 1997) to $`a_{max}`$=0.25 $`\mu `$m. We see that there are little differences between these spectra and those calculated assuming that all the target nuclei are in the interstellar gas (thick, grey curves). In particular, the two profiles of the 6129-keV line are nearly identical, because almost all of the 16O nuclei excited within dust grains of radius $`<`$0.25 $`\mu `$m escape from the solids before decaying.

The dashed curves were calculated for the same power law distribution, but with $`a_{max}`$=2 $`\mu `$m. The extension of the grain radii beyond the MRN limit of 0.25 $`\mu `$m is motivated by the observations with dust detectors aboard the Ulysses and Galileo spacecrafts of relatively large interstellar grains entering the solar system (Landgraf et al.,, 2000), as well as by models of scattering halos observed around X-ray sources (Witt et al.,, 2001). The thin, solid curves were obtained by further extending the MRN distribution to $`a_{max}`$=10 $`\mu `$m. Such large grains are typical of presolar grains of stardust found in primitive meteorites (Zinner,, 1998) and should exist in various circumstellar media. We see that for these two grain size distributions, the very narrow component of the 24Mg\*, 28Si\* and 16O\* lines could in principle be resolved with SPI. The sensitivity of these line shapes to micrometer-sized particles provides a promising method for tracing large dust grains in the ISM.

## 4 Gamma-ray line fluxes from the inner Galaxy

The $`\gamma `$-ray line intensities from the inner Galaxy were estimated by normalizing the emissivity calculations to the observed flux of high-energy $`\gamma `$-rays ($`>`$70 MeV) due to $`\pi ^0`$ production and decay. We derived the latter from Büsching et al., (2001, fig. 1; EGRET data of observing phases 1-4): $`\mathrm{\Phi }_\gamma ^{\pi ^0}`$$``$8.3$`\times `$10-5 photons cm-2 s-1, for -40$`<`$$`\mathrm{}`$$`<`$40 and -6$`<`$$`b`$$`<`$6. We used the model of Dermer, (1986) to calculate the production of $`\pi ^0`$-decay $`\gamma `$-rays by cosmic-rays with the LIS spectrum of figure (3). We assumed that the metal abundances in the inner Galaxy are on the average twice solar. We then obtained for the fluxes of the 847-, 1369-, 1779- and 6129-keV lines excited in proton and $`\alpha `$-particle interactions with both interstellar gas and dust grains: 5.6$`\times `$10-8, 7.5$`\times `$10-8, 4.3$`\times `$10-8 and 3.1$`\times `$10-7 photons cm-2 s-1, respectively. In comparison, the calculated flux of the relatively strong line at 4438 keV from 12C\* (FWHM$``$150 keV) is 7.1$`\times `$10-7 photons cm-2 s-1 and that of the so-called $`\alpha `$$``$$`\alpha `$ line at $``$450 keV (complex line shape with FWHM$``$80 keV; see Tatischeff et al., 2001) is 3.9$`\times `$10-7 photons cm-2 s-1. All these diffuse emission fluxes are far below the SPI sensitivity, such that a near future detection of nuclear interaction $`\gamma `$-ray lines from the inner Galaxy is unlikely, unless there is a large population of Galactic cosmic-rays with kinetic energies below the threshold for $`\pi ^0`$ production (290 MeV for $`p`$+H collisions). Such a distinct cosmic-ray component predominant at low energies has been proposed to account for the quasi linear correlation between Be and Fe abundances for metallicities \[Fe/H\]$`<`$-1 (e.g. Cassé et al.,, 1995).

We acknowledge F. Boulanger for usefull discussions.


text:
# Comparative studies of the magnetic dipole and electric quadrupole hyperfine constants for the ground and low lying excited states of ²⁵𝑀𝑔⁺

## 1 Introduction

The hyperfine interactions in alkali metal atoms and alkaline earth ions have been of interest for quite a long time . A number of theoretical investigations including some based on relativistic many-body theory have been performed  and they compare reasonably well with experiments. Some of the theoretically predicted values could be of experimental interest with the advent of high precision techniques  involving trapped and laser cooled atoms  and ions .

The high precision calculations of different properties of many-electron atoms requires accurate wave-functions in the nuclear region as well as the region far from the nucleus. The study of properties like hyperfine constants requires the former. Since hyperfine interactions are sensitive to electron correlations, the determination of atomic hyperfine constant provides an important test for *ab initio* atomic structure theory .

In this paper we have carried out *ab-initio* calculations of the magnetic dipole and electric quadrupole hyperfine constants and compared the many-body effects for the ground as well as some excited states. Section 2 provides the necessary theoretical background to the magnetic dipole ($`A`$) and the electric quadrupole ($`B`$) hyperfine constants. Section 3 gives an overview of the coupled-cluster theory and its application to this specific problem. Finally in section 4 the results of our calculations are presented and discussed.

## 2 Theory

The interaction between the various moments of the nucleus and the electrons of an atom are collectively referred to as hyperfine interactions . In this paper we shall consider the interactions between the atomic electrons with the nuclear magnetic dipole and electric quadrupole moments. Nuclear spin gives rise to a nuclear magnetic dipole moment and the departure from spherical charge distribution in the nucleus produces an electric quadrupole moment.

The hyperfine interaction is given by

$$H_{hfs}=\underset{k}{}M^{(k)}T^{(k)},$$
(1)
where $`M^{(k)}`$ and $`T^{(k)}`$ are spherical tensors of rank $`k`$, which corresponds to nuclear and electronic parts of the interaction respectively. The lowest $`k=0`$ order represents the interaction of the electron with the spherical part of the nuclear charge distribution. The eigenstates of the atomic Hamiltonian in the presence of a hyperfine interaction are denoted by $`|IJFM_F`$. Here $`𝐈`$ and $`𝐉`$ are the total angular angular momentum for the nucleus and the electron state, respectively, and $`𝐅=𝐈+𝐉`$ with the projection $`M_F`$.

In the first order perturbation theory, the energy corresponding to the hyperfine interaction of the fine structure state $`|JM_J`$ are the expectation values of $`H_{hfs}`$ such that

$$\begin{array}{ccc}W(J)& =& IJFM_F\left|\underset{k}{}M^{(k)}T^{(k)}\right|IJFM_F\\ & =& \underset{k}{}(1)^{I+J+F}\left\{\begin{array}{ccc}I& J& F\\ J& I& k\end{array}\right\}IM^{(k)}IJT^{(k)}J\end{array}$$
(2)
The $`k=1`$ term describes the magnetic dipole coupling of the nuclear magnetic moment with the magnetic field created by the electron at the position of the nucleus. The nuclear dipole moment $`\mu _I`$ is defined (in units of Bohr magneton $`\mu _N`$) as

$$\mu _I\mu _N=II\left|M_0^{(1)}\right|II=\left(\begin{array}{ccc}I& 1& I\\ I& 0& I\end{array}\right)IM^{(1)}I$$
(3)
and the operator $`T_q^{(1)}`$is given by

$$T_q^{(1)}=\underset{q}{}t_q^{(1)}=\underset{j}{}ie\sqrt{\frac{8\pi }{3}}r_j^2\stackrel{}{\alpha _j}𝐘_{1q}^{(0)}(\widehat{r_j}).$$
(4)
Here $`\stackrel{}{\alpha }`$ is the Dirac matrix and $`𝐘_{kq}^\lambda `$ is the vector spherical harmonics. In Eq.(4) the index $`j`$ refers to the $`j`$-th electron of the atom and $`e`$ is the magnitude of the electronic charge. The magnetic dipole hyperfine constant $`A`$ is defined as

$$A=\mu _N\left(\frac{\mu _I}{I}\right)\frac{JT^{(1)}J}{\sqrt{J(J+1)(2J+1)}},$$
(5)
and the corresponding magnetic dipole hyperfine energy $`W_{M1}`$ is given by

$$W_{M1}=AIJ=A\frac{K}{2},$$
(6)
where $`K=F(F+1)I(I+1)J(J+1)`$.

The second order term in the hyperfine interaction is the electric quadrupole part. The electric quadrupole hyperfine constant is defined by putting $`k=2`$ in Eq. (2). The nuclear quadrupole moment is defined as

$$T_q^{(2)}=\underset{q}{}t_q^{(2)}=\underset{j}{}er_j^3C_q^{(2)}(\widehat{r_j}),$$
(7)
Here, $`C_q^{(k)}=\sqrt{\frac{4\pi }{(2k+1)}}Y_{kq}`$, with $`Y_{kq}`$ being the spherical harmonic. Hence the electric quadrupole hyperfine constant $`B`$ is

$$B=2eQ\left[\frac{2J(2J1)}{(2J+1)(2J+2)(2J+3)}\right]^{1/2}JT^{(2)}J,$$
(8)
and the corresponding electric quadrupole hyperfine energy $`W_{E2}`$ is given by

$$W_{E2}=\frac{B}{2}\frac{3K(K+1)4I(I+1)J(J+1)}{2I(2I1)2J(2J1)}.$$
(9)
In Eq.( 4 and 7 ) $`t_q^{(k)}`$are the single particle reduced matrix element for the electronic part. The reductions of the single particle matrix elements into angular factors and radial integral are straightforward by means of using the Wigner Eckart theorem. These single particle reduced matrix elements are given by

$$\kappa t_q^{(1)}\kappa ^{}=\kappa C_q^{(1)}\kappa ^{}(\kappa +\kappa ^{})𝑑r\frac{\left(P_\kappa Q_\kappa ^{}+Q_\kappa P_\kappa ^{}\right)}{r^2}$$
(10)
and

$$\kappa t_q^{(2)}\kappa ^{}=\kappa C_q^{(2)}\kappa ^{}𝑑r\frac{\left(P_\kappa P_\kappa ^{}+Q_\kappa Q_\kappa ^{}\right)}{r^3},$$
(11)
where $`\kappa C_q^{(k)}\kappa ^{}`$ is the reduced matrix element of the spherical tensor and is equal to

$$(1)^{j+1/2}\sqrt{(2j+1)(2j^{}+1)}\left(\begin{array}{ccc}j& k& j^{}\\ \frac{1}{2}& 0& \frac{1}{2}\end{array}\right)\pi (l,k,l^{}),$$

with

$$\pi (l,k,l^{})=\{\begin{array}{c}\begin{array}{cc}1& \mathrm{if}l+k+l^{}\mathrm{even}\\ 0& \mathrm{otherwise}\end{array}\end{array}$$

Here the single particle orbitals are expressed in terms of the Dirac spinors with $`P_i`$ and $`Q_i`$ as large and small components respectively.

## 3 Overview of the coupled cluster theory : method of calculation

We start with an $`N`$ electron closed shell Dirac-Fock (DF) reference state $`|\mathrm{\Phi }`$. The corresponding correlated closed shell state is then

$$|\mathrm{\Psi }=\mathrm{exp}(T)|\mathrm{\Phi },$$
(12)
where $`T`$ is the core electron excitation operator. Then the Dirac-Coulomb eigenvalue equation is

$$H\mathrm{exp}(T)|\mathrm{\Phi }=E\mathrm{exp}(T)|\mathrm{\Phi },$$
(13)
with the Dirac-Coulomb Hamiltonian

$$H=\underset{i}{}\left(c\alpha _ip_i+(\beta _i1)mc^2+V_N\right)+\underset{i<j}{}\frac{1}{r_{ij}}.$$
(14)
This leads to the exact ground state energy $`E`$ of the closed-shell part of the system. Here $`\alpha _i`$ and $`\beta _i`$ are Dirac matrices and $`V_N`$ is the nuclear potential. If we consider the DF state $`|\mathrm{\Phi }`$ as the Fermi vacuum, then the normal ordered Hamiltonian is

$$H_NH\mathrm{\Phi }\left|H\right|\mathrm{\Phi }=HE_{DF}.$$
(15)
If we project $`\mathrm{\Phi }|\mathrm{exp}(T)`$ from the left we obtain the correlation energy ($`\mathrm{\Delta }E`$) and if we project any of the excited determinant $`\mathrm{\Phi }^{}|\mathrm{exp}(T)`$ we additionally get a set of equations which are used to obtain the $`T`$ amplitudes. Using the normal ordered dressed Hamiltonian $`\overline{H_N}=\mathrm{exp}(T)H_N\mathrm{exp}(T)`$ the corresponding equations for correlation energy and amplitudes become

$$\mathrm{\Phi }\left|\overline{H_N}\right|\mathrm{\Phi }=\mathrm{\Delta }E,$$
(16)
and

$$\mathrm{\Phi }^{}\left|\overline{H_N}\right|\mathrm{\Phi }=0.$$
(17)
Here the state $`|\mathrm{\Phi }^{}`$ may be singly excited $`|\mathrm{\Phi }_a^r`$ or double excited $`|\mathrm{\Phi }_{ab}^{rs}`$ and so on. The indices $`a,b,\mathrm{}`$ refer to holes and $`p,q,\mathrm{}`$ to particles. We have considered the coupled cluster single and double (CCSD) approximation, where the cluster operator $`T`$ is composed of one- and two-body excitation operators, *i.e.* $`T=T_1+T_2`$, and are expressed in second quantization form

$$T=T_1+T_2=\underset{ap}{}a_p^{}a_at_a^p+\frac{1}{2}\underset{abpq}{}a_p^{}a_q^{}a_ba_at_{ab}^{pq}.$$
(18)
Contracting the ladder operators  and rearranging the indices, the amplitude equations can be expressed in the form

$$A+B(T)T=0,$$
(19)
where $`A`$ is a constant vector consisting of the matrix elements $`\mathrm{\Phi }^{}\left|H_N\right|\mathrm{\Phi }`$, $`T`$ is the vector of the excitation amplitudes and $`B(T)`$ is the matrix which depends on the cluster amplitudes itself so that Eq. (19) is solved self-consistently. For example, a typical contribution to the term $`\widehat{\widehat{H_NT_2}T_2}`$ is

$$B_{ab}^{pq}=\frac{1}{2}\underset{dgrs}{}V_{dgrs}t_{ad}^{pr}t_{gb}^{sq}.$$
(20)
Here $`V_{dgrs}`$ is the two-electron Coulomb integral and $`t_{ad}^{pr}`$ is the cluster amplitude corresponding to a simultaneous excitation of two electrons from orbital $`a`$ and $`d`$ to $`p`$ and $`r`$ respectively. To obtain a full set of terms which contribute to this specific excitation, diagrammatic techniques are used.

The ground state of $`{}_{}{}^{25}Mg^+`$ contains only one valance electron in the outer most orbital ($`3s_{1/2}`$). To calculate the ground state energy of the system we first compute the correlations for the closed shell system ($`{}_{}{}^{25}Mg^{+2}`$) using the closed shell coupled cluster approach and then use the technique of electron attachment (open shell coupled cluster (OSCC)) method. The energy of the excited state are obtained by the same way. In order to add an electron to the $`k`$th virtual orbital of the DF reference state we define

$$|\mathrm{\Phi }_k^{N+1}a_k^{}|\mathrm{\Phi }$$
(21)
with the particle creation operator $`a_k^{}`$. Then by using the excitation operators for both the core and valance electron the exact state is defined as :

$$|\mathrm{\Psi }_k^{N+1}=\mathrm{exp}(T)\left\{\mathrm{exp}(S_k)\right\}|\mathrm{\Phi }_k^{N+1}.$$
(22)
Here $`\left\{\mathrm{exp}(S_k)\right\}`$ is the normal ordered exponential representing the valance part of the wave operator. Here

$$S_k=S_{1k}+S_{2k}=\underset{kp}{}a_p^{}a_ks_k^p+\frac{1}{2}\underset{bpq}{}a_p^{}a_q^{}a_ba_ks_{kb}^{pq},$$
(23)
where $`k`$ stands for valance orbital. $`S_k`$ contain the particle annihilation operator $`a_k`$, and because of the normal ordering it cannot be connected to any other valance electron excitation operator and then $`\left\{\mathrm{exp}(S_k)\right\}`$ automatically reduces to $`\left\{1+S_k\right\}`$.

Then we can write the Eq.(22) as

$$|\mathrm{\Psi }_k^{N+1}=\mathrm{exp}(T)\left\{1+S_k\right\}|\mathrm{\Phi }_k^{N+1}.$$
(24)
Following the same procedure as in the closed-shell approach, we obtain a set of equations

$$\mathrm{\Phi }_k^{N+1}\left|\overline{H_N}\left\{1+S_k\right\}\right|\mathrm{\Phi }_k^{N+1}=H_{eff}$$
(25)
and

$$\mathrm{\Phi }_k^{{}_{}{}^{}N+1}\left|\overline{H_N}\left\{1+S_k\right\}\right|\mathrm{\Phi }_k^{N+1}=H_{eff}\mathrm{\Phi }_k^{{}_{}{}^{}N+1}\left|\left\{1+S_k\right\}\right|\mathrm{\Phi }_k^{N+1},$$
(26)
where the desired roots can be obtained by diagonalizing $`H_{eff}`$. The Eq.(26) is non-linear in $`S_k`$ because the energy difference $`H_{eff}`$ is itself a function of $`S_k`$. Hence, these equations have to solved self-consistently to determine the $`S_k`$ amplitudes.

Triple excitations are included in our open shell CC amplitude calculations by an approximation that is similar in spirit to CCSD(T) . The approximate triple excitation amplitude is given by

$$S_{abk}^{pqr}=\frac{\widehat{VT_2}+\widehat{VS_2}}{\epsilon _a+\epsilon _b+\epsilon _k\epsilon _p\epsilon _q\epsilon _r},$$
(27)
where $`S_{abk}^{pqr}`$ are the amplitudes corresponding to the simultaneous excitation of orbitals $`a,b,k`$ to $`p,q,r`$ respectively and $`\widehat{VT}`$ and $`\widehat{VS}`$ are the correlated composites involving $`V`$ and $`T`$, and $`V`$ and $`S`$ respectively. $`\epsilon _k`$ is the orbital energy of the $`k`$th orbital. The above amplitudes (some representative diagrams are given in figure 1) are added appropriately in the singles and doubles open shell cluster amplitude equations and these equations are then solved self-consistently. We therefore obtain solutions of $`S_1`$and $`S_2`$ amplitudes taking into consideration the effect of the triple excitations in an approximate way

The expectation value of any operator $`O`$ can be written as the normalized form with respect to the exact state $`|\mathrm{\Psi }^{N+1}`$ as

$$O=\frac{\mathrm{\Psi }^{N+1}\left|O\right|\mathrm{\Psi }^{N+1}}{\mathrm{\Psi }^{N+1}|\mathrm{\Psi }^{N+1}}=\frac{\mathrm{\Phi }^{N+1}\left|\left\{1+S^{}\right\}\mathrm{exp}(T^{})O\mathrm{exp}(T)\left\{1+S\right\}\right|\mathrm{\Phi }^{N+1}}{\mathrm{\Phi }^{N+1}\left|\left\{1+S^{}\right\}\mathrm{exp}(T^{})\mathrm{exp}(T)\left\{1+S\right\}\right|\mathrm{\Phi }^{N+1}}.$$
(28)
For computational simplicity we store only the one-body matrix element of $`\overline{O}=\mathrm{exp}(T^{})O\mathrm{exp}(T)`$. $`\overline{O}`$ may be expressed in terms of uncontracted single-particle lines . The fully contracted part of $`\overline{O}`$ will not contribute as it cannot be linked with the remaining part of the numerator of the above equation.

In the LCCSD approximation Eq. (22) turns out to be

$$|\mathrm{\Psi }_k^{N+1}=\left\{1+T+S_k\right\}|\mathrm{\Phi }_k^{N+1},$$
(29)
and $`\overline{H_N}=H_N+\stackrel{}{H_NT}`$.

The closed and open shell cluster amplitude equations reduce to

$$\mathrm{\Phi }_0^k\left|H_N+\stackrel{}{H_NT}\right|\mathrm{\Phi }_0=0,$$
(30)
and

$$\mathrm{\Phi }_k^{{}_{}{}^{}N+1}\left|\overline{H_N}\left\{1+S_k\right\}\right|\mathrm{\Phi }_k^{N+1}=\mathrm{\Phi }_k^{{}_{}{}^{}N+1}\left|S_v\right|\mathrm{\Phi }_k^{N+1}\mathrm{\Phi }_k^{N+1}\left|H_N\right|\mathrm{\Phi }_k^{N+1}.$$
(31)
The orbitals used in the present work are expanded in terms of a finite basis set comprising of Gaussian type orbitals (GTO)

$$F_{i,k}(r)=r^k\mathrm{exp}(\alpha _ir^2),$$
(32)
with $`k=0,1,2\mathrm{}`$ for $`s,p,d,\mathrm{}`$ type functions, respectively. The exponents are determined by the even tempering condition

$$\alpha _i=\alpha _0\beta ^{i1}.$$
(33)
The staring point of the computation is the generation of the Dirac-Fock (DF) orbitals  which are defined on a radial grid of the form

$$r_i=r_0\left[\mathrm{exp}(i1)h1\right]$$
(34)
with the freedom of choosing the parameters $`r_0`$ and $`h`$. All DF orbitals are generated using a two parameter Fermi nuclear distribution

$$\rho =\frac{\rho _0}{1+\mathrm{exp}((rc)/a)},$$
(35)
where the parameter $`c`$ is the half charge radius and $`a`$ is related to skin thickness, defined as the interval of the nuclear thickness in which the nuclear charge density falls from near one to near zero.

Although we have used a large basis for the generation of the single particle orbitals, the high-lying virtual orbitals (above a certain threshold) are kept frozen as their contributions to the high-lying virtuals in the $`T`$ and $`S`$ amplitudes in the CC equations are negligible. Another advantage of this approximation is that it reduces the memory required to store the matrix elements of the dressed operator $`\overline{H}`$ and the two-electron Coulomb integrals in the main memory, thereby reducing the computational cost. In our calculations, we have included all possible single, double and partial triple excitations from the core.

We have used different $`\alpha _0`$ and $`\beta `$ for different symmetries. The number of basis functions used to generate the even tempered DF states are listed in table 1 and the values of the parameters $`\alpha _0`$ and $`\beta `$ used are also listed. For coupled cluster calculations, we have restricted the basis by imposing an upper bound in energy for single particle orbitals and the convergence of our results are shown in table 3.

## 4 Results and discussions

The values of the magnetic dipole hyperfine constant $`A`$ and electric quadrupole hyperfine constant $`B`$ for different states are given in table 2. Our calculated value of $`A`$ for the ground state $`3s_{1/2}`$ is in good agreement (less than 0.6%) with experiment and it is more accurate than a previous calculation ($``$1%) based on second order relativistic many-body perturbation theory (RMBPT) . This is because unlike the previous work our calculation is based on an approach which is equivalent to all-order MBPT. In particular, we have taken into account all single, double and a subset of triple excitations to all orders in the residual Coulomb interaction. It is, therefore, not surprising that the result of our calculation when carried out by using second order RMBPT is in agreement with the result of the earlier calculation (see table 2) .

In table 2 the results calculated by Safronova *et al*  using linearised coupled cluster in singles doubles approximation (LCCSD) are given. We have performed LCCSD(T) calculations and the results are given in the same table. The reasons for the discrepancy between the two linearised coupled-cluster calculations seem to be the inclusion of partial triple excitations by us and the different choice of basis sets. From our present calculations it is possible to determine the contributions from the non-linear clusters. These contributions vary from 0.3% ($`3s_{1/2}`$) to 0.8% ($`3d_{3/2}`$) for the different states. We find that if we take the linear contributions of the $`T`$ amplitudes but the nonlinear contributions of the combined $`T`$ and $`S`$ amplitudes and then perform the calculation of the hyperfine constants (Eq.28) at the linear level (which we have named as (L)CCSD(T)), the result ($`A(3s_{1/2})=597.45`$ MHz) is in excellent agreement with the experiment (error is $``$0.2%) for the ground state. But theoretically this approach is not complete and hence a proper inclusion of the nonlinear effects is desirable. Table 4 shows the contribution of triples in the final property. Some partial triples effect are taken into account iteratively according to Eq. (27) and the $`S`$ amplitude thus contain some partial triples. The effect of triples in calculation of properties (Eq. 28) thus comes from the triples effect included in $`S`$ amplitude. Figure 1 shows some of the typical diagrams which give rise to the effect due to triples, where figure 1a and 1b correspond to the terms $`\widehat{VS_2}`$ and $`\widehat{VT_2}`$ respectively.

For the present calculations, the number of basis functions actually used are the following : 16$`s_{1/2}`$, 14$`p`$ ($`p_{1/2}`$and $`p_{3/2}`$), 11$`d`$ ($`d_{3/2}`$and $`d_{5/2}`$) and 10$`f`$ ($`f_{5/2}`$and $`f_{7/2}`$). Excitations from all the core orbitals have been considered.

The important contributions to the magnetic dipole hyperfine constants for different states are given in table 5. In particular, we have analyzed the contributions from various many body effects and have demonstrated that the most important contributions come from core polarization and pair correlation effects. The largest contribution comes from $`\overline{O}`$. The next two largest contributions come from $`(\overline{O}S_1+S_1^{}\overline{O})`$ and$`(\overline{O}S_2+S_2^{}\overline{O})`$ which correspond to the pair-correlation (PC) and core-polarization (CP) effects respectively. The contribution from the corresponding MBPT terms are listed in table 6. Figures 2 and 3 represents the pair-correlation and core-polarization diagrams in MBPT and CC respectively.

We have listed the contributions from the different terms containing the dressed operator $`\overline{O}`$ in table 5 and table 6 gives the contributions from the terms containing the operator $`O`$ directly. The results given in tables 5 and 6 show that the CP contribution is larger than the PC in magnitude for all the states, although the ratio of the two effects is not uniform. It is important to note that the former contribution includes the hyperfine interaction of all the core orbitals while only a specific valence orbital is involved in this interaction for the latter (see Figure 3). These individual contributions are presented in table 7 for the ground state of $`Mg^+`$. However, the ground state hyperfine constant $`A`$ for $`Ba^+`$exhibits exactly the opposite behaviour . Even though $`Ba^+`$ has more core electrons than $`Mg^+`$, the relativistic enhancement of the valence ($`6s`$) magnetic dipole hyperfine interaction results in the value of PC exceeding that of CP.

## 5 Conclusion

In this paper we have carried out *ab-initio* relativistic coupled cluster calculations of magnetic dipole ($`A`$) and electric quadrupole ($`B`$) hyperfine constants for the ground and some excited states of $`{}_{}{}^{25}Mg^+`$. Some contributions from partial triples are also taken into account in our calculation. We have shown that the dominant many-body contributions to these properties come from core-polarization and pair-correlation effects.

In addition to comparing with the available experimental data we have also predicted the values of $`A`$ and $`B`$ for a few states which could be of interest in the future. Using ion-trapping and other experimental techniques, it may be possible to measure both the magnetic dipole and electric quadrupole hyperfine constants for different states of $`{}_{}{}^{25}Mg^+`$, thereby checking the accuracy of our calculations. This would constitute an useful test of the validity of the coupled-cluster theory in capturing the many-body effects in hyperfine interactions in light atomic systems with a single valence electron.

> Acknowledgments : This work was supported by the BRNS for project no. 2002/37/12/BRNS. The computation were carried out on our group’s 4 CPU E450 Sun Ultra SPARC machine in IIA and on CDAC, Bangalore’s Teraflop Supercomputer Param Padma.


text:
# Effects of Stellar-Mass and the ISM on the Evolving Morphologies of Planetary Nebulae

## 1. Introduction

Towards the end of the Asymptotic Giant Branch (AGB) phase, stellar evolution predicts episodic mass-loss increases, a consequence of the thermal pulses in the star. Although the occurrence of this modulated mass-loss has played a key role in the interpretation of the different shells found around PNs and AGB stars, hydrodynamic models describing the evolution of the wind in this phase still need to be considered. Because of the inherent difficulties in the mass-loss calculations and in recovering the history of mass-loss from observational studies, the exact evolution of mass-loss during the AGB still remains unknown. For this it is fundamental that the grid used in the numerical computations allow the study of the whole stellar ejecta. It is only then, by comparing the models with the PN structure at large scales that constrains can be placed on the treatment of mass-loss in models stellar evolution models. Previous models have always truncated the computational grids to small scales.

Stellar evolution predicts that stars with main sequence masses in the range of $``$1–5 M will produce PN, whilst PN nuclei and white dwarfs mass distributions peak around 0.6 M. Since most of the mass-loss occurs on the AGB phase, it should be easily observable as ionized mass during the PN stage. However, observations of Galactic PNs reveal on average only 0.2 M of ionized gas.

Another important aspect that needs to be considered is the role played by the stellar progenitor mass in the the PN formation. The central star provides the wind and the radiation field that determines the evolution of the nebular gas during this phase. The details of the post-AGB evolution of the star, and therefore the energy injected in the nebular gas are mainly determined by its core mass. Despite this, most of the numerical studies of PN evolution in the literature have been restricted to a 0.6 M post-AGB evolutionary track.

## 2. The evolution of the stellar ejecta during the AGB phase.

In Villaver, García-Segura, & Manchado (2002) we studied the time-dependent hydrodynamics of the circumstellar gas shells of AGB stars. We directly used the results of stellar evolution (predictions for the wind evolution from Vassiliadis & Wood 1993) as inputs to the simulations.

We find that the wind variations associated with the thermal pulses lead to the formation of transient shells with an average lifetime of 20,000 yr and, consequently, do not remain recorded in the density or velocity structure of the gas. The formation of shells that survive at the end of the AGB phase occurs via two main processes: shocks between the shells formed by two consecutive enhancements of the mass-loss or continuous accumulation of the material ejected by the star in the interaction region with the ISM. We do not find the signature of discrete mass-loss recorded in the density or in the velocity structure of the circumstellar envelope (CSE). Consequently, we argued against the use of the different observed shells as dynamical clocks.

The mass of the AGB stellar progenitors are usually estimated from observing the CSEs, with the assumption that all the observed mass has been ejected by the star. We find, however, that the final mass of the CSE contains a significant fraction of swept-up ISM material. Thus, in order to obtain unbiased AGB-progenitor mass estimates, it is vital that this effect is taken into account. We also predict that the CSEs are a mix of the ISM and the wind material that has been enriched by the stellar interior and brought to the surface of the star by dredge-up processes. This effect should be taken into account for abundance analysis. According to our simulations based in the mass-loss predictions of one particular set of stellar evolution models, large CSEs (up to 2.5 pc) are expected around stars at the tip of the AGB.

## 3. The connection between the stellar progenitor and the PN shell’s evolution

To study PN formation we used a set of computational grids large enough to study the full stellar ejecta, and a set of grids small enough to resolve the processes taking place close to the central star. In Villaver, Manchado, & García-Segura (2002) we considered the evolution of the post-AGB wind and the ionizing radiation field for PN nuclei evolving from progenitors with initial masses between 1 and 5 M (taken from Vassiliadis & Wood 1994). We show the importance of the dynamical effects of ionization on the shell evolution, which can account for the observed disagreement between the kinematical ages and the age of the CSs for different progenitor masses. We find that the evolution of the main shell is controlled by the ionization front rather than by the thermal pressure provided by the hot bubble during the early PN stages.

The halos have sizes up to 2.3 pc (more than twice the size of the main shell) and are formed during the AGB phase. They produce H$`\alpha `$ emissivities between 10 and 5000 times fainter than the main shell and contain most of the total ionized mass lost by the stellar progenitors. Given the low surface brightness of the halos, previous observations have underestimated the ionized masses in PNs.

## 4. The interaction of PNs with the ISM

Several PNs show bow-shock structures, suggesting the interaction of the nebular shell with the ISM while the star is moving. We have approached the problem of PN-ISM interaction through a realistic perspective by considering a low-mass star evolving during the AGB and PN phases while moving through the ISM. In Villaver, García-segura, & Manchado (2003) we showed that even the ejecta of a star with a systemic velocity of 20 km s-1 moving through a low density medium will interact with it and form bow-shock structures qualitatively similar to those observed.

An increase of ram pressure (due to e.g a higher velocity) leads to the development of instabilities that lead to a partial fragmentation of the shell and to a more efficient mixing with the ISM material. This fragmentation allows the UV radiation field to escape from the nebula at certain locations. In Figure 1 we show the logarithm of the gas density during the latest stages of the evolution for a star moving with 80 km s-1 through a low density ISM (no=0.05 cm -3). Figure 1 shows only half of the $`r\theta `$ plane, (where $`r`$ and $`\theta `$ are the radial and polar coordinates respectively) the star is fixed in the grid and the ISM flows in from the top to the bottom. We assume that the ISM moves relative to the star perpendicular to the line-of-sight by fixing the position of the star at the center of the grid and allowing the ISM to flow into the grid at the outer boundary from 0 to 90. In Szentgyorgyi et al. (2003) we considered a stellar motion with a velocity of 85 km s-1 in the study of the PN NGC 246. We qualitatively reproduce the overall shape of the PN, due to the interaction which we find consistent with what is expected in the fast interaction with a rarefied medium at the position of the PN in the Galaxy.

We find that due to ram-pressure stripping, most of the mass ejected during the AGB phase is left downstream of the star in its motion, an effect that might be able to account for the small amount of ionized mass recovered in PN shells. We conclude that the interaction with the ISM plays a major role in the PN formation process even during the early AGB evolution and it’s interaction with the ISM cannot be studied using simple ram pressure balance arguments.

## 5. Summary

We have studied the PN formation by following the evolution of the stellar wind as predicted by stellar evolutionary models and considering the influence of the external ISM. We find that although the mass-loss history during the AGB phase is very different for low- and high-mass progenitors, the final nebular structure is very similar. The mass-loss during the AGB gives rise to the formation of large shells (with sizes up to 3 pc) that contain most of the mass lost by the star plus an additional amount (up to 1 M of ISM material) which is ISM material swept up by the stellar wind. We find that the movement of the CS with respect to its surrounding medium considerably alters the PN formation. The main effects of the interaction, apart from that on the morphology, are that the total size of the outer PN shell (halo) is reduced considerably and most of the mass ejected during the AGB phase is stripped by the ram pressure of the ISM and left in the downstream direction of the stellar movement. The mass stripped away by the ISM when the star is moving might be able by itself to account for the problem of the missing ionized mass in PNs.

## References

Szentgyorgyi, A., Raymond, J., Franco, J., Villaver, E., & López-Martín, L. 2003, ApJ, 594, 874

Vassiliadis, E. & Wood, P. R. 1994, ApJS, 92, 125

Vassiliadis, E. & Wood, P. R. 1993, ApJ, 413, 641

Villaver, E., García-Segura, G., & Manchado, A. 2003, ApJ, 585, L49

Villaver, E., Manchado, A., & García-Segura, G. 2002, ApJ, 581, 1204

Villaver, E., García-Segura, G., & Manchado, A. 2002, ApJ, 571, 880