text:
# Hybridization-induced superconductivity from the electron repulsion on a tetramer lattice having a disconnected Fermi surface

## Abstract

Plaquette lattices with each unit cell containing multiple atoms are good candidates for disconnected Fermi surfaces, which are shown by Kuroki and Arita to be favorable for spin-flucutation mediated superconductivity from electron repulsion. Here we find an interesting example in a tetramer lattice where the structure within each unit cell dominates the nodal structure of the gap function. We trace its reason to the way in which a Cooper pair is formed across the hybridized molecular orbitals, where we still end up with a $`T_c`$ much higher than usual.

The discovery of high-$`T_c`$ superconductors has kicked off renewed interests in electronic mechanisms for superconductivity. In the course of studies it is becoming increasingly clear that superconductivity can arise from repulsive interaction between electrons. The essence there is electrons interact by exchanging spin fluctuations, which can effectively act as an attraction in the gap equation, if we consider anisotropic pairing with nodes in the gap function.

For the Hubbard model, this was first suggested from early calculations by Scalapino et al. . Subsequently a quantum Monte Carlo calculation  indeed indicated an enhancement of the pairing correlation with $`d`$-wave symmetry in the repulsive Hubbard model. The superconducting critical temperature ($`T_c`$) has been estimated with the fluctuation-exchange approximation (FLEX), a kind of renormalized random phase approximation.

One remarkable point is $`T_cO(0.01t)`$, esitmated for the repulsive Hubbard model in the two-dimensional (2D) square lattice, is two orders of magnitudes smaller than the starting electronic energy (i.e., the hopping integral $`t`$), although this gives the right order for the curates’ $`T_c`$. Arita et al. then looked at various lattices (square, trianglar, fcc, bcc, etc) in search of the most favorable case with a FLEX analysis. The best case, as far as these ordinary lattices are concerned, turns out to be the 2D square lattice, so $`T_c<O(0.01t)`$ remains. As discussed in Ref. , there are good reasons why $`T_c`$ is so low, where an important one is the presence of nodes in the superconducting gap function greatly reduces $`T_c`$: While the main pair-scattering, across which the gap function has opposite signs to make the effective interaction attractive, some of the pair scatterings around the nodes have negative contributions to the effective attraction by connecting $`k`$-points on which the gap has the same sign.

So a next important avenue to explore is: can we improve the situation by going over to multiband systems. Kuroki and Arita  have shown that this is indeed the case if we have disconnected Fermi surfaces. In this case $`T_c`$ is dramatically enhanced, because the sign change in the gap function can avoid the Fermi pockets, where all the pair-scattering processes contribute positively.  This has been numerically shown to be the case for the triangular lattice (for spin-triplet pairing)  and a squre lattice with a period-doubling , where $`T_c`$ as estimated with FLEX is as high as $`O(0.1t)`$.

To be more precise, the key ingredinents are: (a) when the Fermi surface is nested, the spin susceptibility $`\chi (𝒒,\omega )`$ has a peak, with a width $`\mathrm{\Delta }𝑸`$. (b) When a multiband system with a disconnected Fermi surface has an inter-pocket nesting (i.e., strong inter-pocket pair scattering and weak intra-pocket one) the gap function has the same sign ($`s`$-wave symmetry) within each pocket, and the nodal lines can happily run in between the pockets.

Here a natural question is: to prepare multiband systems systematically, we can consider lattices comprising some units, or “plaquette lattices”. This is indeed considered by Kimura et al.  who conceived an idea of actually fabricating the structure from arrays of quatum dots. One example is the inset of Fig. 1, where the unit is a square lattice of squres, where the intra-plaquette transfer is stronger than the inter-plaquette one. Now, an important question is: how the structure of the superconducting gap function having higher $`T_c`$’s is determined and whether the $`s`$-wave symmetry on each pocket is compulsory.

In this paper, we study the repulsive Hubbard model on a square lattice of diamonds, where a diamond is a four-site unit rotated by 45 degrees (Fig. 1). One motivation is that this structure is reminiscent of the in-plane lattice of B and C in CaB2C2 , so should not be too unrealistic, although we are not claiming to consider this particular material. We have solved Eliashberg’s equation  for the Hubbard model with the FLEX and obtained a $`T_c`$ that is considerably higher than that of the square lattice. More importantly, however, the present model, despite of its disconnected Fermi surface, has sign changes in the gap function within each pocket, so a different mechanism should be at work. We can trace back its reason in real space that the singlet pairing here results from a hybridization of two molecular orbitals in each plaquette.

So a message of this work is that, in order to search for higher $`T_c`$’s, we should look not only at the Fermi surface but also at the atomic structure in real space.

We have used the four-band version of FLEX, where the Green’s function $`G`$, the spin susceptibility $`\chi `$, the self-energy $`\mathrm{\Sigma }`$, and the superconducting gap function $`\varphi `$ are $`4\times 4`$ matrices, e.g., $`G_{lm}(𝒌,i\epsilon _n)`$, where $`l,m`$ specify the four sites in a unit cell. We can go from the site indeces over to band indeces with a unitary transformation. For the spin susceptibility, we concentrate on its largest eigenvalue denoted as $`\chi `$. From Green’s function and the spin susceptibility we obtain the superconducting gap function (with spin-singlet pairing assumed) and $`T_c`$ by solving the linearized Eliashberg’s equation. In our analysis, we take $`32\times 32`$ $`k`$-point meshes and up to 4096 Matsubara frequencies, where the numerical results are sufficiently converged. The on-site Hubbard repulsion is taken to be a typical strong-correlation value, $`U/W=7/6`$, where $`W`$ is the band width. The band filling (=number of electrons/number of sites) is taken to be a value $`n=0.85`$, which is close, but not too close, to the half-filling where the Mott transition and antiferromagnetic order are expected. In the below we have confirmed that when the eigevalue of Eliashberg’s equation becomes unity prior to the divergence of the spin susceptibility.

The FLEX result for the spin susceptibility, Green’s function, and the superconducting gap function are shown in Fig. 2 for $`t^{}=0.6`$ ($`t=1.2`$) with $`T=0.06`$. Ridges in the Green’s function delineate the Fermi surface. Fermi surface comprises two bands (second and third from the bottom), which we call band A and band B. We can see that each band has a squre-like pocket for the Fermi surface, which is reminiscent of the Fermi surface of Sr2RuO4, with a good reason as we shall see. On the other hand we notice that the spin susceptibility has broad peaks around ($`\pi ,\pi `$). Quite unexpectedly, the gap function has a strange structure in such as way that its amplitude on the Fermi surface is peaked at the corners rather than along the edges.

Figure 3 shows $`T_c`$ as a function of $`t`$. A plaquette lattice is characterized by the intra-plaquette transfer $`t`$ and inter-plaquette one $`t^{}`$. To facilitate comparison, we have fixed the single-particle bandwidth $`W=2(2|t|+|t^{}|)`$ to $`6`$, which is the value of $`W`$ when $`t=t^{}=1`$. Interestingly, $`T_c`$ is larger for smaller $`|t^{}|`$, which is in accord with our expectation that the case of plaquette lattice (tetramirization in the present case) is favorable. $`T_c`$ saturates to a value $`T_c=0.06=0.01W`$, which is two times greater than that ($`<0.03t`$) for the square lattice.

The enhancement of $`T_c`$ in smaller $`|t^{}|`$ can be understood by the nesting between the disconnected Fermi surface (Fig. 4). At half filling ($`n=1`$) the Fermi surface is perfectly nested regardless of the value of $`|t^{}|`$. When less than half-filled, the nesting in the band depicted by $`𝑸`$ in Fig. 3 degrades but not so much for smaller $`|t^{}|`$. Because the effective attraction here is mediated by the antiferromagnetic spin fluctuation (with a large $`|Q|`$), $`T_c`$ is larger for better nesting, i.e., for smaller $`|t^{}|`$ the warping of the quasi-1D Fermi surface. This $`t^{}`$-dependence of $`T_c`$ can be understood also in real space, where the intra-cell spin singlet should be more robustly formed for smaller $`|t^{}|`$.

We can explain the structures of the strange gap function in terms of molecular orbitals (Fig. 5a) and the pairing (Fig. 5b) in real space. As mentioned, the Fermi surface reminds us of the ($`\alpha `$ and $`\beta `$) bands in SrRuO, which consists of two sets of quasi-1D array of orbitals. In the present case the quasi-1D band structure arises from the molecular orbitals in the plaquette. The molecular orbitals that are relevant to bands A and B have $`p_x`$\- and $`p_y`$\- symmetries (Fig. 5a), so that we do have two sets of 1D arrays of orbitals. Due to the hybridization of the two orbtials, the two sets of quasi-1D Fermi surface (two, orthogonal sets of parallel lines) anticross, and we end up with two square-shaped pockets.

Now the interesting point is, the Cooper pair, arising from antiferromagnetic fluctuations, is formed across the adjacent sites (ellipses in Fig. 5b), that is across the two molecular orbitals. This implies that the gap function on the Fermi surface is large only around the positions at which the two ($`p_x`$ and $`p_y`$) bands hybridize, i.e., around the anticrossing points. The pairings depicted in Fig. 5b as solid ellipses and dashed ellipses have opposite signs, so that the symmetry is $`d_{xy}`$. We have thus a curious case of an inter-band pairing (a pairing formed by inter-band pair scatterings). We can in fact confirm that the pairing is interband and $`d_{xy}`$ by looking at the sign of the gap function. The ellipses is Fig. 5 are represented by $`C_{A(B)}C_{A(B)}`$, where $`C_{A(B)\sigma }C_{1\sigma }C_{2(4)\sigma }+C_{4(2)\sigma }C_{3\sigma }`$, where $`C_{i\sigma }`$ annihilates an electron with spin $`\sigma `$ at site $`i`$ as numbered in Fig.5. When the pairing is formed across adjacent sites we can show $`C_AC_A=C_BC_B`$ (opposite signs across the bands) with an odd parity ($`d_{xy}`$) against $`k_xk_x`$ or $`k_yk_y`$.

This explains why the gap function changes sign within each pocket. Hence the $`s`$-wave symmetry on each pocket is not a necessary condition for a high $`T_c`$. Even when there are intra-pocket sign changes, the disconnectivity of the Fermi surface plays an important role, since the inter-pocket sign change associated with $`𝑸`$ does exploit this in avoiding intersecting the Fermi surface.

One might then wonder why the intra-pocket sign changes do not reduce $`T_c`$, since the pair-scattering processes depicted by $`𝑸^{}`$ in Fig. 4 connecting the $`k`$-points across which the gap function has the same sign would normally reduce $`T_c`$. In the present case, however, we have an ingeneous situation where the gap function is strongly peaked at the position (corners in this example) where hybridization occurs, so that the $`𝑸^{}`$ processes, being off the peak, have little effect.

How is the present model compared with the lattice considered in Ref. , where the plaquette is a four-membered ring in either case. Here the consideration in real space comes in handy. The latter case is a 2D system rather than quasi-1D, where the pockets arise due to a band folding due to the tetramerization,. The pairing in real space is formed within each molecular orbital (Fig. 6) unlike the present case, so the superconducting gap function has an $`s`$-wave symmetry on each Fermi surface pocket (but opposite signs across the two pockets with the nesting $`𝑸0`$ due to the band folding (inset of Fig. 4). Despite of this distinction, the resulting $`T_c`$’s are similar between the two cases.

Another interesting comparison is with Sr2RuO4 , which also suggests an importance of the real space structure. Although the Fermi surface is quasi-1D in Sr2RuO4 as well, the oxide has strong pieces of evidence for triplet pairing, and even when one considers the spin-singlet pairing, $`T_c`$ (for the $`\alpha ,\beta `$ bands for this material) estimated with the FLEX is very small. Due to the quasi-one dimensionality the spin susceptibility $`\chi (𝒒,\omega )`$ has linear ridges in k-space. As a result, the pair scatterings have large contributions all over the ridges, so that the (extended) $`s`$-wave gap function involves unfavorable pair scatterings across which the gap has the same sign, resulting in the reduced $`T_c`$. This makes the spin-triplet, $`p`$-wave pairing more favorable, for which the resulting gap function has large amplitudes all over the Fermi surface, except for dips at the corners . By contrast, the present case has a built-in antiferromagnetic structure within the unit cell in real space, so that the spin susceptibility has a well-defined peak around ($`\pi ,\pi `$). This makes a specific $`𝑸`$ to be relevant in the pair scattering, which makes the gap function peaked at specific points, which in turn gives rise to an enhancement in $`T_c`$.

To summarize, we should question not only the structure of the Fermi surface but also the structure of the molecular orbitals and the superconducting pairings in real space in understanding the superconductivity from the electron repulsion. One interesting tendency is that high $`T_c`$’s are found in plaquette systems where unit cells having multiple sites are connected with a relatively small inter-plaquette hopping , which also suggests the importance of the real space picture. In fact, the relation between real-space and the momentum-space pictures has been discussed for systems consisting of dimers recently. It is an interesting future problem to see how the real-space and the momentum-space pictures are related with each other for wider class of systems in realizing high $`T_c`$’s.

We wish to thank Professor J. Akimitsu for valuable discussions. This work was supported in part by a Grant-in-Aid for Scientific Research and Special Coordination Funds from the Ministry of Education of Japan. Numerical calculations were performed at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo.


format: text

text:
# University of St. Andrews Open Cluster Survey for Hot Jupiters

## 1. Introduction

The technique of discovering planets from their photometric transits enables us to probe fainter stars for planetary companions than are accessible to radial velocity surveys. This allows us to probe a wider variety of stellar environments. The technique also offers us the possibility of discovering statistically significant numbers of hot Jupiter planets from the ground within a relatively short timescale. Such a sample could be used to establish the relationships between planetary formation and key properties such as metallicity, age, etc. Janes (1996) suggested that clusters would be ideal targets for a transit search. They provide large numbers of stars within a relatively small area, the ages and metallicities of which can be determined. Radial velocity (RV) surveys have indicated that stellar metallicity in particular seems to be correlated with the presence of planets. Brown et al. (2001) reported surveying the globular cluster 47 Tuc with the Hubble Space Telescope, but interestingly found no transits. In the light of the RV results, this raises the question of whether and how this deficit is related to the low metallicity of the cluster and/or star/radiation density. To address this question, we are surveying a number of open clusters with a range of ages/metallicities. Open clusters do not suffer from the dense crowding of stars present in globulars and so can be studied with ground-based telescopes. They are also generally younger and higher in metallicity.

## 2. Observations and Data Reduction

We were awarded 30 nights of bright time on 2.5m Isaac Newton Telescope (INT) in 1999/06-07 and 2000/09. The $`0.3\mathrm{}^{}`$ field of view of the Wide Field Camera (WFC) is ideal for transit searches. We observed 3 clusters for the first 20 nights, and one cluster continuously for a further 10 nights. PSF photometry was carried out using DAOphot (Stetson 1987) while post-processing software was written in-house. Further details of these observations and our data reduction procedures can be found in Street et al. (2002) and Street (2002). The 20 nights of data on NGC 6819 have been fully reduced for all 4 CCDs, including over 38,000 stars. Figure 1 illustrates that we have reached the photometric precision required to detect hot Jupiter transits.

## 3. Transit Identification

We have developed a matched-filter algorithm which can reliably detect transit events in a reasonable time scale. Our algorithm initially generates single-transit model lightcurves for a range of transit durations and epochs. These are least-squares fitted to all lightcurves, optimising for transit depth and out-of-transit magnitude in the process. The parameters and the $`\chi _{1transit}^2`$ of best-fitting model are recorded and compared with the $`\chi _{const}^2`$ of the fit of a constant model. The transit-finder index, $`\mathrm{\Delta }\chi ^2`$, is then given by: $`\mathrm{\Delta }\chi ^2=\chi _{const}^2\chi _{1transit}^2`$.

Candidate transit events are identified by plotting $`\mathrm{\Delta }\chi ^2`$ .vs. $`\chi _{1transit}^2`$. A cut-off is established by fitting a straight line through the main body of points, then raising it by a specified number of sigma. Genuine transits will lie to the upper left of this plot. All stars lying above the cut-off are visually examined. Multi-transit model lightcurves are then fitted to the selected candidates using a range of periods, and the parameters and $`\chi ^2`$ of the new best-fitting model are recorded. The candidate selection procedure is then repeated.

## 4. Results and Future Work

In this way we have identified 7 stars showing transit-like events (eclipses of similar amplitude and duration to those of a planetary transit). The lightcurves of these stars are shown in Figure 2 and their details are given in Table 1.

Colour information indicates that in most cases, the companion object is likely to be a very low-mass star or brown dwarf - interesting objects in their own right. However, the eclipses these stars show are very similar to those of planetary transits, indicating that we can detect such events.

Our current analysis has some limitations: the algorithm is susceptible to false alarms caused by blending, stellar activity and scattered data taken in poor conditions. We are currently refining our algorithm and investigating improvements in our reduction to eliminate these effects.

## References

Brown, T. M., 2001, ApJ, 552, p. 699-709.

Janes, K., 1996, JGR, 101, iss. E6, p. 14853-14860.

Street, R. A. et al., 2002, MNRAS, 330, p. 737-754.

Street, R. A., 2002, PhD thesis., University St. Andrews.

Stetson, P. B., 1987, PASP, 99, p. 191.


text:
# Dark Energy and Matter Evolution from Lensing Tomography

## I Introduction

The weak gravitational lensing of faint galaxies weak  provides the most direct probe of mass distribution in the universe (e.g. BarSch01 ). Moreover the evolution of clustering in the mass distribution is arguably the best theoretically-grounded probe of the dark energy and dark matter deprobes . Although observations of weak lensing on large scales weaklss  are still in the discovery phase weakdet , future wide-field surveys have the potential to rival the statistical precision and cosmological utility of luminosity distance measures from supernova surveys Hut02 ; Hu01c . Even in the context of a precisely-determined homogeneous cosmology, lensing measurements are unique in that they probe the clustering properties of the dark matter and energy. These are fixed by the homogeneous cosmology only under particular assumptions of the particle constituents (e.g. cold dark matter and scalar field dark energy) CalDavSte98 ; Hu98 .

Much of the critical cosmological information lies in the temporal or radial direction. A potential obstacle for weak lensing is that the observables are inherently two-dimensional. All of the matter along the line-of-sight to a distant source contributes to the lensing. For a family of cosmological models that is described by a handful of parameters, this is not a serious drawback. Lack of radial information is largely compensated by a large angular dynamic range and external cosmological information.

Given the lack of compelling models for the dark energy and controversies surrounding the phenomenology of the dark matter on small scales, it is interesting to consider a more model-independent approach. Indeed recent studies of alternate parameterizations of the dark energy have revealed potential ambiguities in the interpretation of luminosity distance measurements WanGar01 ; MaoBruMcMSte02 ; WelAlb02 ; Teg01 ; HutSta02 . To address these issues with weak lensing, recovery of the temporal dimension becomes critical.

With future surveys that possess source photometric redshift information, recovery of the lost information is possible in principle through tomography. Photometric redshift techniques are already being applied and tested on current lensing data Witetal01 . The full two-point statistical information can be regained by cross-correlating the lensing observables on all source redshift planes Hu99 . This method utilizes both the angular clustering and the temporal evolution of the density field but obscures the nature and hence the model-dependence of the information. Additionally, the joint observables are survey dependent and computationally cumbersome to analyze.

In this paper, we instead isolate the temporal information by applying recently developed techniques to reconstruct the radial density field itself Tay01 ; HuKee02 . We will further focus solely on the linear regime where predictions are well-understood. Even utilizing only this theoretically-clean subset of information in the data, future surveys can potentially provide interesting model-independent constraints on the properties of the dark energy and matter.

The outline of the paper is as follows. In §II, we discuss the method for reconstruction and statistical forecasts. In §III, we study constraints on the growth function for a fixed homogeneous cosmology and in §IV the dark energy density evolution assuming pure cold dark matter. We discuss these results in §V.

## II Tomographic Reconstruction

We begin in §II.1 by briefly reviewing the tomographic reconstruction of the dark matter density field in a fixed background cosmology as studied in HuKee02 . We then generalize to the case where the cosmology and the density field must be jointly recovered from the data in §II.2 and review Fisher techniques for statistical forecasts in §II.3. In §II.4, we outline the fiducial cosmology and survey parameters used for illustrative purposes in the following sections.

### II.1 Known Homogeneous Cosmology

We consider the data to be the lensing convergence $`\kappa `$ in an angular pixel discretized into bins of source redshift composed into a data vector $`𝐝_\kappa `$. Weak lensing dictates that the data is a linear projection of an underlying density field plus noise in the convergence measurement $`𝐧_\kappa `$

$$𝐝_\kappa =𝐏_{\kappa \mathrm{\Delta }}𝐬_\mathrm{\Delta }+𝐧_\kappa ,$$
(1)
with

$`[𝐏_{\kappa \mathrm{\Delta }}]_{ij}=\{\begin{array}{cc}\frac{3}{2}H_0^2\mathrm{\Omega }_m\delta D_j\frac{(D_{i+1}D_j)D_j}{D_{i+1}}\hfill & \text{ }D_{i+1}>D_j\text{ },\hfill \\ 0\hfill & \text{ }D_{i+1}D_j\text{ },\hfill \end{array}`$ (2)
where $`D`$ is the comoving distance in a flat universe

$$D(z)=_0^z\frac{dz^{}}{H(z^{})},$$
(3)
with $`H^2=8\pi G\rho _{\mathrm{tot}}/3`$ defining the Hubble parameter and $`\delta D_j`$ is the width of bin $`j`$. Here and throughout subscripts on matrices are labels, whereas matrix elements are denoted as $`[]_{ij}`$. Here $`\mathrm{\Delta }=(\delta \rho /\rho )/a`$ is the density fluctuation in the bin with the growth rate in a matter-dominated universe scaled out. Note that these distances depend on the assumed cosmology. We assume for now that $`D(z)`$ has been fixed by other observations, e.g. future supernovae surveys, but relax this assumption in the following section.

The minimum variance estimator of the underlying density field is given by HuKee02

$$\widehat{𝐬}_\mathrm{\Delta }=𝐑_{\mathrm{\Delta }\kappa }𝐝_\kappa ,$$
(4)
where the reconstruction matrix

$`𝐑_{\mathrm{\Delta }\kappa }`$ $`=`$ $`𝐍_{\mathrm{\Delta }\mathrm{\Delta }}𝐏_{\kappa \mathrm{\Delta }}^t𝐍_{\kappa \kappa }^1.`$ (5)
Here

$$𝐍_{\mathrm{\Delta }\mathrm{\Delta }}=[𝐏_{\kappa \mathrm{\Delta }}^t𝐍_{\kappa \kappa }𝐏_{\kappa \mathrm{\Delta }}]^1,$$
(6)
is the noise covariance of the estimator. Note that $`𝐑_{\mathrm{\Delta }\kappa }𝐏_{\kappa \mathrm{\Delta }}=𝐈`$ so that the estimator is unbiased.

The statistical properties of the recovered density field contain cosmological information. The recovered density field is an average of the density fluctuation over a window (or mask) $`W_i(𝐱)`$ defined by the angular pixel and redshift binning. The signal covariance of these density averages

$`\left[\overline{𝐒}_{\mathrm{\Delta }\mathrm{\Delta }}\right]_{ij}`$ $`=`$ $`{\displaystyle \frac{\varphi _i}{\varphi _0}}{\displaystyle \frac{\varphi _j}{\varphi _0}}{\displaystyle \frac{d^3k}{(2\pi )^3}W_i(𝐤)W_j^{}(𝐤)P(k)},`$
where $`P(k)`$ is the linear power spectrum today, $`\varphi _i=D_{\mathrm{grow}}(z_i)/a_i`$ is the linear decay rate of the potential field, with $`D_{\mathrm{grow}}`$ the linear growth rate of the density field, normalized so that $`\varphi _i=1`$ in the matter dominated regime, and $`W_i(𝐤)`$ are the Fourier transforms of the windows. The two-point statistical properties of the reconstruction therefore contain information on the growth rate and underlying power spectrum of the density field.

### II.2 Unknown Homogeneous Cosmology

The inversion of Eqn. (4) requires an assumption of a distance-redshift relation in $`𝐏_{\kappa \mathrm{\Delta }}`$. If this relation is not fixed by external constraints, then the reconstructed density field will be a biased measure of the true density field. In this case, both the distance-redshift and growth rate must be fit to the data.

If the assumed $`D(z)`$ is close to the true $`D(z)`$, then the reconstruction of the previous section still serves as a useful representation of the data. The reconstruction matrix $`𝐑_{\mathrm{\Delta }\kappa }`$ employs a slightly incorrect assumption of the projection matrix so that $`𝐑_{\mathrm{\Delta }\kappa }𝐏_{\kappa \mathrm{\Delta }}𝐈`$. The statistical properties of the estimated density field are encapsulated in the noise matrix (6), which remains unchanged, and a new signal covariance matrix

$$𝐒_{\mathrm{\Delta }\mathrm{\Delta }}=𝐑_{\mathrm{\Delta }\kappa }𝐏_{\kappa \mathrm{\Delta }}\overline{𝐒}_{\mathrm{\Delta }\mathrm{\Delta }}𝐏_{\kappa \mathrm{\Delta }}^t𝐑_{\mathrm{\Delta }\kappa }^t.$$
(7)
The two-point statistics now also contain information about the distance-redshift relation $`D(z)`$.

### II.3 Statistical Forecasts

The two-point statistical information in the reconstruction can be exposed through the familiar Fisher approach (e.g. TegTayHea97 ). If the parameters that underly the two-point correlation are given by a vector $`𝐬_p`$, the Fisher matrix is given by

$$[𝐅_{pp}]_{ij}=\frac{N_{\mathrm{pix}}}{2}\mathrm{tr}[𝐂_{\mathrm{\Delta }\mathrm{\Delta }}^1𝐂_{\mathrm{\Delta }\mathrm{\Delta },i}𝐂_{\mathrm{\Delta }\mathrm{\Delta }}^1𝐂_{\mathrm{\Delta }\mathrm{\Delta },j}]$$
(8)
with the covariance matrix

$$𝐂_{\mathrm{\Delta }\mathrm{\Delta }}=𝐍_{\mathrm{\Delta }\mathrm{\Delta }}+𝐒_{\mathrm{\Delta }\mathrm{\Delta }}.$$
(9)
Here we have assumed that the convergence is measured in $`N_{\mathrm{pix}}`$ independent pixels. We will often write this factor as the total sky coverage in independent pixels $`N_{\mathrm{pix}}=4\pi f_{\mathrm{sky}}/A_{\mathrm{pix}}`$ where $`A_{\mathrm{pix}}`$ is the pixel area in steradians. A complete treatment would track the small correlations between neighboring pixels on contiguous patches of sky HuKee02 .

The inverse of the Fisher matrix $`𝐅_{pp}^1`$ gives an estimate of the covariance matrix $`𝐂_{pp}`$ of the measured parameters $`\widehat{𝐬}_p`$. Note that under a re-parameterization of the space, the Fisher matrix transforms as a covariant tensor

$$𝐅_{\stackrel{~}{p}\stackrel{~}{p}}=𝐉_{p\stackrel{~}{p}}^t𝐅_{pp}𝐉_{p\stackrel{~}{p}},[𝐉_{p\stackrel{~}{p}}]_{ij}\frac{p_i}{\stackrel{~}{p}_j},$$
(10)
We will use this fact to go from model-independent parameterizations of the underlying functions to model-dependent ones.

### II.4 Fiducial Model and Survey

We take as a fiducial cosmology a flat $`\mathrm{\Omega }_{\mathrm{tot}}=1`$ universe with $`\mathrm{\Omega }_c=0.3`$ in cold dark matter, $`\mathrm{\Omega }_b=0.05`$ in baryons, $`\mathrm{\Omega }_{\mathrm{DE}}=0.65`$ in dark energy; an equation of state of the dark energy of $`w(z)=1`$ corresponding to a cosmological constant; a dimensionless Hubble constant, $`h=0.65`$, scalar spectral index $`n=1`$, and amplitude of the initial curvature power spectrum $`\delta _\zeta =4.8\times 10^5`$ ($`\sigma _8=0.92`$, see Hu01c  for specific definitions of parameters). Since tomographic reconstruction will mainly be useful for next-generation surveys, it is reasonable to assume that CMB experiments will by that time have determined many of the underlying cosmological parameters to high accuracy. For simplicity we will here assume that $`\mathrm{\Omega }_ch^2`$, $`\mathrm{\Omega }_bh^2`$, $`n`$, $`\delta _\zeta `$ and $`\mathrm{\Omega }_{\mathrm{tot}}`$ are completely fixed to their fiducial values. We will return to this point in §V. These parameters fix the shape and high redshift normalization of the potential power spectrum and so we will work in the context that only the growth function and distance-redshift relation need be determined through lensing.

For definiteness, we will take the fiducial survey to be defined with circular pixels of area $`1\mathrm{deg}^2`$ so that the reconstructed density field is in the linear regime JaiSel97 . Larger pixels would reduce the number of independent pixels and hence increase the sample variance in the signal dominated regime. In the figures that follow, we take $`f_{\mathrm{sky}}=0.1`$ ($``$ 4000 deg2) but note that all lensing errors may be rescaled as $`(f_{\mathrm{sky}}/0.1)^{1/2}`$. For the redshift binning, we take $`\mathrm{\Delta }z=0.1`$ out to $`z=3`$. With this binning, the signal covariance matrix $`𝐒_{\mathrm{\Delta }\mathrm{\Delta }}`$ is nearly diagonal and the statistics reduce to the evolution of the density variance in bins.

We will assume a convergence noise spectrum of the form

$$𝐍_{\kappa \kappa }=\mathrm{diag}[\gamma _{\mathrm{rms}}^2/N_i],$$
(11)
as appropriate for random intrinsic galaxy ellipticities. We take $`\overline{n}=3.6\times 10^5`$ gal. deg-2 and $`\gamma _{\mathrm{rms}}=0.3`$ as an estimate of the usable galaxies and the shear noise per galaxy measured from a space-based platform (A. Refregier, private communication) and form the number of galaxies per bin $`N_i`$ from a redshift distribution Kai98

$$\frac{dN}{dz}\frac{dD}{dz}D\mathrm{exp}[(D/D_{})^4],$$
(12)
where $`D_{}`$ is set to reproduce a median redshift $`z_{\mathrm{med}}=1`$. These fiducial survey specifications are chosen to represent the upper range of the capabilities of surveys in the foreseeable future.

## III Growth Rate

For a fixed distance-redshift relation and high redshift power spectrum, the remaining degrees of freedom in the two-point statistical properties of the reconstructed density field are contained in the growth function $`\varphi (z)`$. To study how lensing constrains this function, we begin with a model-independent approach.

Consider the set of binned growth rates $`\varphi _i`$ as the parameters to be estimated. With fine-binning of the density field $`\mathrm{\Delta }z=0.1`$, this leads to estimates that have large correlated errors since the reconstruction effectively takes differences of noisy data. The long-time scale evolution of the density field is faithfully preserved in the reconstruction HuKee02 . Since structure grows in linear theory on the expansion time scale in gravitational instability models, this information is sufficient to constrain cosmology.

To better understand the information contained therein, consider the principle component or eigenvector decomposition of the Fisher matrix $`𝐅_{\varphi \varphi }=𝐒\mathrm{\Lambda }𝐒^t`$ and the linear combinations of the data they define

$`\widehat{𝐬}_\lambda =𝐒^t\widehat{𝐬}_\varphi ,𝐂_{\lambda \lambda }=\mathrm{\Lambda }^1.`$ (13)
In other words, the eigenvectors are the redshift representation of a new basis that is complete and yields uncorrelated, orthogonal measurements with variance given by the inverse eigenvalue $`1/\lambda _i`$. The largest eigenvalues correspond to the minimum variance directions and are shown in the upper panel of Fig. 1. The first mode has a broad single peak between $`0<z<z_{\mathrm{med}}=1`$ corresponding to the bell-shaped weight in the lensing projection of Eqn. (2). The higher modes exhibit oscillatory structure and capture information on the low order derivatives of the growth function around this intermediate redshift as well as the region $`z>z_{\mathrm{med}}`$. The spectrum of eigenvalues $`\lambda _i^{1/2}`$, scaled to represent inverse rms noise and normalized for $`f_{\mathrm{sky}}=0.1`$ are shown in the bottom panel. Although the eigenvalues reflect only the noise properties and not the signal-to-noise, the growth rate in the fiducial cosmology is sufficiently flat so that $`[𝐬_\lambda ]_i\lambda _i^{1/2}`$ shares the same form (bottom panel). Most of the signal comes from the first few eigenmodes but the signal-to-noise in the first 10-15 eigenmodes remains substantial for a large survey.

Although this principle component analysis is ideal for exposing the nature of the information, the oscillating windows makes it somewhat difficult to visualize its impact for model testing. For this purpose, it may be preferable use more localized linear combinations that retain the uncorrelated property at the expense of having overlapping (non-orthogonal) windows. Consider the linear combinations defined by rows of the “square root” of the Fisher matrix HamTeg00

$$𝐆𝐒\mathrm{\Lambda }^{1/2}𝐒^t,$$
(14)
and with a normalization chosen so that the window elements sum to unity. As shown in Fig. 2, they yield well-localized windows and provide a visualization of the data with error boxes whose width is determined as that enclosing the central $`60\%`$ of the window. We show here the errors appropriate for a survey with $`f_{\mathrm{sky}}=0.1`$.

Any model-independent constraints on the growth rate can be translated into a model-dependent one by examining the $`\chi ^2`$ of the model fit. In terms of Fisher forecasts, this is equivalent to a re-parameterization through the transformation law (10). As an example consider a family of growth functions that represent a rescaling and pivot around a fixed $`z_p`$ from the fiducial model

$$\varphi (z)=A\left(\frac{1+z}{1+z_p}\right)^p\varphi _{\mathrm{fid}}(z).$$
(15)
The pivot point can be chosen to decorrelate the errors between $`A`$ and $`p`$ by examining Fisher re-parameterization of $`z_p`$ EisHuTeg99b . This choice then has the nice property that the error $`\sigma (A)`$ in the two parameter model ($`A,p`$) are also those of the single parameter family of models ($`A`$). For the fiducial model and survey $`z_p=0.64`$, $`\sigma (A)=0.0023f_{\mathrm{sky}}^{1/2}`$, $`\sigma (p)=0.0089f_{\mathrm{sky}}^{1/2}`$. Physically, such constraints would limit the fraction $`f`$ of the dark matter in a smooth component, for example a light neutrino below its free-streaming scale. A smooth component induces a change in the growth rate in the matter dominated epoch of $`p=3f/5`$ and a consequent change in the amplitude compared with the initial conditions of $`\delta A4f`$ (e.g. HuEisTeg98 ).

We show a neutrino model with a total mass $`m_{\mathrm{tot}}=0.2`$eV distributed equally into three species BeaBel02  in Fig. 2. This test is potentially substantially more powerful than its galaxy clustering analogue due to the lack of an unknown bias HuEisTeg98 ; Elgetal02 . While most of the constraint would come from the amplitude $`A`$, information on the growth index $`p`$ is useful for distinguishing such effects from those of the dark energy (§IV) and uncertainties in the initial conditions (§V). By varying the pixel size one can in principle test the scale dependence of growth rate predicted in such models. Likewise, deviations would occur if the dark energy is not effectively smooth on the scale of the pixels. Lensing tomography offers a unique opportunity to test the clustering properties of both the dark energy and the dark matter on intermediate scales.

## IV Dark Energy Density Evolution

If the distance-redshift relation $`D(z)`$ is uncertain due to the dark energy, then these degrees of freedom must be incorporated into the statistical forecasts as well. For simplicity, we will in this context assume that the dark matter is composed solely of cold dark matter.

Fortunately for a smooth dark energy component, both the distance-redshift relation (3) and the growth rate are fixed by the dark energy density evolution $`\rho _{\mathrm{DE}}(z)`$ deprobes . Here the growth rate obeys

$`{\displaystyle \frac{d^2\varphi }{d\mathrm{ln}a^2}}`$ $`+`$ $`\left[{\displaystyle \frac{5}{2}}{\displaystyle \frac{3}{2}}w(z)\mathrm{\Omega }_{\mathrm{DE}}(z)\right]{\displaystyle \frac{d\varphi }{d\mathrm{ln}a}}`$ (16)
$`+`$ $`{\displaystyle \frac{3}{2}}[1w(z)]\mathrm{\Omega }_{\mathrm{DE}}(z)\varphi =0,`$
where the initial conditions are $`\varphi =1`$ and $`d\varphi /d\mathrm{ln}a=0`$ and the equation of state

$$w(z)\frac{p_{\mathrm{DE}}}{\rho _{\mathrm{DE}}}=\frac{1}{3}\frac{d\mathrm{ln}\rho _{\mathrm{DE}}}{d\mathrm{ln}a}1.$$
(17)
As in the case of the growth rate, we chose a model-independent parameterization as the primary representation. Consider the dark energy density in redshift bins WanGar01 ; Teg01 , specifically

$$d_i=\mathrm{ln}\left(\frac{\rho _{\mathrm{DE}}(z_i)}{\rho _{\mathrm{cr0}}}\right).$$
(18)
We choose $`\mathrm{\Delta }z`$ of the bins to be the same as the density reconstruction $`\mathrm{\Delta }z=0.1`$. Since $`\rho _{\mathrm{cr0}}`$ is the critical density today, in the limit of fine binning $`d_1=\mathrm{ln}\mathrm{\Omega }_{\mathrm{DE}}`$. We also require a finite $`w^{}(z)dw/dz`$ to ensure that the dark energy remains smooth. Therefore we choose $`d(z)`$ as a spline interpolation of $`d_i`$.

Again, the recovery of finely-binned parameters is noisy and correlated across neighboring bins. The information content is best revealed through the Fisher principle component analysis. Shown in Fig. 3 (top) are the first 5 eigenfunctions. The qualitative difference between the eigenfunctions of the growth and that of the density is the presence of substantial low redshift information in the latter. The dark energy density at low redshift affects the distance to all higher redshifts. In Fig. 3 we show the eigenvalues as $`\lambda _i^{1/2}`$ and the signal-to-noise ratio for the fiducial model and $`f_{\mathrm{sky}}=0.1`$.

It is again desirable to find an uncorrelated but more localized representation of the data. Unfortunately, the Fisher square-root technique of Eqn. (14) does not yield localized windows for the dark energy density. Instead we choose a close analogue, the Cholesky decorrelation HamTeg00  where the windows are the columns of $`𝐋`$ and $`𝐅_{dd}=\mathrm{𝐋𝐋}^t`$ again normalized to sum to unity. The windows are shown in Fig. 4 (lower panel). The windows also have the interesting property that they are strictly zero below some minimum redshift.

In Fig. 4 (upper panel) we show the projected constraints on these localized modes compared with the predictions for a $`w=0.8`$ model (points) and the actual function $`\rho _{\mathrm{DE}}(z)`$ (curve) in this model. The deviation between the points and the curve reflects the non-locality of the windows. Note that for a cosmological constant model ($`\rho _{\mathrm{DE}}=`$ const.) there is no deviation by definition (straight line) and that the expectation value of all points is $`\mathrm{ln}\mathrm{\Omega }_{\mathrm{DE}}`$. Any statistically significant difference in the values of the points in this reconstruction represent a detection of a dark energy component that is not a cosmological constant WanGar01 . This is true in spite of the non-locality of the windows and independently of the model-dependent parameterization of the dark energy.

Again one may always test specific models for the dark energy from the model-independent parameterization. Many dark energy models can be parameterized by $`\mathrm{\Omega }_{\mathrm{DE}}`$, $`w(z_{\mathrm{eff}})`$, and $`w^{}=dw/dz|_{z_{\mathrm{eff}}}`$. The pivot point $`z_{\mathrm{eff}}`$ can be chosen to be the best constrained redshift or “sweet spot” by decorrelating the errors in $`w`$ and $`w^{}`$. As was the case for the growth rate, the resulting errors on $`w(z_{\mathrm{eff}})`$ are the same as in the case of a two parameter model $`(\mathrm{\Omega }_{\mathrm{DE}},w)`$. For the fiducial model and survey, this is $`z_{\mathrm{eff}}=0.33`$ and $`\sigma (w)=0.016f_{\mathrm{sky}}^{1/2}`$. Note that this constraint is marginalized over $`(\mathrm{\Omega }_{\mathrm{DE}},w^{})`$ without prior assumptions to their values.

As in the case of supernovae luminosity distance measures WelAlb02 ; MaoBruMcMSte02 , there remains a degeneracy between $`\mathrm{\Omega }_{\mathrm{DE}}`$ and $`w^{}`$ in that they may both be adjusted upward to keep the dark energy density at the well-constrained low redshifts fixed. In Fig. 5 shows, we show the $`68\%`$ confidence region with various assumptions of prior knowledge on $`\mathrm{\Omega }_{\mathrm{DE}}`$: none, $`\sigma (\mathrm{\Omega }_{\mathrm{DE}})=0.03`$, $`=0.01`$. With $`\mathrm{\Omega }_{\mathrm{DE}}`$ fixed, the errors become $`\sigma (w^{})=0.046f_{\mathrm{sky}}^{1/2}`$. With no prior, the errors degrade to $`\sigma (w^{})=0.069f_{\mathrm{sky}}^{1/2}`$.

Conversely if $`w^{}`$ is fixed, the constraints on the dark energy density improve to $`\sigma (\mathrm{\Omega }_{\mathrm{DE}})=0.003f_{\mathrm{sky}}^{1/2}`$. This two parameter family ($`\mathrm{\Omega }_{\mathrm{DE}},w`$) represents the amplitude and slope of the dark energy density itself at a normalization point of $`z=0`$. The remaining degeneracy with $`w`$ can be removed by again going to the “sweet spot”, here $`z=0.13`$ where the errors on the dark energy density improve by a factor of 1.5. Again this also represents the errors on $`\mathrm{\Omega }_{\mathrm{DE}}`$ in the single parameter family of dark energy models $`\sigma (\mathrm{\Omega }_{\mathrm{DE}})=0.002f_{\mathrm{sky}}^{1/2}`$. Errors on $`w`$ of course remain unchanged.

## V Discussion

We have shown that lensing surveys covering more than a few percent of the sky with good photometric redshift information can probe the time evolution of the linear growth function and distance-redshift relation, both of which are sensitive to properties of the dark energy and dark matter. Specifically, we have tested a model-independent parameterization of the linear growth rate and/or dark energy density discretized into bins in redshift. Deviations in the growth rate would indicate a component of the dark matter that is not effectively cold or dark energy that is not smooth on the lensing scale. Deviations in the constancy of the dark energy density would rule out a cosmological constant model.

In this exploratory study, we have made several simplifying assumptions that would need to be addressed in a concrete implementation. Perhaps the primary one is that future CMB measurements will completely fix the high redshift cosmology. The most uncertain piece involves the amplitude of the initial fluctuations on the scales relevant to the lensing pixels, $`k0.05`$ Mpc-1 for degree scales. Fortunately, the pivot point of CMB anisotropy experiments with several arcminute scale resolution is sufficiently close to the lensing scale that the slope and shape of the initial power spectrum do not cause much ambiguity Hu01c .

To fully utilize the lensing information, the initial amplitude must be fixed to an accuracy better than the amplitude of the growth function, which we have found to be $`0.002f_{\mathrm{sky}}^{1/2}`$, i.e. percent level accuracy for surveys of several thousand square degrees. For CMB anisotropies, this precision requires that the optical depth during reionization must be determined to $`\sigma (\tau )0.01`$ to resolve the amplitude degeneracy. If determined from CMB observations alone, this will require polarization measurements with a precision comparable to the Planck satellite planck , which can in principle achieve $`\sigma (\mathrm{ln}\delta _\zeta )=0.0044`$ at $`k=0.05`$ Mpc-1 Hu01c . Direct measurements of the reionization epoch can also resolve the ambiguity Hu01c . Even in the absence of this information, the evolution of the growth and luminosity-distance relation are still constrained. These issues are best addressed through joint parameter estimation.

On the lensing side, the most important assumption is that the noise in the convergence measurements is well-calibrated and not significantly larger than the projections based on intrinsic ellipticities. Even aside from the demanding requirements for control of systematic errors, there may be intrinsic correlations in the ellipticities intrinsic  that need to be modeled or avoided by increasing the pixel scale and redshift bin widths. The recovered information is largely insensitive to the redshift bin width since the high signal-to-noise modes are all low frequency. Errors scale roughly as $`A_{\mathrm{pix}}^{1/2}`$ due to the loss of independent modes in a fixed survey area.

We have also neglected sample covariance between the pixels but note that we have correspondingly neglected the information contained in such correlations. Indeed, we have completely neglected the information contained in the non-linear regime which in fact contains the majority of the information from lensing tomography HuKee02 . Clearly, future studies will be required to see how best to mine the model-independent information contained in lensing tomography.

Acknowledgments: I thank D. Huterer and C.R. Keeton for useful conversations and E. Linder for pointing out a typo in the draft. This version includes errata from a coding bug, affecting constraints involving the present dark energy density, pointed out by K. Abazajian and S. Dodelson. WH is supported by NASA NAG5-10840 and the DOE OJI program.


text:
# Instability Heating of Sympathetically-Cooled Ions in a Linear Paul Trap

## I introduction

Laser cooling  has allowed atomic physicists to study gas-phase atomic systems at ultra-low temperatures, largely removing thermal effects and allowing many atomic systems to be studied at their most fundamental quantum level (e.g. atom interferometry ). However, although a range of neutral and ionic atomic systems has been cooled to sub-Kelvin temperatures with direct laser cooling , there is a limitation of this approach. Conventional laser cooling techniques, both Doppler  and sub-Doppler , require isolated two- or three-level systems. If a system has additional low-energy meta-stable states or multiple nearly degenerate ground states (e.g. rotational structure in molecular ions) then the laser-cooling cycle terminates when the system falls into a non-laser-cooled level. The solution to this problem generally requires one additional re-pumping laser system for each meta-stable state (e.g. Ref. ). Cooling of molecules is especially challenging due to their large number of accessible rotational states. Understandably, there only exist a small number of very complicated theoretical proposals (e.g. Ref. ) to directly cool molecular ions. Therefore, many atomic systems such as those lacking the necessary energy level structure and virtually all molecular systems cannot be directly laser cooled. However, a much wider range of atomic and molecular systems can be taken to very low temperatures through a process known as sympathetic cooling.

In this technique, a collection of two interacting species, the sample species (which can not be directly laser cooled) and the laser-cooled species, are confined to the same region of space and isolated from their surroundings through the use of neutral particle or ion traps. As the laser-cooled species interacts with the light field its temperature falls. When these cold atoms interact or collide with the sample species, they extract energy from the sample species thus sympathetically reducing the sample species temperature. Recently, a variation of this approach has been applied to neutrals whereby evaporatively-cooled atoms sympathetically cool a sample species to produce Quantum Degenerate Gases from species whose self-interaction (the interaction between one atom with another identical atom) is too weak to permit efficient evaporative cooling on its own . This type of sympathetic cooling in neutrals depends fundamentally on the short range interactions (collisions) between the directly-cooled and sample species, and thus the effectiveness of sympathetic cooling depends fundamentally on the particular pair of species selected and their associated inter-atomic or atom-molecule interaction potentials. This makes it rather difficult to produce simple and broadly general rules regarding the effectiveness of sympathetic laser cooling that would apply to all neutral atomic and molecular species. However, when considering sympathetic cooling between ions, the dominant interaction force is Coulomb interactions, which are both long range and quite strong. This means that ions do not tend to approach each other close enough that the short-range interaction potentials, which differ from atom to atom and molecule to molecule, become important. Therefore, the effectiveness of sympathetic cooling between charged particles tends to be a function of only the mass and charge of each of the two ions, allowing for rather generalised rules regarding sympathetic cooling of charged particles to be deduced.

This work marks the beginning of a theoretical and experimental program to study the processes associated with the sympathetic cooling of ions. The longer term goals are to develop the tools and expertise to produce low-temperature, gas-phase samples in as wide a range of atomic and molecular species as is physically possible for use in such areas as molecular spectroscopy, spectroscopy of stable and unstable atomic isotopes, and the physical implementation of Quantum Information Science protocols. The following pages outline a computational/theoretical study of the apparently simple question: What is the range of ion masses that a particular laser cooled ionic species can sympathetically cool in a quadrupole ion trap?

Sympathetic cooling of molecular ions by laser-cooled atomic ions within a linear Paul trap has been experimentally demonstrated  using ions of similar mass; reaching translational temperatures below 100 mK. Further experimental and molecular dynamics (MD) simulations have shown axial translational temperatures as low as 10 mK using a single ion species to sympathetically cool over a wider mass range , but no particular or general limits on this range were determined. More recently, work related to a novel ion trap in-situ mass spectrometry technique under development by the Baba and Waki group  has lead to an MD simulation study, which illustrated that for particular trap conditions there is a lower bound on the sample ion mass ($`m_s`$) that can be sympathetically cooled. This limit is 0.54 times the mass of the laser cooled ion ($`m_c`$. However, no upper bound on the masses that could be sympathetically cooled was examined. This work investigates both the cooling of higher mass ions, indicating a lack of a concrete upper bound, and examines the question of sympathetic cooling of light sample ions in a more general setting. The latter point lead to techniques to modify the cooling threshold value, and provided a general theoretical framework for the heating processes in terms of the concept of instability heating.

## II Sympathetic Cooling of Trapped Ions: The Basics

As mentioned in the previous section, sympathetic cooling of ions operates through Coulomb interactions with laser cooled ions. Theoretically, such interactions between ions that are completely isolated from their surroundings, other than via the laser-cooling beam, will eventually bring all of the ions into thermal equilibrium at the laser-cooled ion temperature. However, trapped ions are not completely isolated because the electric fields (and magnetic fields for Penning Traps) used to suspend the ions in space can themselves produce a number of heating mechanisms (e.g. rf-heating ) and it is the balance of laser cooling and trap heating that leads to an equilibrium temperature. Therefore, depending on the relative efficiency of energy extraction through collisions with laser cooled ions and trap related heating effects a given ion can either increase or decrease its energy. It is this balance that leads to the heating-cooling threshold observed in Ref. . In this work, we will concentrate on ions stored in a linear-geometry, rf quadrupole-electric field trap (commonly referred to as a linear Paul trap). In this type of trap, the ions are contained in the z-direction by a DC trapping field, while a rf-oscillating saddle-shaped potential traps the ions in the xy-plane .

The stability of a single charged particle in an oscillating electric quadrupole field is discussed in many references . Highlighting the relevant points, we start with a single ion in a linear Paul Trap, where the confining potential along the z-axis of the quadrupole is approximated as harmonic. The ion is subject to the following potential in x and y coordinates

$$\mathrm{\Phi }(x,y)=\frac{U_{rf}\left(y^2x^2\right)\mathrm{cos}(\mathrm{\Omega }t)}{2r_0^2},$$
(1)
where $`U_{rf}`$ and $`\mathrm{\Omega }`$ are the amplitude and angular frequency of the applied rf-electric potential, respectively, and $`r_0`$ is the distance from trap centre to one of the electrodes. Given that $`𝐅=e\mathrm{\Phi }`$, the equations of motion for an ion of mass $`m`$ and charge $`e`$ in the x-y plane are given by

$`\ddot{x}+(e/mr_0^2)U_{rf}\mathrm{cos}(\mathrm{\Omega }t)x`$ $`=`$ $`0,`$ (2)
$`\ddot{y}(e/mr_0^2)U_{rf}\mathrm{cos}(\mathrm{\Omega }t)y`$ $`=`$ $`0.`$ (3)
Recasting these equations in unitless parameter form $`(\xi =\mathrm{\Omega }t/2,q_u=q_x=q_y=2eU_{rf}/m\mathrm{\Omega }^2r_0^2)`$ we get

$$\frac{^2u}{\xi ^2}[2q_u\mathrm{cos}(2\xi )]u=0,$$
(4)
where $`u`$ represents either x or y. Equation (4) is the Mathieu equation in its canonical form , its general solution being

$$u=\alpha ^{}e^{\mu \xi }\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}C_{2n}e^{2in\xi }+\alpha ^{\prime \prime }e^{\mu \xi }\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}C_{2n}e^{2in\xi }$$
(5)
where $`\alpha ^{}`$ and $`\alpha ^{\prime \prime }`$ are integration constants determined by the initial conditions, (i.e. $`u(\xi _0)`$, $`\dot{u}(\xi _0)`$, and $`\xi _0`$), while $`C_{2n}`$ and $`\mu `$ depend solely on the value of $`q_u`$ and not the initial conditions . If $`\mu `$ is purely imaginary then the ion will have an oscillatory motion (trapped particle). However, if $`\mu `$ is real or complex the motion will be unbound (unstable particle). The determination of $`\mu `$ from $`q_u`$ in general is a difficult problem and is best handled numerically. It is well known that if $`|q_u|<0.908`$, then $`\mu `$ is purely imaginary and if $`|q_u|>0.908`$ $`\mu `$ becomes complex, other than for some very localised values of $`q_u`$ . For this reason $`q_u`$ is called the stability parameter.

For stable trapping parameters, the motion of an ion in x and y, as described by Eq. (5), can be broken down into a superposition of the ion’s micromotion and its secular motion. The micromotion is an oscillatory motion at the angular frequency of the trapping field that increases in amplitude the further the ion is from the trap centre. It is the direct result of the oscillation of the trapping field driving the ion back and forth as the field reverses. The lower frequency secular motion is a simple harmonic oscillation in the time-averaged pseudo-potential generated by the trapping field. The secular frequency of this motion, for small values of $`q_u`$, is given by

$$\omega _{\mathrm{sec}}=\frac{q_u}{\sqrt{8}}\mathrm{\Omega },$$
(6)
while the amplitude can be expressed in terms of the temperature based on the assumption of $`\frac{1}{2}k_BT`$ of energy in the relevant mode:

$$u_0=\frac{4}{q_u\mathrm{\Omega }}\sqrt{\frac{k_BT}{m}}.$$
(7)
Sympathetic cooling requires that multiple interacting ions are simultaneously stored in the trap. We cannot deal with such a many-body system analytically, but rather it is possible to calculate the motion of each ion by integrating the equations of motion for each ion in the trap. Additional terms have to be added to Eq. (3) to account for the Coulomb repulsion between the ions, as well as the laser cooling . The equation of motion for z is similar to the x and y equations, although it is simpler because the electric field in this direction is DC.

Our simulation numerically integrated the motion of 35 laser cooled Mg<sup>+</sup> ions interacting with 5 sample ions using the Euler-Picard predictor-corrector method with 0.1 to 1 ns time-steps. The ions are confined radially by the potential given by Eq. (1) and axially by a harmonic well with a force constant of 1 mV/mm. The initial conditions for the simulations started the ions in random positions with the sample ions at a temperature of 40 K and the laser cooled ions at a temperature of 15 K. The laser intensity was selected so that in the absence of sample ions, rf-heating balanced the laser cooling and the laser-cooled ion temperature would remain the same. In addition, the laser intensity was maintained at a low enough level to prevent the formation of Wigner or ion crystals. These conditions were chosen to be equivalent to those used in Ref. . The simulation output provided the average sample ion and laser-cooled ion temperatures as a function of time. The temperature was calculated by averaging the energy of each ion over 4 cycles of secular motion and then calculating the average ion energy for the 5 sample ions and the 35 laser cooled ions, independently. The ion temperatures were determined from these average energies using

$$E=\frac{5}{2}k_BT.$$
(8)
Integration was generally carried out for 30 ms, which was usually sufficient to determine whether the temperature of the sample ions was increasing or decreasing with time. First, the laser-cooled ions were interacted with a range of more massive sample ions (up to 8 times the mass of the laser-cooled ions). This confirmed that sample ions heavier than the laser-cooled ions generally cooled. However, as the mass difference increased, the rate of cooling decreased.

Next the case in which the sample ions are less massive than the laser-cooled ions was examined. For a given $`U_{rf}`$ and $`\mathrm{\Omega }`$, heating or cooling of the sample ion was observed to depend on its mass. We confirmed that for the same trap conditions as in Ref. , there is the same lower bound on the sample ion mass that can be sympathetically cooled, i.e. $`m_s=0.54m_c`$. At this point, tests were made to determine if sample ions with $`m_s<0.54m_c`$ could be cooled by changing the trap parameters. It was found that there always exist values of $`q_u`$ for which this is possible. A large number of simulations for a range of sample ions with mass below $`0.54m_c`$ were carried out. In all cases the results were more or less the same as those illustrated in Figs. 1 and 2. Figure 1 is for ($`m_s=8`$ amu, $`m_c=24`$ amu) and Fig. 2 is for ($`m_s=12`$ amu, $`m_c=24`$ amu). Both of these figures show the typical result that if the value of $`q_u`$ is reduced sufficiently, the sample ions can be cooled. To further test this hypothesis, maps of the $`U_{rf}`$-$`\mathrm{\Omega }`$ space were made to determine the trapping field parameters for which the sample ions heated, cooled, or remained basically unchanged in temperature. Figure 3 shows a pair of typical maps. The $`q_u=0.908`$ curve indicates the single sample ion trapping stability threshold. As can be seen, sympathetic cooling does not occur close to this threshold and thus the ions are heated. As we move away from this threshold we enter a region where there appears to be a balance thus the ions neither heat nor cool. Finally, for low enough values of $`q_u`$ sympathetic cooling begins to occur. This begs the question of why does cooling occur away from the single ion stability threshold and can we predict the threshold for sympathetic cooling? The answer to these questions are given in the next two sections.

## III Instability Heating Theory

Although we have not yet discussed the precise heating mechanism responsible, the results summarised in Fig. 3 clearly suggest the existence of a second mass threshold for the behaviour of trapped ions. The first mass threshold, dictated by the $`|q_u|<0.908`$ trap stability criterion, determines if the ion stays in the trap. The second higher-mass threshold, which is also a function of the trapping field parameters, determines whether or not the sample ion can be sympathetically cooled.

From the results of the previous section, it is reasonable to conclude that the two mass thresholds have similar functional dependences on the trapping field parameters and ion mass, suggesting that the physical causes of these two thresholds might be related. These considerations prompted us to investigate if the heating process is related to transient instabilities in the ion motion due to interactions with the other ions in the trap, an effect not considered in the single ion treatment. We will refer to this process as instability heating.

The basic physical mechanism for instability heating is as follows. Equation (5) gives the general form for motion of an ion in an oscillating quadrupole field. For a single trapped ion, $`q_u`$ is a constant and for $`|q_u|<0.908`$, $`\mu `$ is purely imaginary thus the orbit is oscillatory. When $`|q_u|>0.908`$, $`\mu `$ will contain a real part which will cause the trajectory of the ion to grow exponentially. However, when the ion in question is in the presence of other ions, $`q_u`$ can no longer be considered a constant but the additional randomly varying Coulomb forces require it to be replaced with a time-varying effective stability parameter $`q_{\mathrm{eff},u}`$. Viewed as such, $`|q_{\mathrm{eff},u}|`$ can now exceed the threshold value of 0.908 for short periods of times without the ion escaping the trap. These time intervals in which the ion trajectory grows exponentially are followed by intervals for which $`|q_{\mathrm{eff},u}|`$ drops below 0.908 thus re-stabilizing the ion at a larger maximum radius and a higher energy orbit.

To calculate $`q_{\mathrm{eff},u}`$, we will consider our physical system at a particular point in time and replace the sample ion of interest, which is exposed to the trap forces ($`F_{T,u}`$) and Coulomb force due to other ions ($`F_{i,u}`$) present in the trap, with an equivalent ion. Furthermore, we will assume that this equivalent ion is positioned at the same point in the trap as the ion of interest, having the same velocity and acceleration, and is subject to the same trap force but is not exposed to the Coulomb force due to the other ions, thus requiring it to have a time-varying effective mass. The effective mass of the equivalent ion in terms of the mass of the sample ion is given by

$$m_{\mathrm{eff},u}=\frac{m}{\left(1+\frac{F_{i,u}}{F_{T,u}}\right)}.$$
(9)
Since the effective mass is the mass for an equivalent single ion in a trap, Eq. (4) can be used to obtain $`q_{\mathrm{eff},u}`$ given by

$$q_{\mathrm{eff},u}=q_u\left(1+\epsilon _u\right),$$
(10)
where

$$\epsilon _u=\frac{F_{i,u}}{F_{T,u}}.$$
(11)
Note that $`q_{\mathrm{eff},u}`$ is both time varying and can have different values for the x and y directions. We can now write the instantaneous stability criterion as

$$|q_{\mathrm{eff},u}|=\left|q_u\left(1+\epsilon _u\right)\right|0.908.$$
(12)
When $`\left|q_{\mathrm{eff},u}\right|>0.908`$ (i.e. when the inter-ionic forces sufficiently exceed the trapping forces) we would expect the possibility of exponential growth of ion trajectory, leading to increased ion temperatures. Figures 4a and 4b illustrate that this is indeed the case. These simulations show a short snapshot of the temporal evolution of the x-coordinate of a sample ion in parallel with the computed value of $`q_{\mathrm{eff},x}`$. They were obtained with a single sample ion, a single laser-cooled ion, and with $`q_u`$ ($`=0.9`$) set very near the single ion instability threshold in order to maximize the frequency of instability points. The evolution of the particle motion shows a sequence of steps in the amplitude of the ion’s motion, corresponding to steps in the temperature of the ion. For each of these steps, we observe a corresponding spike in $`q_{\mathrm{eff},x}`$, where $`|q_{\mathrm{eff},x}|`$ exceeds $`0.908`$. In some cases, the amplitude steps down, which is expected since Eq. (5) contains both exponentially growing and decaying terms. The energy drops only when the decaying term is much larger than the growth term, while it increases under all other conditions. On the other hand, there are definitely many points where $`q_{\mathrm{eff},x}`$ exceeds $`0.908`$, but no significant change in the motional amplitude occurs. The explanation for this curious behaviour is given below.

What is perhaps more important than $`q_{\mathrm{eff},u}`$ merely exceeding 0.908 is the time interval over which $`q_{\mathrm{eff},u}`$ has an unstable value. As one can see from Eq. (11), even in the presence of a distant laser cooled ion $`\epsilon _u`$ undergoes a singularity each time the trapping force on the sample ion becomes zero. This occurs at least twice per cycle of the rf-trapping potential. However, when the time period over which $`q_{\mathrm{eff},u}`$ takes an unstable value is short compared to the time interval for the rest of the cycle, the ion trajectory does not significantly grow, hence no measureable heating occurs. Since the ion motion evolves at the secular frequency, a singularity lasting for such a short time interval should not significantly effect the energy of the sample ion or its temperature.

More significantly, $`\epsilon _x`$ (or $`\epsilon _y`$) also undergoes a singularity as the ion crosses the y (or x) axis. This occurs for each half cycle of secular motion. The important difference between this and the trapping field singularities is that the instability periods can last over a much longer time period as illustrated in Fig. 5. Here, the solid curve represents the trap force on a moving ion as a function of time and the horizontal dotted lines represent the magnitude of the slowly varying Coulomb force due to another ion. For times less than 0.8 $`\mu `$s, the amplitude of the trap force is large (i.e. the sample ion is not near the y axis) and thus the inter-ionic force dominates only for very short time-intervals as illustrated by A, which is too short to significantly effect the trajectory. However, when the ion is close to the y axis (i.e. times near 1 $`\mu `$s), the trapping force is weaker and the inter-ionic force dominates over a longer time period, possibly over several cycles of the trapping field if the ion spacing is small. This is illustrated by event B. This results in a measureable heating of the sample ion. The importance of the period over which the ion remains unstable is shown in part (c) of Fig. 4. We see that the quantized steps in the amplitude of the ion motion occur when $`q_{\mathrm{eff},x}`$ exceeds $`0.908`$ and does so for a period exceeding 30 to 40 ns. To illustrate this point, three events are highlighted. At point A, $`q_{\mathrm{eff},x}`$ spikes to well above 1 for well over 100 ns, resulting in a clear step in the motional amplitude. At point B, there is also an amplitude step, but here $`q_{\mathrm{eff},x}`$ barely exceeds 1. However, $`q_{\mathrm{eff},x}`$ stays above the stability threshold for a relatively long time ($`>100`$ ns). On the other hand, at point C there are several spikes all exceeding 1. However, none of the spikes last longer than 20 ns, and the ion motion shows no significant change.

## IV Calculation of Heating-Cooling Threshold

To generate a predicted shape for the cooling-heating threshold, we need to calculate a balance between the heating and cooling rates. To calculate the heating rate, we first need the variation of the mean instability time with the trapping field parameters. From Eq. (11) and 12, the time interval for which $`|q_{\mathrm{eff},u}|`$ is larger than 0.908 is given by the time interval when

$$|F_{T,u}|\frac{|F_{i,u}|q_u}{0.908\pm q_u}.$$
(13)
is satisfied. The upper sign is used when $`F_{T,u}`$ and $`F_{i,u}`$ have opposite sign and the lower sign when they are pointed in the same direction. $`F_{i,u}`$ is the average inter-ionic force during a collision. An estimate of this time interval can be obtained as a function of trapping field parameters. The trapping force can be determined by taking the gradient of Eq. (1). Using Eqs. 6 and 7, the force on an ion near the trap centre can be approximated by (see Fig. 5)

$$F_{T,x}=\frac{ex}{r_0^2}U_{rf}\mathrm{sin}\mathrm{\Omega }t^{}=\frac{eu_0\mathrm{sin}\omega _{\mathrm{sec}}t}{r_0^2}U_{rf}\mathrm{sin}\mathrm{\Omega }t^{},$$
(14)
where $`t`$ and $`t^{}`$ are used to indicate that the secular motion and micromotion may be out of phase. We neglect the small effect of micromotion on the x-coordinate of the ion near the trap axis. Given that the time scale of interest is the secular period, much longer than the trapping field period, the trapping field oscillation can be averaged out. Strictly speaking this is zero, but stability is dependent on $`|q_{\mathrm{eff},u}|0.908`$. Thus it is really the absolute value of the forces, not their signs which is of interest. Therefore, it is reasonable to replace $`\mathrm{sin}\mathrm{\Omega }t^{}`$ with its rms value. Assuming that we are only interested in times close to the point where the ion crosses the trap axis, the small angle approximation $`\omega _{\mathrm{sec}}t1`$ can be used. Thus

$$F_{T,x}=\frac{eu_0\omega _{\mathrm{sec}}t}{2r_0^2}U_{rf}.$$
(15)
Using Eqs. 6 and 7 , we find

$$F_{T,x}=\frac{eU_{rf}}{2r_0^2}\sqrt{\frac{2k_BT}{m_s}}t.$$
(16)
Substituting Eq. (16) into Eq. (13), and solving for the instability time $`\tau `$, we find:

$$\tau =\frac{0.908|F_{i,u}|q_u}{0.908^2q_u^2}\frac{r_0^2}{eU_{rf}}\sqrt{\frac{m_s}{2k_BT}}.$$
(17)
One can see immediately that the instability increases in duration as $`q_u`$ approaches 0.908, and eventually becomes infinity at the point that $`q_u`$ reaches the single ion instability threshold.

Given Eq. (17), we can calculate the heating rate by also considering the collision rate, the frequency at which these axis-crossing instabilities happen, and the rate at which the energy grows during the instability. The collision rate is simply a function of the number of trapped ions and trapping volume, while axis crossings occur at twice the secular frequency. To determine the rate of energy growth during the instability, we must calculate how rapidly the radius of the ion trajectory grows. The parameter $`\mu `$ in Eq. (5) governs this growth. Since $`q_u=0.908`$ marks the transition from complex to real values for $`\mu `$, this coefficient, at least near the transition point, should take on the form

$$\mu =\eta \sqrt{|q_u|0.908},$$
(18)
where $`\eta `$ is a constant of order unity. This cannot be shown analytically for the Paul trap. However, for the closely related rotating saddle trap that can be solved analytically, $`\mu `$ takes the form of Eq. (18. Converting back from unitless parameters, the rate of exponential growth in the motion during the instability is given by

$$\frac{1}{2}\mu \mathrm{\Omega }=\frac{1}{2}\eta \mathrm{\Omega }\sqrt{|q_{\mathrm{eff},u}|0.908},$$
(19)
while the rate of energy increase is simply twice this amount due to the harmonic nature of the potential. Here, $`q_u`$ has been replaced with $`q_{\mathrm{eff},u}`$ since it is the quantity that reaches unstable values. We will approximate Eq. (19) using a mean value for $`q_{\mathrm{eff},u}`$ by assuming that $`F_{i,u}`$ is roughly constant over the instability period and that $`F_{T,u}`$ can be replaced by the average of its maximum magnitude during instability and its minimum value of zero (see Eq. (13)). When $`q_{\mathrm{eff},u}`$ for parallel and anti-parallel forces are averaged, the result is independent of $`q_u`$ leaving an average value for the exponential energy growth rate of

$$\mu \mathrm{\Omega }=\eta \mathrm{\Omega }\sqrt{0.908}.$$
(20)
Finally, taking the heating rate to be proportional to the the instability time, energy growth rate, the collision rate ($`\sigma _{\mathrm{coll}}`$), and the axis crossing frequency ($`2\sigma _{sec}`$) yields

$$R_h=\frac{K_h\sigma _{\mathrm{coll}}0.908|F_{i,u}|q_u}{0.908^2q_u^2}\frac{\eta }{4\pi }\sqrt{\frac{0.908}{m_sk_BT}}.$$
(21)
To avoid overall heating of the ions, $`R_h`$ must be balanced by the cooling effect, which occurs via elastic collisions between slow and fast moving particles of different masses. Simple classical mechanics arguments show that the efficiency of elastic collision energy transfer decreases as the mass difference between the particles increases. The energy transfer efficiency is given by $`4\rho /(1+\rho )^2`$, where $`\rho =m_s/m_c`$ . Taking the cooling rate to be proportional to the product of this efficiency and the collision rate gives

$$R_c=K_c\sigma _{\mathrm{coll}}\frac{m_s/m_c}{\left(1+\frac{m_s}{m_c}\right)^2}.$$
(22)
When the two rates balance out, we should see a transition from heating to cooling of the sample ions. Equating Eqs. 21 and 22, we obtain

$$\frac{q_u}{0.908^2q_u^2}=\frac{4\pi K}{0.908|F_{i,u}|\eta }\sqrt{\frac{k_BT}{0.908}}\frac{\sqrt{m_s^3}/m_c}{(1+m_s/m_c)^2},$$
(23)
where the constant K, which is independent of the trapping parameters, is defined to be $`K_c/K_h`$. The collision rates in Eqs. 20 and 21 have been taken to be the same. This is valid if the number of laser-cooled ions significantly exceeds the number of sample ions. In general, cooling occurs when sample ions collide with laser-cooled ions, while instability heating is caused by interactions between any two ions. However, if most of the ions in the trap are laser-cooled ions, then most of the instability heating is a result of interactions between sample ions and laser-cooled ions. This means that both the heating and cooling rates are proportional to the same collision rate, and thus this rate cancels out. Therefore, the rate of collisions determines the rate of heating or cooling and not the sign of the temperature change.

Equation (23) provides the theoretical justification that lines of constant $`q_u`$ define the heating-cooling boundaries in U<sub>rf</sub>-$`\mathrm{\Omega }`$ space. Although we cannot calculate threshold values for $`q_u`$ from first principles, because of the unknown scaling factors, it can be done numerically. This was done to produce the $`q_u=0.509`$ curve in Fig. 3b, showing that this line effectively separates the regions of heating and cooling. Once the threshold value of $`q_u`$ is determined for one mass combination, Eq. (23) can be used to determine the value for other mass combinations. This is illustrated by the $`q_u=0.408`$ line in Fig. 3a. This line was not fit to the data, but rather it was calculated from the constant value for Fig. 3b. It should be noted that this model has been applied only to a limited mass range and that perhaps other mass effects have to be included if very large mass differences are to be considered.

Up to this point, we have concentrated mainly on the effect of the mass of the sample ion. However, the laser-cooled ion mass is also significant. This is because the efficiency of energy transfer between two colliding particles decreases as the ratio of the masses deviates from unity, becoming effectively zero when one mass is negligible as compared to the other. However, the instability-heating rate is expected to be relatively unchanged as it results from the Coulomb interactions between the ions and with the trapping field. Thus for larger mass ratios, the heating-cooling threshold should move further from the single-ion instability threshold in order to reduce the amount of instability heating that must be overcome. This reduces the overall collision rate which decreases the rate of heating or cooling of the sample ion. The effect of laser-cooled ion mass was examined by running simulations with 18 amu sample ions interacting with either 24 amu or 36 amu laser cooled ions for a range of $`q_u`$. The results are summarized in Table I which illustrates the expected slight shift in the cooling-heating threshold for different laser-cooled ion masses. These results are in good agreement with the heating-cooling transition values predicted by Eq. (23) based on the 12 amu/24 amu threshold. The predicted heating-cooling threshold values are 0.585 for 18 amu/24 amu and 0.56 for 18 amu/36 amu.

## V Discussion and Future Directions

The results of the preceding sections provide us with a physical insight into the competing processes that lead to sympathetic heating or cooling of ion mixtures. Generally speaking, we can conclude that if sympathetic heating is occurring for light sample ions, the trap parameters can be adjusted to move away from the single ion instability threshold, thereby allowing the sample ions to cool. However, there is a limit to this effect. As one moves far away from the instability threshold, as would be needed to cool very light ions, the trapping field becomes weaker, the trapping volume grows, and the collision rate becomes very small. Although an ion still cools, the rate of cooling may be so small as to be of little practical use.

Instability heating is to be distinguished from rf-heating . This transient effect is caused by random interactions between ions, temporarily driving them into unstable motion. In contrast, rf-heating of theions near the quadrupole axis is caused by continuous Coulomb interactions with the ions which are constantly driven by rf field at the periphery. Rf-heating cannot be suppressed solely by changing the trap parameters, only by reducing the numbers of ions and cooling the sample to the point where all the ions crystallize along the central axis of the trap where micromotion is not present. However, instability heating, which is a major source of heating in the low-density gas phase regime, can be quenched by choosing favourable trap parameters so that the ions are less susceptible to heating. This is achieved by minimising the duration of potentially unstable motion and the frequency at which these potential instabilities occur, so that in the event of a collision the ions is more likely to decrease its energy.

Further theoretical work is needed to study the effect of changing the relative numbers of sample and laser cooled ions. When the number of sample ions increases, the interactions between them can no longer be ignored. This will ultimately cause more heating, since such interactions cannot result in a decrease in the mean energy. Although, Eq. (16) suggests that instability heating should decrease in the higher energy regime this has yet to be tested. Such collisions at high temperatures ($``$300K) require much more computational time and other techniques may be developed to model them. At low energies ($``$1K) Eq. (16), which is applicable for ions in the gas phase, becomes less valid since the ions are in the crystalline phase. Instability heating is expected to be insignificant here since the axial separation is always large enough to prevent the transverse forces, which cause instabilities, from becoming too large. The dominant heating mechanism in this low energy crystalline phase would be rf-heating when the number of ions is large.

In this work, we have ignored the internal degrees of freedom of the sample ions, which is reasonable if the ions are all atomic, but less so in the case of molecular sample ions. The distinction between atomic sample ions and molecular sample ions was deliberately overlooked, since the focus of this work is the translational degrees of freedom. In ignoring the molecular ions internal structure, we have a first order approximation as to whether or not the sample molecular ions will translationally heat or cool. Translational cooling in this approximation is seen as a significant component of a complete model of the sympathetic cooling. This is because a major energy transfer process is expected to be intra-molecular transfer of energy from rotational and vibrational modes to translational degrees of freedom followed by inter-ionic translational energy transfer to the laser-cooled atoms, as modeled in this work.

## VI Conclusions

The mechanism through which sample ions are heated in a linear Paul trap is a direct result of collisions, the same collisions that would otherwise bring about thermal equilibrium at a lower temperature. A collision can cause an ion to become briefly unstable and thus rather than energy being transferred to the laser cooled ion for dissipation in the radiation field, kinetic energy is pumped into the ion from the electric fields if its $`q_{\mathrm{eff},u}`$ is greater than 0.908. We have shown that instability heating is minimised when $`q_u`$ is small. Since $`q_u`$ is a function of mass, trap frequency, and trap depth, although a lower mass limit exists below which sample ions will heat, this threshold can be changed by modifying the trap parameters. This mass limit corresponds to a second characteristic threshold value for $`q_u`$, above which heating will occur . More generally, there exists a three dimensional space spanned by mass, trap frequency, and trap amplitude which is divided into regions of sympathetic heating-cooling. Constraining any two of these variables will result in an upper or lower limit in the third.

## VII Acknowledgments

The Natural Sciences and Engineering Research Council of Canada supported this research. One of the authors (TH) would like to acknowledge the financial support of the Faculty of Graduate Studies in the University of Calgary. D. Hobill and H. Graumann provided very useful assistance related to the challenges of computational physics. The authors thank D. Feder for his helpful criticism and assistance in preparation of this manuscript.


text:
# Search for Long-lived States in Antiprotonic Lithium

## I Introduction

One of the most impressive success-stories of the last decade in few-body physics are the high-precision experimental and theoretical studies of long-lived states in antiprotonic helium exp ; kor ; elan ; kam . In view of this fact it is natural to pose the question, whether such long-lived antiprotonic states can exist in other systems, too. There were some experimental attempts to observe delayed components in annihilation products after stopping of slow antiprotons in different media wid , however, except for helium, no evidence of long-lived states was found. Theoretical predictions concerning the possible existence of such states could largely increase the willingness of experimentalists to go on with their attempts.

There are two possible directions in which further candidates for long-lived antiprotonic states could be searched for. First, one could consider atoms with electron structure similar to helium, that is, noble gases with closed outermost electron shells. This possibility was examined in some detail in Briggs  with an essentially negative answer concerning the possibility of formation of long-lived states after antiproton capture in noble gases.

The second possibility is to consider the next simplest atom, the lithium, which has three electrons, one of which could be replaced by the antiproton. In order to get an idea of the possibility of the occurrence of long-lived states in this four-body system, we have performed a semi-quantitative calculation of its level structure.

## II Calculation Method

The Hamiltonian of the $`(Li^{3+}+\overline{p}+2e)`$ four-body system can be written as

$$\widehat{H}=\frac{1}{2M}\mathrm{\Delta }_𝐑\frac{3}{R}+\widehat{h}^{(2)}(𝐫_1,𝐫_2;𝐑),$$
(1)
with the two-electron Hamiltonian

$$\widehat{h}^{(2)}(𝐫_1,𝐫_2;𝐑)=\widehat{h}^{(1)}(𝐫_1;𝐑)+\widehat{h}^{(1)}(𝐫_2;𝐑)+\frac{1}{|𝐫_1𝐫_2|}$$
(2)
The single-electron two-center Hamiltonian $`\widehat{h}^{(1)}`$ corresponds to the electron motion in the field of $`Li^{3+}`$ and $`\overline{p}`$ fixed at a distance $`𝐑`$:

$$\widehat{h}^{(1)}(𝐫_n;𝐑)=\frac{1}{2}\mathrm{\Delta }_{𝐫_i}\frac{3}{|𝐫_i\alpha 𝐑|}+\frac{1}{|𝐫_i+(1\alpha )𝐑|}$$
(3)
In eqs.(1-3) R is the vector pointing from $`Li^{3+}`$ to $`\overline{p}`$, while the $`𝐫_i`$ are the electron coordinates measured from the $`Li^{3+}\overline{p}`$ center of mass. 111 These are not exactly the Jacobian coordinates of the system, however, the extra terms in the Hamiltonian arising from the difference are very small and the accuracy of our calculation does not necessitate their consideration $`M`$ is the reduced mass of $`Li^{3+}`$ and $`\overline{p}`$:

$$\frac{1}{M}=\frac{1}{m_{Li^{3+}}}+\frac{1}{m_{\overline{p}}},$$

while $`\alpha `$ is defined as

$$\alpha =\frac{m_{\overline{p}}}{m_{Li^{3+}}+m_{\overline{p}}}$$

To calculate the spectrum of this four-body system or, at least, a part of it, we have used a Born-Oppenheimer-like (BO) approximation, in which the solution of the Schrödinger-equation

$$(\widehat{H}E_n^J)\mathrm{\Psi }_n^J(𝐫_1,𝐫_2,𝐑)=0$$
(4)
is attempted in two successive steps. First, the equation

$$(\widehat{h}^{(2)}(𝐫_1,𝐫_2;𝐑)\epsilon _{n\mu }^{(2)}(R))\mathrm{\Phi }_{n\mu }(𝐫_1,𝐫_2;𝐑)=0,$$
(5)
describing the motion of two electrons in the field of $`Li^{3+}`$ and $`\overline{p}`$ separated by a fixed vector R, has to be solved. The solutions $`\mathrm{\Phi }_{n\mu }`$ are characterized by the conserved quantum number $`\mu `$ – the sum of the electron angular momentum projections on the R direction. Next, the total wave function $`\mathrm{\Psi }_n^J(𝐫_1,𝐫_2,𝐑)`$ is approximated as a single product:

$$\mathrm{\Psi }_n^J(𝐫_1,𝐫_2,𝐑)=\frac{u_{n\mu }^J(R)}{R}D_{M_J,\mu }^J(\varphi ,\mathrm{\Theta },0)\mathrm{\Phi }_{n\mu }(𝐫_1,𝐫_2;𝐑),$$
(6)
where the Wigner’s $`D`$-functions are needed to ensure correct angular momentum quantum numbers $`JM_J`$ for the total wave function. The energy eigenvalues $`E_{n\mu }^J`$ are calculated from the radial equation

$$\left(\frac{1}{2M}\frac{d^2}{dR^2}+\frac{J(J+1)2\mu ^2}{2MR^2}\frac{3}{R}+\epsilon _{n\mu }^{(2)}(R)E_{n\mu }^J\right)u_{n\mu }^J(R)=0,$$
(7)
which is obtained by substituting Eq.(6) into Eq.(4), multiplying by $`D_{M_J,\mu }^J(\varphi ,\mathrm{\Theta },0)^{}\mathrm{\Phi }_{n\mu }(𝐫_1,𝐫_2;𝐑)^{}`$ and integrating over the electron coordinates and the angular variables $`(\varphi ,\mathrm{\Theta })`$ of R. In deriving Eq.(7), according to the adiabatic approximation, terms containing the derivatives of $`\mathrm{\Phi }_{n\mu }(𝐫_1,𝐫_2;𝐑)`$ with respect to R are neglected. In BO type calculations usually the lowermost electron configurations are used, for which $`\mu =0`$ (so called $`\sigma `$-term) in which case the $`D`$-function in Eq.(6) reduces to a spherical harmonics $`Y_{JM_J}`$.

For a one-electron problem Eq.(5) is replaced by

$$(\widehat{h}^{(1)}(𝐫_i;𝐑)\epsilon _{n\mu }^{(1)}(R))\phi _{n\mu }(𝐫_i;𝐑)=0$$
(8)
and it can be solved exactly by separation of variables in spheroidal coordinates $`\phi _{n\mu }(𝐫_i;𝐑)`$ being the well-known two-center wave functions. In our case, however, the problem of two interacting electrons in the field of two fixed Coulomb-centers is not solvable exactly and the determination of eigenvalues and eigenfunctions necessitates a highly non-trivial calculation. One could think of a variational approach of the type ahlr ; another possibility could be the diagonalization of the electron-electron interaction on the basis of two-center functions. In both cases the calculations are quite cumbersome and slowly converging when the size of the basis is increased. Since we believe, that the question of possible existence of metastability in antiprotonic lithium can be studied by a semi-qualitative exploration of its spectrum, we have chosen a simpler, elegant and efficient method of calculation for $`\epsilon _{n\mu }^{(2)}(R)`$ proposed in Briggs .

The main idea of the method is to represent the two-electron energy eigenvalue as a sum of two single-electron energies:

$$\epsilon ^{(2)}(R)=\epsilon _{Z_1,Z_2}^{(1)}(R)+\epsilon _{z_1,z_2}^{(1)}(R),$$
(9)
where $`\epsilon _{Z_1,Z_2}^{(1)}(R)`$ and $`\epsilon _{z_1,z_2}^{(1)}(R)`$ are energy eigenvalues of two-center equations of the type (8) with effective charges $`(Z_1,Z_2)`$ and $`(z_1,z_2)`$ instead of the physical charges $`(3,1)`$. The effective charges $`(z_1,z_2)`$ are chosen in such a way, that $`\epsilon _{z_1,z_2}^{(1)}(R=0)`$ and $`\epsilon _{z_1,z_2}^{(1)}(R=\mathrm{})`$ should reproduce the experimental values of the first ionization potentials of $`He`$ atom and $`Li^+`$ ion, respectively. Thus we get $`z_1=2.3578`$ and $`z_2=1.0135`$. As for $`(Z_1,Z_2)`$, the corresponding two-center eigenvalues in the $`R=0`$ and $`R=\mathrm{}`$ limits should reproduce the second ionization potentials of $`He`$ and $`Li^+`$, and the second electron is in this case the last one, therefore the physical values $`Z_1=3`$ and $`Z_2=1`$ were taken. The electron-electron repulsion is taken into account in this method by the deviation of the effective charges $`(z_1,z_2)`$ from their physical (integer) values. In this way the approximate $`\epsilon ^{(2)}(R)`$ of Eq.(9) reproduces the experimental two-electron binding for the two limiting cases $`R=0`$ and $`R=\mathrm{}`$, while for intermediate $`R`$-values the solution of the corresponding (effective) two-center problems seems to provide a reasonable interpolation prescription. This approach has been checked in the case of $`(He^{2+}+\overline{p}+2e)`$ and $`(Li^{3+}+\overline{p}+2e)`$ systems; its results were compared with those of a detailed variational calculation of ahlr . Results of the calculations are presented in Fig.(1). The agreement of $`\epsilon ^{(2)}(R)`$ obtained from Eq. (9) with the variational values is amazingly good in a wide range of $`R`$. The same procedure was applied to calculate the energy of the first excited electron configuration, where the limiting cases were adjusted to reproduce the energies of the first excited $`(1s2s)`$ states of $`He`$ and $`Li^+`$.

Having obtained the electronic energies $`\epsilon ^{(2)}(R)`$ the effective potentials of Eq.(7)

$$v_{\text{eff}}^{Jn}(R)=\frac{J(J+1)}{2MR^2}\frac{3}{R}+\epsilon _n^{(2)}(R),$$
(10)
can be calculated. Here the index $`\mu =0`$ was omitted, while the electron configuration label $`n`$ can take the values $`n=(1s1s),(1s2s)`$. For some values of the total angular momentum $`J`$ the effective potentials are shown in Fig.(2) for the ground- and first excited electron configuration.

The energy eigenvalues $`E_n^{J\nu }`$ are then calculated by solving Eq.(7) with these effective potentials, and the ”vibrational” quantum number $`\nu `$ is introduced to distinguish among the states with the same $`J`$-value.

## III Results and Discussion

The resulting spectrum of the $`(Li^{3+}+\overline{p}+2e)`$ system is shown in Fig.(3). Apart from the energy levels of the initial system (full black circles) we have shown also the energies of daughter states which can be formed after Auger-emission of one or two electrons. The open circles correspond to the energies of the $`(Li^{3+}+\overline{p}+e)`$ three-particle system, while the open squares are the hydrogen-like two-particle energies of the $`(Li^{3+}+\overline{p})`$ system.

Comparing the spectra of Fig.(3) with the well-known spectra of the $`He`$ atomcules (see e.g. shim ), we can find some apparent similarities and differences. The basic similarity can be formulated as follows: there are many states in the spectrum, from which Auger-emission is possible only with large electron orbital momentum and therefore is strongly suppressed. This could be one reason for metastability of these states; of course, this is only a necessary condition and by no means a sufficient one.

The basic difference, on the other hand, is the much higher density of states in the expected capture region (around the $`Li`$ atom ground state energy) which is due to the essential difference in the electron structure of $`He`$ and $`Li`$: the last electron is strongly bound in $`He`$, while very loosely in $`Li`$.

It can be noted, that the spectrum of the $`(Li^{3+}+\overline{p}+e)`$ system (open circles) strongly resembles the $`He`$ atomcule spectrum, therefore long-lived states in an isolated $`(Li^{3+}+\overline{p}+e)`$ system could be certainly expected. However, in contrast to the $`He`$ case, this system is charged and thus its interaction with atoms of the surrounding medium might be more violent, leading to a faster collisional de-excitation of these states.

## IV Conclusions

We have calculated the spectra of ($`Li^{3+}+\overline{p}+2e`$) and ($`Li^{3+}+\overline{p}+e`$) four- and three-body systems. Although the structure of the obtained spectra allows the existence of long-lived states, our calculations do not put us in a position to make definite statements about their experimental observability. This latter depends on several further factors, as well. One is the formation mechanism: for the time being we have no information about the population rate of the huge amount of states in the vicinity of expected antiproton capture energy. The physics of formation of the $`(Li^{3+}+\overline{p}+2e)`$ system in the reaction

$$(Li^{3+}+3e)+\overline{p}(Li^{3+}+\overline{p}+2e)+e$$

is quite different from the analogous process in $`He`$ due to the large difference in the binding energies of the outermost electron. In the $`Li`$ atom the first ionization potential is only $`0.198a.u.`$, which means that even adiabatic ionization is possible: when the distance between the antiproton and $`Li`$ atom becomes less than $`5a.u.`$, the binding energy of the last electron in their common field becomes zero and the electron is ”pushed” into a continuum state.

Another unknown factor is the way, how the eventually formed $`(Li^{3+}+\overline{p}+2e)`$ systems interact with the media atoms and what is the role of collisional de-excitation in their life-time. In any case, in order to reduce the undesired effect of this factor, probably, the experiments looking for long-lived states should be performed in dilute vapors of $`Li`$.

###### Acknowledgements.

One of the authors (VB) wishes to thank the NATO senior fellowship grant 1004/NATO/01, which made possible his visit to Hungary, while the other author (JR) is grateful for the OTKA grants T 026244 and T 029440.