text:
# From the triangular to the kagome lattice: Following the footprints of the ordered state.

## I Introduction

Geometrically frustrated antiferromagnets are the historical candidates for the realization of a spin liquid ground state. Indeed, the spin-half Heisenberg antiferromagnet on the triangular lattice was the first model to be proposed by Anderson and Fazekas an ; fa-an  in 1973 as a system where geometric frustration and quantum fluctuations could prevent zero-temperature magnetic ordering in two dimensions, stabilizing instead a ground state with gapped spin excitations and exponentially decaying correlations. Since then, a good amount of work has been devoted to investigate the nature of the ground state of the triangular Heisenberg model huse ; tri1 ; bernu1 ; bernu2 ; tri2 ; tri3 ; tri4 ; tri5 ; aza ; luc ; luct ; boni  which remained an open question until quite recently. At present, however, there is a general consensus on the existence of long-range antiferromagnetic order following a $`120^{}`$ Nèel pattern in the ground state of this modelbernu1 ; bernu2 ; luc ; boni : frustration and quantum fluctuations on the two-dimensional triangular lattice are not strong enough to stabilize a non-magnetic ground state.

A more promising candidate for a disordered ground state can be obtained through a ‘dilution’ of the triangular lattice, leading to the so-called kagome net (Fig. 1). In fact, on this geometry, due to the lower coordination ($`z=4`$ compared with $`z=6`$ in the triangular case), frustration is much stronger and even in the classical limit it gives rise to an infinite number of classical ground states, with ordered and disordered configurations degenerate in energy.kagoc1 ; kagoc2 ; kagoc3  Due to the extensive entropy of the classical ground state, the so-called order from disorder mechanism – usually stabilizing, among degenerate manifolds, long-range ordered configurations – is much less effective than in other frustrated models. In particular, while harmonic fluctuations select planar configurations, they turn out to be completely insensitive to their degree of order, and only non-linear effects eventually stabilize a classical ground state with Nèel correlations with a $`\sqrt{3}\times \sqrt{3}`$ pattern (see Fig. 1). Whether such classical minimum energy configuration possesses true long-range order kagoc1  or it is a critical point with power-law correlations kagoc2 ; kagoc3  is still an open issue. In any case, due to the particularly delicate mechanism leading to the classical antiferromagnetic order, the latter is expected to be easily destabilized by quantum fluctuations.

The investigation of the quantum Heisenberg antiferromagnet on the kagome lattice has been capturing increasing attention kago1 ; kago2 ; kago3 ; kago4 ; kago5 ; kago6 ; kago7 ; lech ; wal ; mam ; t2 ; lhul  for some time now. Mainly on the basis of numerical work, evidence has accumulated supporting a spin-liquid ground state, even though a very peculiar one: this, in fact, would be characterized by a small gap ($`J/20`$) to spin excitations, and by an exponentially large number of singlets contiguous to the ground state.lech ; wal ; mam ; t2 ; lhul  The classification of spin liquid of Type II has been recently proposed t2  to classify this particular behavior.

On the experimental side, the triangular geometry provides the scenario for interesting physical phenomena taking place in the superconducting compounds NaxCoO2-yH2O trisu  and the organic materials $`\kappa `$-(BEDT-TTF)2X, with X being I3, Cu\[N(CN)2\]Br or Cu(SCN)2 orga . Two spin-1/2 kagome-like materials have been also recently reported, the Volborthite Cu3V2O7(OH)2 2H2O vol  and the kagome-staircase compounds Ni3V2O8 and Co3V2O8 kast .

In the present work, we consider the spin-$`1/2`$ Heisenberg antiferromagnet on a lattice that interpolates between the triangular and the kagome ones. The Hamiltonian is

$$\widehat{H}=J\underset{ij}{}\widehat{𝐒}_i\widehat{𝐒}_j+J^{}\underset{ij^{}}{}\widehat{𝐒}_i\widehat{𝐒}_j^{},$$
(1)
where $`\widehat{𝐒}_i`$ are spin-half operators, $`ij`$ denotes nearest-neighbor bonds belonging to the kagome lattice, and $`ij^{}`$ are the remaining bonds connecting kagome and non-kagome sites. A scheme is indicated in Fig. 1. In this way, $`J^{}=J`$ corresponds to the usual triangular lattice, while $`J^{}=0`$ defines the kagome one. Our aim is to start from the triangular limit and to investigate the stability of the ordered state as $`J^{}/J`$ decreases. This is sensible because the classical Nèel state on the triangular lattice is compatible with the expected $`\sqrt{3}\times \sqrt{3}`$ classical ordering on the kagome antiferromagnet. The first to study this model were Zeng and Elser kago5 , who performed a spin-wave analysis and concluded that, for spin-1/2 particles, the ordered state could be stable for $`J^{}`$ down to $`J^{}/J0.2`$. Very recently, this model has been investigated with a coupled cluster treatment farn . This technique is based on the three-sublattice structure characterizing the $`120^o`$-Néel order of the triangular lattice, which is found to break down very close to $`J^{}=0`$, indicating the instability of the magnetic ordered state very close to the kagome limit.

In this paper, we tackle this problem using variational approaches, and exact diagonalization of small clusters. In particular, in Section II, we employ the so-called fixed-node (FN) technique fn  in order to improve the accuracy of a wave function with long-range antiferromagnetic order previously introduced in the pure triangular case huse ; luc . This technique, is usually used in the context of quantum Monte Carlo simulations as a method to approximate the Hamiltonians affected by sign-problem instabilities fn  and obtain exact ground-state properties of the corresponding “effective Hamiltonian”, no longer affected by the sign problem. Here we use the FN method to define a variational state with long-range antiferromagnetic order, and we check its accuracy in describing the ground state of the model as $`J^{}/J`$ is reduced through a direct comparison with exact diagonalization results on the $`6\times 6`$ cluster. In order to detect any indication of a change in the nature of the ground state approaching the kagome limit, in Section III we analyze the structure of the low-energy spectra as a function of $`J^{}/J`$, using the exact diagonalization of the Hamiltonian on small clusters. Section IV is finally devoted to summary and conclusions.

## II Variational approaches

A fairly accurate representation of the ground state of the spin-$`1/2`$ Heisenberg antiferromagnet on the triangular lattice can be obtained starting from a $`120^{}`$ Néel ordered state and including Gaussian fluctuations by means of a Jastrow factor containing two-spin correlationshuse ; luc

$$|\psi _v=\widehat{P}_0\mathrm{exp}\left(\frac{1}{2}\underset{i,j}{}v(ij)\widehat{S}_i^z\widehat{S}_j^z\right)|N,$$
(2)
where $`\widehat{P}_0`$ is the projector onto the $`S^z=0`$ subspace, $`|N`$ is the classical Néel state in the $`xy`$ plane,

$$|N=\underset{x}{}\mathrm{exp}\left[\frac{2\pi i}{3}\left(\underset{i\mathrm{B}}{}S_i^z\underset{i\mathrm{C}}{}S_i^z\right)\right]|x,$$
(3)
and $`|x`$ is an Ising spin configuration specified by assigning the value of $`S_i^z`$ for each site. On the square lattice case, the classical Néel state reproduces exactly the phases of the ground state of the Heisenberg Hamiltonian according to the Marshall theorem.luc  On the triangular lattice, instead, the exact phases of the ground state are unknown, and the classical part of the wave function (2) does not reproduces them accurately. However, as originally suggested by Huse and Elser, huse  a very accurate ansatz of the ground-state phases can be obtained by including three-spin correlation factors of the form:

$$T(x)=\mathrm{exp}\left(i\beta \underset{i,j,k}{}\gamma _{ijk}S_i^zS_j^zS_k^z\right),$$
(4)
defined by the coefficients $`\gamma _{ijk}=0,\pm 1`$, appropriately chosen so as to preserve the symmetries of the classical Néel state, and by the variational parameter $`\beta `$. In particular the sum in Eq. (4) runs over all distinct triplets of sites $`i,j,k`$ where both $`i`$ and $`k`$ are nearest neighbors of $`j`$, and $`i`$ and $`k`$ are next-nearest neighbors to one another. The sign factor $`\gamma _{ijk}=\gamma _{kji}=\pm 1`$ is invariant under rigid translations and rotations in real space by an angle of $`120^{}`$ of the three-spin cluster $`i,j,k`$, but changes sign under rotations by $`60^{}`$. The resulting wave function reads therefore:

$$|\psi _v=\widehat{P}_0\underset{x}{}\mathrm{\Omega }(x)\mathrm{exp}\left(\frac{1}{2}\underset{i,j}{}v(ij)S_i^zS_j^z\right)|x,$$
(5)
with the phase factor given by

$$\mathrm{\Omega }(x)=T(x)\mathrm{exp}\left[\frac{2\pi i}{3}\left(\underset{i\mathrm{B}}{}S_i^z\underset{i\mathrm{C}}{}S_i^z\right)\right].$$
(6)
Since the Hamiltonian is real, a better variational wave function is defined by the real part of Eq.5.

The two-body Jastrow potential in (5) contains in principle as many variational parameters as the independent distances on the lattice. However, the same level of accuracy can be obtained by optimizing separately the nearest-neighbor and next nearest-neighbor distances and adopting for the longer-range correlations an expression based on the consistency with linear spin-wave theory: luc

$$v(r)=\frac{\eta _{\mathrm{}}}{N}\underset{q0}{}e^{iqr}v_q$$
(7)
with

$$v_q=1\sqrt{\frac{1+2\gamma _k}{1\gamma _k}},$$
(8)
$`\gamma _k=2(\mathrm{cos}k_x+2\mathrm{cos}k_x/2\mathrm{cos}\sqrt{3}k_y/2)`$ and $`\eta _{\mathrm{}}`$ is a variational parameter. For the anisotropic triangular lattice (see Fig.1) we have optimized separately the nearest-neighbor bonds connecting two kagome sites ($`\eta _1`$), and a kagome site and a non-kagome site ($`\eta _1^{}`$), as well as the next-nearest-neighbor bonds ($`\eta _2`$). As a result, the total number of variational parameters for the present variational wave function is 5. Their optimal values, and the corresponding variational energies, are reported for $`N=36`$ in Tab. 1, for various values of $`J^{}/J`$.

The accuracy of the present long-range ordered wave function in the triangular lattice limit ($`J^{}/J=1`$) has been analyzed in detail in Ref.luc ; luct . In this limit the wave function is known to provide a qualitatively correct representation of the ground-state correlations.huse ; luc ; boni ; luct  In order to check the accuracy within a larger range of $`J^{}/J`$ we have compared several variational properties with the exact ground-state values calculated by exact diagonalization on the largest cluster presently accessible, $`N=36`$. As shown in Fig. 2, for $`0.4J^{}/J1`$ both the accuracy on the ground-state energy and the overlap with the exact ground state, remain approximatively constant and equal to the values in the triangular limit ($`J^{}/J=1`$). In addition, the wave function (2) provides an accurate representation of the phases of the actual ground state up to $`J^{}/J0.2`$, as it can be checked by measuring the average-sign

$$s=\underset{x}{}\mathrm{sgn}[\psi (x)\psi _0(x)]|\psi _0(x)|^2,$$
(9)
with $`|\psi _0=_x\psi (x)|x`$ (Fig. 2). This remarkable feature is due to the presence of the triplet term (4), allowing to adjust the phases in a non-trivial way, without changing the underlying Néel order. For instance, in the triangular limit the average sign (overlap) of the wave function is 0.733 (0.562) without the triplet term and 0.932 (0.779) with it.

Since the variational wave function (2) reproduces accurately the phases of the ground state its quality can be improved by adopting the FN scheme of Ref.fn . This allows one to obtain a new variational wave function, $`|\psi _{\mathrm{fn}}`$, defined as the ground state of the so-called FN effective Hamiltonian, whose matrix elements, $`H_{x,x^{}}^{\mathrm{eff}}`$, can be constructed starting from the original Hamiltonian and a variational guess on the ground-state phases given, in our case, by the wave function (5), $`|\psi _v=_x\psi _v(x)|x`$:

$`H_{x^{},x}^{\mathrm{eff}}=\{\begin{array}{ccc}H_{x^{},x}& \mathrm{if}& \overline{H}_{x^{},x}0\hfill \\ 0& \mathrm{if}& \overline{H}_{x^{},x}>0\hfill \\ H_{x,x}+𝒱(x)& x=x^{}& \end{array}`$ (13)
where $`\overline{H}_{x^{},x}=\psi _v(x^{})H_{x^{},x}/\psi _v(x)`$, and

$$𝒱(x)=\underset{\{\overline{H}_{x^{},x}>0,x^{}x\}}{}\overline{H}_{x^{},x}.$$
(14)
Indicating with $`E_0^{\mathrm{fn}}`$ the ground-state energy of the FN Hamiltonian, with $`E_v^{\mathrm{fn}}=\psi _{\mathrm{fn}}|\widehat{H}|\psi _{\mathrm{fn}}/\psi _{\mathrm{fn}}|\psi _{\mathrm{fn}}`$ and $`E_v=\psi _v|\widehat{H}|\psi _v/\psi _v|\psi _v`$ the energy expectation values on $`|\psi _{\mathrm{fn}}`$ and $`|\psi _v`$, respectively, it is possible to show fn  that the following chain of inequalities holds:

$$E_vE_0^{\mathrm{fn}}E_v^{\mathrm{fn}}E_0,$$

where $`E_0`$ is the ground-state energy of $`\widehat{H}`$. Hence, the FN procedure is granted to produce a wave function with a better variational energy than $`|\psi _v`$. In addition, also the lowest eigenvalue of the FN Hamiltonian, $`E_0^{\mathrm{fn}}`$, gives an upper bound of the ground-state energy better than the variational energy $`E_v`$. This is the quantity which is usually considered in the quantum Monte Carlo application since it is the most directly accessible.

The FN Hamiltonian is explicitly defined in such a way that the matrix $`\overline{H}_{x^{},x}^{\mathrm{eff}}=\mathrm{sgn}[\psi _v(x^{})]H_{x^{},x}^{\mathrm{eff}}\mathrm{sgn}[\psi _v(x)]`$ has all negative off-diagonal matrix elements. By the Perron-Frobenius theorem, the amplitudes of the ground state of this matrix, $`\overline{\psi }_{\mathrm{fn}}(x)`$, have the same sign for all the configurations $`|x`$. This implies in turn that the FN ground state, $`|\psi _{\mathrm{fn}}`$, and the starting variational wave function, $`|\psi _v`$, have exactly the same phases. In fact, the amplitudes of the ground state of $`H_{x^{},x}^{\mathrm{eff}}`$, $`\psi _{\mathrm{fn}}(x)`$, are related to the amplitudes of the ground state of the transformed matrix $`\overline{H}_{x^{},x}^{\mathrm{eff}}`$, $`\overline{\psi }_{\mathrm{fn}}(x)`$, by the simple relation $`\overline{\psi }_{\mathrm{fn}}(x)=\psi _{\mathrm{fn}}(x)\mathrm{sgn}[\psi _v(x)]`$. For this reason, the FN wave function is expected to provide an accurate description of the ground state only when it is constructed with a good variational ansatz of the ground-state phases. In the present case for $`0.2J^{}/J1`$.

In order to thoroughly check the accuracy of the FN wave function, we have exactly diagonalized with the Lanczos algorithm the FN Hamiltonian for several values of $`J^{}/J`$ on the $`N=36`$ cluster. This has allowed us to calculate not only the FN energies but also the overlap of the FN wave function with the exact ground state. As shown in Fig. 2, both the FN upper bounds to the ground-state energy, $`E_v^{\mathrm{fn}}`$ and $`E_0^{\mathrm{fn}}`$, are sizably more accurate than the simple variational estimate $`E_v`$. In particular, the FN wave functions has a much higher overlap than $`|\psi _v`$, and its accuracy is almost constant and comparable to the one in the triangular limit down to values of $`J^{}/J`$ as small as $`0.2`$.

We have finally compared the variational, FN, and exact estimates of the antiferromagnetic order parameter for a $`120^{}`$ Nèel order,

$$m^2=36\frac{^2}{N(N+6)},$$
(15)
where $`^2`$ is the sublattice magnetization squared bernu2 . Interestingly, by decreasing $`J^{}/J`$ the exact ground-state order parameter remains approximatively constant down to $`J^{}/J0.2`$ thus indicating a possible destabilization of the antiferromagnetic order only very close to the kagome limit. In addition, though the variational and FN estimates of the order parameter are approximatively $`10\%`$ higher than the exact one the same degree of agreement is observed in all the range $`0.2J^{}/J1`$. For the $`N=36`$ cluster, the FN wave function, based on the variational wave function (5), provides a good quantitative description of the exact ground state for $`0.2J^{}/J1`$. In this range, the variational and the exact expectation values of antiferromagnetic order parameter remain constant and equal to their values in the triangular limit. Below $`J^{}/J0.2`$, the exact value of the order parameter begins to decrease and the accuracy of the our Nèel ordered wave function quickly degrades, indicating a change in the ground-state correlations only very close to the kagome limit.

In order to support further the stability of the antiferromagnetic phase for very small values of $`J^{}/J`$ we have extended the variational and the FN calculations of the order parameter $`m^{}`$ to much larger sizes by using quantum Monte Carlo techniques.calandra  As shown in Fig. 3, even at a value of $`J^{}/J=0.2`$, the lowest coupling ratio when the variational wave function is expected to be accurate (see Fig.2), the order parameter $`m^{}`$ remains sizably larger than the triangular $`J^{}/J=1`$ case. This indicates that, as we increase the size at fixed ratio $`J^{}/J<1`$, the stability of the ordered phase, already evident in the exact diagonalization in the $`N=36`$ cluster (see Fig.2), becomes more and more clear. We expect that this qualitative behavior, confirmed both by the exact diagonalization, and even more strongly by the variational and the FN approaches on larger sizes, is a genuine feature of the model, even though the quantitative results that we have obtained by quantum Monte Carlo may be affected by a sizable error. This feature may appear rather surprising, as the quantum fluctuations should increase for small $`J^{}/J`$ and should tend to destabilize the ordered phase, as expected for instance within spin-wave theory.kago1  However, the wave function that we have used is consistent with spin wave theory in the large spin limit, and therefore, since also at the variational level the value of $`m^{}`$ increases, we conclude that the quantum fluctuations are not very accurately described within a method that is not controlled by the variational principle (the large spin limit), at least in the region of small $`J^{}/J`$.

A similar size-scaling analysis can be carried out to estimate the ground-state energy per spin in the thermodynamic limit as illustrated in Fig. 4. The extrapolated energies are listed in Tab. 2. Since the FN energy error is known exactly up to the $`6\times 6`$ cluster, we can estimate the ground-state energy in the thermodynamic limit, by adding the $`6\times 6`$ correction to the infinite-size FN estimate. The corresponding extrapolated values of the ground-state energy shown in the Tab.2 represent reasonable benchmark values for this quantity (or at least good upper bounds), as the FN error is expected to increase weakly for larger sizes even for a good variational ansatz. This behavior has been verified up to $`6\times 6`$, and is also very reasonable to expect in general for an approximate variational calculation.

In the triangular limit, $`J^{}/J=1`$, our estimated ground-state energy coincides with the corresponding one obtained with the coupled cluster method farn  by extrapolating in the size of the clusters ($`m\mathrm{}`$, with the notations of Ref. farn ). This value is instead slightly lower than the extrapolated results based on small clusters ($`N36`$) obtained in Ref.bernu2 , giving $`e_0=0.5445`$. A good agreement between the coupled cluster for $`m=6`$ and the FN method is also seen for values of $`J^{}/J`$ down to $`0.2`$ (Tab. 2). Furthermore, our extrapolated energies, based on the error for the $`6`$ cluster and the FN energy, remain lower than the $`m=6`$ coupled cluster result by a similar amount, suggesting that our variational approach remains enough accurate also in this region (unfortunately the extrapolations $`m\mathrm{}`$ are not given in Ref.farn  for $`0<J^{}/J<1`$). Instead in the kagome limit our variational ansatz should be clearly less accurate, as the ground state is believed to be a spin liquid with no magnetic order in the thermodynamic limit, i.e. qualitatively different from our initial variational guess given by Eq. (2). This may explain why in this case our energy estimate is slightly higher than the $`m\mathrm{}`$ coupled cluster result. We note, however, that our values reported in Tab. (2) represent in general reasonable upper bounds for the energy, as they are obtained by a rigorous variational method such as the FN one.

## III Low-energy excitations

An effective method to investigate the possibility of magnetic order in spin systems is to analyze the structure of the spectrum of finite-size samples, following the strategy of Refs. bernu1 ; bernu2 ; lech . The tendency toward antiferromagnetic order in the thermodynamic limit manifests itself in finite-size clusters through the fact that the spin excitations with the lowest energies can be described by the effective Hamiltonian of a “quantum top”. Within such a description, the energy of the lowest levels in the different subspaces labeled by total spin $`S`$, can be approximated by,

$$E_N(S)E_N(0)=\frac{S(S+1)}{2I_N},$$
(16)
where $`I_N`$ is the inertia of the top, which is an extensive quantity. Hence, the plots of the lowest energy levels as functions of $`S(S+1)`$, have the appearance of a “Pisa tower” with a slope that decreases as $`N`$ increases. An important property of the states of the “Pisa tower” (indicated in Ref. bernu2  quasi degenerate joint states) is that they belong to irreducible representations of the point group that are compatible with the symmetry of the ordered state. In the case of the $`\sqrt{3}\times \sqrt{3}`$ order, these are the representations labeled as $`\mathrm{\Gamma }_1,\mathrm{\Gamma }_2,\mathrm{\Gamma }_3`$ of the $`C_{3v}`$ group, which correspond to $`[𝐤=0,_\pi \mathrm{\Psi }=\mathrm{\Psi },_{2\pi /3}\mathrm{\Psi }=\mathrm{\Psi },\sigma _x\mathrm{\Psi }=\mathrm{\Psi }];[𝐤=0,_\pi \mathrm{\Psi }=\mathrm{\Psi },_{2\pi /3}\mathrm{\Psi }=\mathrm{\Psi },\sigma _x\mathrm{\Psi }=\mathrm{\Psi }];[𝐤=𝐐,_\pi \mathrm{\Psi }=\mathrm{\Psi },_{2\pi /3}\mathrm{\Psi }=\mathrm{\Psi },\sigma _x\mathrm{\Psi }=\mathrm{\Psi }]`$, respectively. We have followed the same notation of Refs. bernu1 ; bernu2 ; lech , $`_\varphi `$ denoting a rotation of $`\varphi `$, $`\sigma _x`$, a reflection with respect to a mirror plane with the normal pointing along $`x`$, and being $`𝐐=(2\pi /3,2\pi /3)`$, the corner of the Brillouin zone of the kagome lattice.

There are not so many finite periodic clusters, which can be exactly diagonalized and that interpolate between the triangular and kagome lattices without frustrating the antiferromagnetic order. The smallest ones are $`N=12`$ and $`N=36`$, which evolve, respectively, toward $`N=9`$ and $`N=27`$ at the pure kagome limit ($`J^{}=0`$). The unit cell of the interpolating lattice contains $`2\times 2`$ unit cells of the pure triangular lattice, implying a reduction of the translation operations of the periodic cluster of a factor $`4`$ with respect to the number of translations at the pure triangular limit. This implies a much bigger Hilbert space to be dealt with in the numerical procedure. In concreteness, even exploiting all the available symmetries, the number of states is $`42035724`$ for $`N=36`$ spins. In this cluster, a complete study of several excited states within different spin sectors for several values of $`J^{}`$ is prohibitive from the computational point of view, while the $`N=12`$ cluster maybe too small. We have, therefore, included in the analysis the clusters with $`N=16`$ and $`N=28`$ sites (which evolve toward kagome clusters with $`N=12`$ and $`N=21`$, respectively), by introducing twisted boundary conditions as explained in Refs. bernu2 ; lech . The latter are equivalent to rotate the local frame in $`(\pm 2\pi /3,2\pi /3)`$ at each translation along a unitary lattice vector. This procedure restores the otherwise frustrated antiferromagnetic order in the $`x,y`$ spin plane but obviously breaks the spin-rotational symmetry as well as some symmetry operations of the point group. For this reason, in this case, a law like that expressed in (16), with the replacement $`S(S+1)S_z^2`$, should be obeyed by the quasi degenerate joint states, as the total $`S^z`$ along the z axis, cannot couple with the remaining total spin components in the clusters with twisted boundary conditions. Results are shown in Figs. 56, and 7 for $`N=12,16`$, and $`28`$, respectively. Figures 5 and 6 show the full spectra, obtained by diagonalizing all the blocks of the Hamiltonian matrix. In the case of $`N=28`$ (Fig. 7), we have not diagonalized the full Hamiltonian but obtained the lowest eigenvalues within all the subspaces using the Lanczos algorithm.

We have plotted the states with symmetries compatible (non-compatible) with the antiferromagnetic $`\sqrt{3}\times \sqrt{3}`$ order with crosses (circles). The first striking feature is the fact that only a subset of the states forming the “Pisa tower” in the triangular limit remains aligned as $`J^{}/J`$ decrease. Those states build up, so to say, a “small Pisa tower”, whose slope decreases with $`J^{}`$.

In order to analyze in more detail the behavior of the states along the “small Pisa tower” we show zooms of the low-energy and low-spin sector of the spectra for the $`N=12`$, and $`N=28`$ clusters in Figs. 8 and 9, respectively. In the case of $`N=12`$, shown in Fig. 8, this small set contains levels that belong to representations compatible with the $`\sqrt{3}\times \sqrt{3}`$ order and they remain aligned within all the range $`0.1J^{}/J1`$. Signatures of departure from that behavior are observed for very low $`J^{}/J=0.05`$, where the lowest level with $`S=2`$ deviates from the line of the “small Pisa tower” while a state with a symmetry not compatible with the antiferromagnetic order becomes quasi degenerate with it. In Fig. 8 the data for the related $`N=9`$ kagome lattice are also shown. In that case, a similar deviation of the state with $`S=5/2`$ from the line connecting the lowest levels in the $`S=1/2`$ and $`S=3/2`$ sectors is observed but within a scale which is an order of magnitude larger in the energy axis. Another contrasting feature that comes from the comparison between $`J^{}=0`$ and $`J^{}0`$ is that in the pure kagome limit, there are two states with $`S=1/2`$ (one of them belonging to the manifold of the degenerate ground state) whose symmetries are not compatible with the $`\sqrt{3}\times \sqrt{3}`$ order that have energies within the spin gap. Instead for $`J^{}/J=0.05`$ the spin gap is clean from such states. In addition, the gap to the lowest excitation with $`\mathrm{\Delta }S=1`$ does not evolve with continuity for $`J^{}/J0`$ to the corresponding gap in the kagome limit (see also Fig. 10).

Similar remarks apply to the 16- and 28- sites clusters shown in Figs. 6 and 7. In these cases, see e.g., Fig. 9, the “small Pisa Tower” contains a larger number of states but the deviation from a perfect alignment is observed already at $`J^{}/J0.2`$, where states within subspaces not compatible with the $`\sqrt{3}\times \sqrt{3}`$ order appear at low energies. In particular, for $`S_z=1`$ one of those states is almost degenerate with the one belonging to the Pisa tower while for $`S_z=2`$ it is well below it. As in the $`N=12`$ cluster, also for $`N=16`$ and 28 the gap to the first spin excitation monotonically decreases with $`J^{}`$ (Fig. 10).

To complete the analysis, we have computed the ground-state energy ($`E_0`$) and the lowest eigenvalue in the subspace with $`S=1`$ within the representation $`\mathrm{\Gamma }_1`$ of $`C_{3v}`$ ($`E_{S=1,\mathrm{\Gamma }_1}`$), in the periodic lattice with $`N=36`$ sites. The latter does not actually correspond to the lowest energy $`S=1`$ excitation in the pure triangular limit, which is in the subspace corresponding to $`\mathrm{\Gamma }_2`$. In any case, the energy difference $`\mathrm{\Delta }=E_{S=1,\mathrm{\Gamma }_1}E_0`$ is, in general, an upper bound for the spin gap of the cluster in all the range $`0<J^{}/J1`$. The behavior of $`\mathrm{\Delta }`$ as a function of $`J^{}`$ is shown in Fig. 10. As in the clusters previously considered, the gap decreases rapidly as $`J^{}`$ decreases, and, in particular, it does not evolve continuously for $`J^{}0`$ towards the value of the corresponding kagome lattice ($`N=27`$, indicated with an open square in Fig. 10). Instead, it tends to a value which is even smaller than the magnitude of the spin gap for the $`N=36`$ kagome lattice (indicated by a cross in the same figure). This is in contrast with the behavior exhibited by the ground-state energy, which evolves smoothly toward the value of the $`N=27`$ kagome lattice (see the inset of Fig. 10).

The behavior of the spin gap and the spectra strongly suggests the closing of the spin-gap in the thermodynamic limit within the whole range of $`0<J^{}/J1`$. In fact, such gap is known to close in the thermodynamic limit at the triangular point ($`J^{}/J=1`$) bernu1 ; bernu2 ; luc  and it decreases systematically by decreasing $`J^{}`$ in all the cluster considered as the low-energy scale for the lowest-spin excitations is set by $`J^{}`$, i.e., by the states of the “small Pisa tower”. This can be understood in a simple uncorrelated framework where for $`J^{}<J`$ the lowest spin-excitations are clearly obtained through a spin-flip on a non-kagome site, with an energy cost $`6J^{}`$, compared to a cost $`4J+2J^{}`$ for a spin flip on a kagome site. As the number of non-kagome spins is just a fraction of the total number of sites, on a finite cluster such a mechanism would apply only to the lowest spin excitations: hence the reduction of the number of states belonging to the “Pisa tower” with respect to the triangular case.

This simple picture immediately suggests that the nature of the spin excitations is intrinsically different for finite $`J^{}`$ and in the pure kagome limit. In the latter case, in fact, the spins on the non-kagome sites do not belong to the Hilbert space of the model and such low-energy spin excitations are not possible. For this reason by turning on the $`J^{}`$, the spin gap has a finite discontinuity and it is reasonable to expect the model to be unstable under the perturbation introduced by $`J^{}`$ bonds. This is also confirmed by the analysis of the low-energy excitations for small values of $`J^{}/J`$ (see Fig. 8)

One step further in the analysis of the nature of the ground state leads to the question whether or not the latter remains ordered down to the kagome limit. The structure of the spectra of the small $`N=12`$ cluster indicates that could well be the case, while those of $`N=16`$ and $`N=28`$ suggests some change in the nature of the ground state around $`J^{}/J0.2`$. It is worth recalling, however, that the results for the two latter clusters have been obtained by using twisted boundary conditions which could lead to weaker signatures of ordering. On the other hand, in two dimensions, it seems difficult to conceal the possibility of the closing of the spin gap with the absence of some kind of magnetic order. A simple geometrical analysis reveals that the non-kagome spins, which we have already identified as the responsible of the lowest-energy excitations, form themselves a triangular lattice with a cell parameter twice the size of the one of the usual triangular lattice. Therefore, a possible scenario for the evolution of the ground state as $`J^{}`$ decreases could be that at some point a crossover takes place from a magnetic order in the usual $`\sqrt{3}\times \sqrt{3}`$ pattern to a magnetic order with a similar pattern but in the triangular lattice of the non-kagome spins. Both kinds of order are commensurate and the ground state could undergo a smooth evolution from one to the other.

## IV Summary and Conclusions

We have investigated the low energy properties of a $`J^{}J`$ triangular lattice that interpolates between the usual triangular ($`J^{}/J=1`$) and kagome ($`J^{}=0`$) lattices. To this end, we have used a variational approach based on a FN wave function accurately describing the ground state in the triangular ($`J^{}/J=1`$) limit, and exact diagonalization techniques.

We have analyzed the quality of the approximation to the exact ground state, provided by the FN technique in a periodic cluster with $`N=36`$ sites and then extended the calculation up to $`N=144`$ sites by using quantum Monte Carlo. We have found that such a wave function, describing a state with a $`\sqrt{3}\times \sqrt{3}`$ Nèel ordered phase, is very close to the exact one for $`0.2J^{}/J1`$ and that in this range its accuracy is almost independent of the $`J^{}/J`$ ratio. Consistently, the ground-state expectation value of the antiferromagnetic order parameter remains approximately constant down to $`J^{}/J0.2`$. Below this value the order parameter is suppressed and the quality of the FN wave function degrades, suggesting a change in the nature of the ground-state. This is also confirmed by the analysis of the low-energy spectra on small clusters, showing some signatures of instability of Nèel ordering for $`J^{}/J0.2`$. However, at the same time, such analysis also indicates quite clearly that the low-energy scale for spin excitations is set by $`J^{}`$ and that in particular the spin-gap has a finite discontinuity at $`J^{}/J=0`$, the triplet excitations being gapless for any non zero value of $`J^{}/J`$. This would imply that the kagome disordered phase is unstable against a slight perturbation tending to restore the $`z=6`$ coordination number of the triangular lattice. A possible scenario is that within the large amount of singlets that are quasi degenerate with the ground state in the kagome clusters, the one corresponding to the $`\sqrt{3}\times \sqrt{3}`$ state is favored by some kick produced by $`J^{}`$.

Our results are in disagreement with the predictions of the spin-wave theory of Ref. kago1 , indicating a progressive reduction of the antiferromagnetic ordering for $`J^{}/J<1`$ and a complete melting of the $`\sqrt{3}\times \sqrt{3}`$ order for $`(J^{}/J)_c0.2`$. Instead, our conclusions are closer to those of the recent coupled cluster treatment of Farnell et al. farn  providing evidence for instability of the antiferromagnetic order very close or possibly at the kagome point ($`(J^{}/J)_c=0.0\pm 0.1`$), and the existence for small $`J^{}/J`$ of a regime whose correlations are very different to those of the triangular antiferromagnet. Further investigations are necessary to clarify the nature of this regime.

## V Acknowledgments

L. A. thanks Prof. Fulde for his hospitality as well as the the Alexander-von-Humboldt Stiftung and CONICET, Argentina for the support. This work was partially supported by INFM-PRA MALODI, and by NSF under Grant No. DMR02-11166 (L.C.).


format: text

text:
# Radiative collisional heating at the Doppler limit for laser-cooled magnesium atoms

## I Introduction

Laser cooling and trapping methods are an important ingredient in the recent achievements in the low-temperature physics of gaseous atoms and molecules Metcalf03 . In alkali atoms one can reach very low temperatures with Sisyphus and polarization gradient techniques, which surpass the Doppler cooling method in efficiency, and appear in magneto-optical traps without additional efforts Lett89 . The drawback is that it becomes impossible to test the basic two-state Doppler cooling theory with these systems.

The appearance of sub-Doppler cooling is based on the hyperfine structure of the alkali atoms. The same hyperfine structure is reflected in the complicated molecular state structure of the quasimolecule formed by two colliding atoms. In dilute atomic gases binary collisions dominate the atomic interactions, and the quasimolecule states couple with the cooling and trapping light field. This leads to inelastic light-assisted collisions that cause either loss or heating of the trapped atoms Suominen96 ; Suominen98 ; Weiner99 . Thus in alkali atoms the hyperfine structure leads to difficulties in testing the Doppler theory, as well as in modelling the collisional processes.

The major isotopes of alkaline earth atoms, on the other hand, have no nuclear spin and thus the hyperfine structure is missing. This makes the Doppler limit for temperatures the true laser cooling limit, and allows for testing the basic Doppler theory. Experimental studies have been published recently for 88Sr Xu02 , and have been obtained for 24Mg as well Thomsen03 . They indicate an intensity-dependent heating rate which prevents reaching the theoretical Doppler limit.

Inelastic laser-assisted collisions are expected to depend strongly on intensity even when the atomic transition is saturated ($`I_s=`$ 0.444 W/cm2 for the 1S0-1P1 laser cooling transition of 24Mg) Suominen98 . In this region the problem is that most theoretical approaches are applicable only to the weak field situation, and can not handle the energy exchange when the energy change spectrum is continuous (as is the case for radiative collisional heating). An important exception is the Monte Carlo wave function method (MCWF) Holland94 ; Piilo01 , which we apply here to the case of 24Mg, to obtain the average energy increase per collision near the Doppler limit (at 1 mK for 24Mg). Although simple in theory, the collisions between the alkaline earth atoms at strong fields are complicated since the partial waves and rotational states of the ground and excited states (respectively) of the quasimolecule form a network of two infinite sets of coupled states Machholm01 .

In this article we report the results of Monte Carlo simulations of collisional radiative heating, based on a truncated set of partial waves and rotational states. We derive a heating rate due to these collisions, and compare it to the photon scattering heating rate that sets the Doppler limit in the standard collisionless theory. We also outline a channelling mechanism which is possible due to the specific structure and selection rules of the alkaline earth dimer quasimolecule, and which could lead to the suppression of inelastic processes, and test it with the simulations. Here we concentrate on the results provided by the simulations, and report the multitude of technical simulation details elsewhere.

## II Magnesium quasimolecule

The lasers for cooling and trapping in magneto-optical traps (MOT) are normally detuned a few atomic linewidths ($`\gamma _{\mathrm{at}}`$) below the atomic transition, which means that the attractive quasimolecule states have a finite internuclear distance, the Condon point $`R_C`$, where they become resonant with the laser field. For 24Mg and other alkaline earth atoms we have two such states, a strongly coupled $`{}_{}{}^{1}\mathrm{\Sigma }_{u}^{+}`$ state, and a weakly coupled $`{}_{}{}^{1}\mathrm{\Pi }_{g}^{}`$ state Stevens77 ; Czuchaj01 .

In the radiative heating mechanism, the colliding atoms on a partial wave $`l`$ are excited at $`R_C`$ to an attractive quasimolecule state with rotational quantum number $`J`$. The kinetic energy of the relative atomic motion increases as the atoms approach each other further, until they return to the ground state via spontaneous emission. The energy exchange due to such a process forms a continuous distribution below the trap depth energy (above which the radiative heating becomes radiative escape).

Since the $`{}_{}{}^{1}\mathrm{\Pi }_{g}^{}`$ state is only weakly coupled to the light field, and provides mainly a channel for trap loss, we can concentrate in our modelling of heating on the $`{}_{}{}^{1}\mathrm{\Sigma }_{u}^{+}`$ state only (see Ref. Machholm01  for a more detailed discussion). The ground state partial waves and the excited rotational states are coupled by the selection rules $`lJ=l\pm 1`$, i.e., the $`Q`$ branch for molecular excitation is missing. In the weak field limit one can consider each partial wave as paired independently to a rotational state, forming a set of independent two-state systems. At strong fields multiple couplings are possible, and one can move back and forth on the “angular momentum ladder” shown in Fig. 1. Note that for the bosonic 24Mg atoms the exchange symmetry allows only even values for $`l`$.

Channelling. One should note, as Fig. 2 illustrates, that each pair ($`l,J`$) is coupled resonantly at a different value of $`R_C`$. Thus at modest field strengths one could consider a situation where the separation of the Condon points is valid while the excitation probability is close to unity. Due to the missing $`Q`$ branch one can see that except for the $`s`$-wave, we always have a sequence $`lJ=l1l^{}=J1=l2`$ so that an incoming ground state partial wave $`l`$ is channelled via two crossings into another ground state partial wave $`l^{\prime \prime }=l2`$, then reflected by the centrifugal barrier, returning the same route back to the original state $`l`$. This process would reduce the collisions practically elastic (the $`s`$-wave still remains ”unchannelled”, though).

## III Collisional heating rate

The energy increase from a single collision connects to a heating rate and the appropriate rate coefficient in the following way. The collisional heating rate $`\kappa _H(n)`$ describes the kinetic energy change per unit volume and unit time. For identical particles it will be equal to $`\frac{1}{2}K_Hn^2`$, where $`K_H`$ is the rate coefficient for collisional heating, and $`n`$ is the atomic density. The factor of $`1/2`$ removes the doubling in collision counting for identical particles. The total collisional heating rate is obtained by integrating $`\kappa _H`$ over the trap volume.

For trap loss rate coefficient $`K_{\mathrm{loss}}`$ one normally calculates the collision frequency times the probability for loss, i.e., $`v\sigma _{\mathrm{loss}}`$, where $`v`$ is the relative velocity, and $`\sigma _{\mathrm{loss}}`$ is the cross-section for collisions, including the probability for a loss event to occur. This simple classical picture is connected to thermodynamics by assuming a thermal equilibrium distribution of velocities, $`f(v)dv`$, over which we take an average, and obtain $`K_{\mathrm{loss}}=v\sigma _{\mathrm{loss}}`$. The dynamics of the two-body collision on the microscopic level then enters in calculating $`\sigma _{\mathrm{loss}}`$.

To calculate the heating rate coefficient $`K_H`$ is slightly more complicated, because in addition to the probability of an inelastic collision event to happen, we also need to estimate the amount of kinetic energy increase associated with it. Technically, we would have a continuous distribution of final energies corresponding to each initial value of $`v`$, due to the randomness of the spontaneous emission events. For practical reasons we consider an averaged rate, i.e., we calculate the average change in kinetic energy per unit time. Since we can perform the averaging over final energy states before the averaging over initial states, we define an inelastic (heating) cross section $`\sigma _H(v)`$ (units energy$`\times `$length2) that gives the difference of the average final relative kinetic energy and initial relative kinetic energy.

In the partial wave approximation we can write the quantum mechanical cross-section for identical atoms in a three-dimensional trap as Julienne89

$$\sigma (v)=\frac{\pi }{k^2}\underset{l=0(\mathrm{even})}{\overset{\mathrm{}}{}}(2l+1)P_l(v),$$
(1)
where $`k`$ is the wave number related to $`v`$ and $`P_l(v)`$ is the event probability. Thus, as a generalization, we write for heating

$$\sigma _H(v)=\frac{\pi }{k^2}\underset{l=0(\mathrm{even})}{\overset{\mathrm{}}{}}(2l+1)\mathrm{\Delta }E_{\mathrm{sc}}(v,l),$$
(2)
where $`\mathrm{\Delta }E_{\mathrm{sc}}(v,l)`$ is the average single-collision energy increase related to the initial partial wave $`l`$.

In principle one should take the thermal average over an isotropic Maxwell-Boltzmann distribution of relative velocities $`v`$. Due to the complexity of the Monte Carlo simulations we are, in practice, limited to calculating $`\mathrm{\Delta }E_{\mathrm{sc}}(v,l)`$ for a rather narrow range of initial relative velocities $`v`$. As in Ref. Machholm01 , we define a single-energy rate coefficient as ($`E=\frac{1}{2}\mu v^2`$)

$$K_H(E)=\frac{E}{hQ_E}\underset{l=0(\mathrm{even})}{\overset{\mathrm{}}{}}(2l+1)\mathrm{\Delta }E_{\mathrm{sc}}(v,l),$$
(3)
with $`Q_E=(2\pi \mu E/h^2)^{3/2}`$, where $`\mu `$ is the reduced two-particle mass. The simulations have shown this to be a reasonable approximation.

## IV Simulation results

We have performed a series of computationally tedious simulations at the detuning $`\delta =3\gamma _{\mathrm{at}}`$ for various laser intensities, and the main results are given in Table 1. We have used the semiclassical trajectory argument (turning point should be at $`R<R_C`$ for any involved partial wave) to truncate the sum in Eq. (3). Clearly a strong increase in the energy change with intensity continues above the saturation intensity $`I_s`$. The semiclassical Landau-Zener theory also indicates that in this intensity region excitation probabilities should be equal to unity. Thus the continued increase can be attributed to the re-excitation of the decayed population, as suggested in Ref. Suominen98 ; the delayed decay model predicts (roughly) $`\mathrm{\Delta }E_{\mathrm{sc}}\sqrt{I/I_s}`$.

Additional simulations showed that pre- and post-collision steady state formation is fast and strong for $`I5I_s`$, and makes the energy exchange independent of the initial partial wave, i.e., $`\mathrm{\Delta }E_{\mathrm{sc}}(v,l)\mathrm{\Delta }E_{\mathrm{sc}}(v)`$ in Eq. (3). In other words, in a system where the energy states are strongly coupled and remain so even asymptotically, the definition of an initial state is obscured, and we can actually make meaningful studies because the spontaneous emission always leads to the same steady state, independent of initial $`l`$, well before $`R_C`$ is reached. For $`I5I_s`$ there is some initial state dependence, but on a scale that allows one to use the values given in Table 1 as a reasonable estimate of magnitude.

As for the channelling effect, there was no sign of it in the simulations. The explanation is rather clear. The Condon points are close to each other both in position and energy, and the true picture of the excitation process is more like adiabatic following of the field-dressed (adiabatic) quasimolecule states, shown in Fig. 3. However, this may not be the case for other values of detuning, or weaker fields than those used in our simulations. For this reason we want to emphasise the general possibility for channelling, and also because it is a good test for the validity of the reduction into independent two-state models Machholm01 . We have confined our computationally demanding simulations to large intensities because then there are suitably many quantum jumps for having reasonable statistics for small ensembles of 64 or 128 members.

Finally it should be noted that any truncation of the partial wave/rotational state manifolds must be done so that the total number of states is even. Otherwise one would introduce an artificial dark state as one of the field-dressed eigenstates, into which the strong and fast steady state formation would trap the system well before the actual collision begins.

## V Comparison of rates

The heating rate prediction of the simulations must be compared to the photon scattering heating rate that enters the calculation of the Doppler limit. In this we follow Ref. Lett89  and obtain the approximate rate

$$\left(\frac{dE}{dt}\right)_{\mathrm{heat}}=\frac{E_r\gamma _{\mathrm{at}}}{2}\frac{I/I_s}{1+I/I_s+(2\delta /\gamma _{\mathrm{at}})^2},$$
(4)
where $`E_r=\mathrm{}^2k_r^2/(2\mu )`$ is the photon recoil energy, and $`I_s=\pi \gamma _{\mathrm{at}}hc/(3\lambda ^3)`$ is the saturation intensity Metcalf03 . Adding the numbers for 24Mg, we get the heating rate per particle ($`II_s`$)

$$\left(\frac{dE}{dt}\right)_{\mathrm{heat}}\frac{E_r\gamma _{\mathrm{at}}}{80}60\frac{\mathrm{K}}{\mathrm{s}}.$$
(5)
The above expression has been derived for a one-dimensional two-beam setup, and its extension to the full six-beam three-dimensional situation is an interesting and open question, and one of the strong motivations to study Doppler cooling with alkaline earth atoms.

The collisional heating rate coefficient becomes, using Eq. (3),

$$K_H(E)=4\pi ^{3/2}\mathrm{\lambda ̄}^3\frac{E_r^2}{h}\left(\frac{E}{E_r}\right)^{1/2}\left(\frac{3\gamma _{\mathrm{at}}}{2|\delta |}\right)^{2/3}\frac{\mathrm{\Delta }E_{\mathrm{sc}}}{E_r}.$$
(6)
Here $`\mathrm{\lambda ̄}=k_r^1=\lambda /(2\pi )`$, with $`\lambda =285.21`$ nm for 24Mg. This result is obtained if we note that by truncating the partial wave series to the maximum classically allowed $`l`$, denoted by $`l_{\mathrm{max}}`$, we have $`E=\mathrm{}^2l_{\mathrm{max}}(l_{\mathrm{max}}+1)/(2\mu R_C^2)`$. The Condon point $`R_C`$ is determined by the $`{}_{}{}^{1}\mathrm{\Sigma }_{u}^{+}`$ state potential, $`U(R)=3\mathrm{}\gamma _{\mathrm{at}}/[2(k_rR)^3]`$, being equal to $`\mathrm{}\delta `$ at $`R_C`$. Finally, $`_{l=0(\mathrm{even})}^{l_{\mathrm{max}}}(2l+1)=(l_{\mathrm{max}}+1)(l_{\mathrm{max}}+2)/2l_{\mathrm{max}}(l_{\mathrm{max}}+1)/2=\mu R_C^2E/\mathrm{}^2`$. To compare with the single-particle rate (5) we write (for $`\delta =3\gamma _{\mathrm{at}}`$)

$$\left(\frac{dE}{dt}\right)_{\mathrm{coll}}=\frac{1}{2}K_H(E)nk_B\times 1.3\times 10^{14}\frac{\mathrm{K}}{\mathrm{s}}\times \frac{\mathrm{\Delta }E_{\mathrm{sc}}}{E_r}\times n.$$
(7)
For $`n10^{11}`$ atoms/cm3 and $`II_s`$ this gives about 0.1 K/s. This result implies that radiative collisional heating can not explain the problems in reaching the Doppler limit for 24Mg, because the densities are usually in the range of $`10^9`$ to $`10^{11}`$ atoms/cm3. We have earlier performed two-state studies for Cs and found that there the heating rate present in Sisyphus cooling should match the collisional heating rate at $`n=10^{12}`$ atoms/cm3 Holland94 .

As for making a comparison to the 88Sr case, without performing simulations, we can make a rough estimate. In Eq. (3) we have written $`K_H(E)`$ in a scaled form, so we can put in the corresponding values for 88Sr as well ($`E/k_BT_D=0.77`$ mK, $`\lambda =460.73`$ nm, $`E_r=1\mu `$K, and taking again $`|\delta |=3\gamma _{\mathrm{at}}`$). The remaining problem is to estimate $`\mathrm{\Delta }E_{\mathrm{sc}}`$. In a toy model that assumes total excitation during a collision, and ignores re-excitation, we can write $`\mathrm{\Delta }E_{\mathrm{sc}}U^{}(R_C)\tau v`$, where $`U^{}(R_C)`$ is the slope of the excited state (ignoring $`J`$) potential at $`R_C`$, and $`\tau 1/\gamma _{\mathrm{at}}`$ is the average survival time on the excited molecular state, and $`v`$ is the initial collision velocity. This leads to an estimate of $`\mathrm{\Delta }E_{\mathrm{sc}}/E_r=23(E/E_r)^{1/2}`$, which at least for 24Mg gives 230, i.e., the correct magnitude near $`I=I_s`$. With all this, we get for 88Sr at the Doppler temperature the result $`K_H(E)2\times 10^{12}`$ Kcm3/s, to be compared with $`E_r\gamma _{\mathrm{at}}/802.5`$ K/s. These results meet only at densities of $`n10^{12}`$ atoms/cm3, which are a few magnitudes higher than those used in the experiment of Ref. Xu02 . This toy model also predicts a detuning dependence of $`K_H(E)`$ of $`|\delta |^{2/3}`$, i.e., heating decreases with detuning, which means that while the number of involved partial waves increases with diminishing $`|\delta |`$, the effect is more than countered by the reduction in the steepness of the excited state potential.

## VI Conclusions

This study is the first attempt to describe quantitatively the radiative collisional heating by taking into account the inherent multistate structure due to the partial waves and rotational states. The truncation of the partial wave manifold has required careful testing of various possibilities. Although it appears that collisional radiative heating does not play a role in the thermodynamics of the trap at realistic densities for 24Mg (and other alkaline earth atoms), our work does not rule it out completely. Although assumed to be irrelevant here, the weakly coupled $`{}_{}{}^{1}\mathrm{\Pi }_{g}^{}`$ state may play a role, as well as the two quasimolecule states corresponding to the repulsive potentials. Experimentally the role of collisional processes can be distinguished from other processes by looking at the scaling with density; this has already been used to separate the cold collision contributions from the background collisions in trap loss. While we aim to improve our modelling and to use our heating results in estimating radiative trap loss, we also expect that more experimental data will become accessible in the near future. It should be pointed out that other effects may also lead to heating. The methods to study e.g. radiative trapping (reabsorption of scattered photons) are, however, very different from the ones employed here.

## VII Acknowledgments

The authors acknowledge the EU network CAUAC (contract HPRN-CT-2000-00165) and the Academy of Finland (projects 50314 and 206108) for financial support. JP has also been supported by the EU network COCOMO (contract HPRN-CT-1999-00129) and the Academy of Finland (project 204777). Discussions with E. Rasel, J. Thomsen and J. Ye are gratefully acknowledged, as well as useful comments by I. Mazets and J. Vigué. The Finnish IT Center for Science (CSC) is gratefully acknowledged for providing the supercomputer environment used in the simulations.


text:
# Tests of silicon sensors for the CMS pixel detector

## 1 Introduction

The CMS experiment at the Large Hadron Collider (LHC) at CERN (Geneva, Switzerland) will be equipped with a silicon pixel detector providing high resolution 3D coordinates of the tracks. It will consist of three layers of pixel sensors (100 m $`\times `$ 150 m pitch size) in the barrel region and two discs in the forward regions . The innermost barrel layer is expected to be exposed to a fluence111All particle fluences are normalized to 1 MeV neutrons ($`\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$). of $`3\times 10^{14}\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$ per year at full luminosity. The second and third layer will be exposed to $`1.4\times 10^{14}\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$ and $`0.6\times 10^{14}\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$ per year, respectively. The charge sharing induced by the Lorentz force in the 4 T CMS detector magnetic field can improve the spatial resolution. Hence it is very important to study the effects of charge collection efficiency and Lorentz drift in irradiated silicon devices. Here we present measurements of Lorentz angle and charge collection efficiency, and compare two sensor designs at different irradiation doses.

## 2 Experimental setup

The measurements were performed in the H2 beam line of the CERN SPS in June and September 2003 using 150-225 GeV pions. A silicon reference telescope  was used to allow a precise determination of the particle hit position in the pixel detector. Both, the telescope modules and the pixel front-end, were mounted onto a common frame (see fig. 1). The beam telescope consists of 4 modules. Each consists of two 300 m thick single-sided silicon detectors ($`32\times 30`$ mm2) with a strip pitch of 25 m (readout pitch of 50 m), which are oriented perpendicularly. The resulting intrinsic resolution of the beam telescope is around 1 m.

The pixel sensor with the readout chip is mounted on a rotating support positioned between the second and the third module. A trigger signal is generated by a silicon PIN diode. The data acquisition system, including slow control (temperature and bias voltage of the pixel sensor) is written in LabView and LabWindows CVI (National Instruments) running on a PC. The analog signals are digitized in a VME based readout system by two CAEN (V550) and one custom built FADCs . The whole setup was placed in a open 3 T Helmholtz magnet with magnetic field parallel to the beam. The pixel sensors were cooled down up to $`30^{}`$ C by means of water cooled Peltier elements.

## 3 Pixel sensors and front-end electronics

The CMS pixel sensors are manufactured in the ”n-on-n” technique, consisting of $`\mathrm{n}^+`$ structures on n bulk silicon. This allows a partially depleted operation of highly irradiated sensors after type inversion, but requires inter-pixel isolation. Two isolation techniques were considered in our latest prototype designs: p-spray , where a uniform medium dose of p-impurities covers the whole structured surface, and p-stop, where higher dose rings individually surround the $`\mathrm{n}^+`$-implants. Fig. 2 shows the layouts of these designs. Due to possible failures in the pixel contacts (bump-bonding) an intrinsic biasing method has to be introduced. This is realized by a bias grid and punch-trough structures  for the p-spray detectors and by openings in the rings for the p-stop detectors.

All our test devices had $`22\times 32`$ pixels, with a sensitive area of $`2.75\times 4`$ mm2, and a thickness of 280 m. The readout pitch was 125 m $`\times `$125 m.
Some of the devices were irradiated in a 24 GeV proton beam at the CERN PS. The sensors received total particle fluences of $`3.3\times 10^{14}\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$, $`8.1\times 10^{14}\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$ and $`1.1\times 10^{15}\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$ respectively. During the irradiation the sensors were kept at room temperature and without bias, while after irradiation we stored and operated them at $`20^{}`$ C. $`IV`$-curves were measured  before and after irradiation, all measured sensors showed a breakdown voltage above 600 V. Finally the sensors were bump bonded to non irradiated readout chips of the type PSI30/AC30  allowing for a non zero-suppressed analog readout of all 704 pixel cells. A fast ($`20`$ ns) external signal for the internal sample and hold mechanism had to be provided from the trigger PIN diode.

## 4 Beam test measurements

### 4.1 Lorentz angle measurements

The Lorentz angle is obtained by the direct measurement of the charge drift in the magnetic field using the grazing angle method .

The particle beam (parallel to the magnetic field) hits the pixel surface at a shallow angle $`\alpha =15^{}`$ (see fig. 3). The deposited charge drifts according to the combined electric and magnetic forces, resulting in a deflection of the particle track projection on the surface by the angle $`\beta `$. The Lorentz angle is obtained from the relation $`\mathrm{tan}\mathrm{\Theta }_L=\mathrm{tan}\beta /\mathrm{tan}\alpha `$.
The deflection measured at 3 T magnetic field is shown in fig. 4 for the non-irradiated p-spray sensor.

The angle $`\beta `$ is measured by slicing the histogram perpendicularly to the $`y`$-axis. The position of the center of each slice is shown in fig. 5 as a function of the position along the $`y`$-axis. A measurement without magnetic field is used to correct for the detector misalignment with respect to the beam (bottom line in fig. 5). For the irradiated device in fig. 5 there are two regions with a different slope, and hence, two different values of the Lorentz angle. The difference in the Lorentz angle can be explained by nonuniform distribution of the electric field in the irradiated devices. Different models of the electric field distribution in the irradiated silicon devices are discussed in Refs. .

The evidence of the non-uniform electric filed will be shown also in the next section.
As most of the signal charge is collected from the region close to the pixel implant, this region was used to determine the Lorentz angle. Figure 6 shows the measured values for both non-irradiated and irradiated sensors extrapolated to 4 T magnetic field.

We observe a strong dependence of the Lorentz angle on the bias voltage (electric field), while it is only weakly affected by the irradiation or sensor design. For the non-irradiated sensors a Lorentz angle of $`26^{}`$ can be reached at a bias voltage of $`100`$ V, while irradiated sensors have to be operated at higher bias voltages, where the Lorentz angle drops to roughly $`10^{}`$. Since the electron mobility increases with decreasing temperature the Lorentz angle measured at the lower temperature $`20^{}`$ C is $`1^{}2^{}`$ larger. Our values for the Lorentz angle are in a good agreement with measurements and simulations in Refs. .

### 4.2 Charge collection in irradiated sensors

After irradiation the collected charge decreases due to charge trapping and partial depletion of the sensor. The measurements of the charge collection efficiency as a function of the sensor depth were performed using the grazing angle method without magnetic field. The averaged charge cluster profiles in non-irradiated and irradiated p-spray sensors are shown in fig. 7. In the unirradiated sensor the charge is collected uniformly across the whole sensor depth, while the irradiated devices shows two regions with different collected charge. At low bias voltage in the irradiated detectors some charge is also collected from the side opposite to the pixel implants. This can indicate that the depletion starts from both sides of the detector, because of the nonuniform distribution of the electric filed.

The integral of the collected charge from the particle entry point up to a point along the $`y`$-coordinate is shown in fig. 8. There are two regions with different charge collection, represented by two different slopes of the curve for the irradiated device.

The total collected charge for different designs and irradiation doses as a function of bias voltage is shown in fig. 9.

### 4.3 Signal-to-noise ratio

The charge collection depends on the position of the incident particle with respect to the pixel implant, caused by the regions with reduced sensitivity. The average charge collected in the hit pixel as a function of the hit position is shown in fig. 10 for non-irradiated and irradiated sensors. One can see that after irradiation the average pixel signal decreases, and the area with the reduced charge collection increases in the case of the p-stop design.

The signal-to-noise ratio (S/N) averaged over the whole pixel area is shown in fig. 11. The measured signal-to-noise ratio is above 30 even after the irradiation. We observe lower S/N ratio in the unirradiated p-stop device, because of the charge spread over several pixels caused by the p-stop openings. A more detailed study of charge collection for these sensors can be found in Ref. .

## 5 Conclusions

Two different designs of the CMS prototype pixels sensors were tested up to 600 V bias voltage and after exposure to particle fluences up to $`1.1\times 10^{15}\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$. For the irradiated devices we observe two regions across the sensor with different charge collection and value of the Lorentz angle. This result is an indication of a nonuniform electric field in the irradiated devices.
The Lorentz angle mainly depends on the bias voltage and it shows a weak dependence on the irradiation dose or sensor design. The Lorentz angle at 4 T magnetic field reaches 26 for the non-irradiated devices at a bias voltage of 100 V and it drops to 8.3 for the ones irradiated at $`1.1\times 10^{15}\mathrm{n}_{\mathrm{eq}}/\text{cm}^2`$ and at a bias voltage 600 V.
After irradiation and at high bias voltage (600 V) the charge collection is about 60% of the value observed for the non-irradiated devices, and it is slightly larger for the p-spray design.
The signal-to-noise ratio decreases from 65 to about 35 after irradiation giving us confidence in operating the CMS pixel detector up to the maximum expected irradiation dose.
For the irradiated devices the charge collection, and therefore the signal height, increases with the bias voltage but the Lorentz angle decreases. For this reason the bias voltage must be optimized for the best performance of the CMS pixel detector.

## Acknowledgments

We gratefully acknowledge Silvan Streuli from ETH Zurich and Fredy Glaus from PSI for their immense effort on the bump bonding, Maurice Glaser and Michael Moll from CERN for carrying out the irradiation, Kurt Bösiger from the Zürich workshop for the mechanical construction, György Bencze and Pascal Petiot from CERN for the H2 beam line support and finally the whole CERN-SPS team. The authors would like to thank Dr. Danek Kotlinski from PSI for very useful discussions.


text:
# Circuit with small-capacitance high-quality Nb Josephson junctions

(September 19, 2003)

## Abstract

We have developed a fabrication process for nanoscale tunnel junctions which includes focused-ion-beam etching from different directions. By applying the process to a Nb/(Al–)Al2O3/Nb trilayer, we have fabricated a Nb single-electron transistor (SET), and characterized the SET at low temperatures, $`T=0.0440`$ K. The superconducting gap energy and the transition temperature of the Nb SET agree with the bulk values, which suggests high quality Nb junctions. The single-electron charging energy of the SET is estimated to be larger than 1 K.
Appl. Phys. Lett. 84, 410 (2004) \[DOI: 10.1063/1.1640798\]

Circuits with Josephson tunnel junctions are one of the most promising candidates for quantum bits (qubits) for quantum computation in solid-state electronic devices.Makhlin et al. (2001); Fitzgerald (2002) The Josephson junctions (JJs) for qubits should have large subgap resistance and an appropriate $`E_J/E_C`$ ratio, where $`E_J`$ is the Josephson energy and $`E_C`$ is the charging energy of the junction. The desired range of $`E_J/E_C`$ depends on the degree of freedom (charge, flux, or phase) to be controlled in the circuit. For cases of charge qubitNakamura et al. (1999) and flux qubit,Chiorescu et al. (2003) a typical junction size is on the order of $`0.1\times 0.1`$ $`\mu `$m2, for which the fabrication technique of Josephson tunnel junctions has been well established only for Al/Al2O3/Al junctions. However, materials with larger superconducting gap energy $`\mathrm{\Delta }`$, e.g., Nb ($`\mathrm{\Delta }_{\mathrm{Nb}}/\mathrm{\Delta }_{\mathrm{Al}}8`$), are more attractive especially for flux qubits, where one wants to have large $`E_J/E_C`$ and keep $`E_C`$ sufficiently larger than the thermal energy $`k_BT`$. Note that $`E_J`$ is proportional to $`\mathrm{\Delta }`$.

There have been a number of attempts to fabricate small-capacitance Al/Al2O3/Nb or Nb/(Al–)Al2O3/Nb JJs.Martinis and Ono (1990); Harada et al. (1994); Pavolotsky et al. (1999); Kim et al. (2002); Dolata et al. (2002, 2003) Conventional nanoscale fabrication processes based on $`e`$-beam lithography and multiangle shadow evaporation, which work fine for Al/Al2O3/Al JJs, tend to deteriorate the quality of Nb, i.e., $`\mathrm{\Delta }`$ and the superconducting temperature $`T_c`$ are reduced considerably (a summary of this can be found in Fig. 1 of Ref. Dolata et al., 2003). Rather large values of $`\mathrm{\Delta }_{\mathrm{Nb}}/e`$ have been obtained in single junctions fabricated by a sloped-edge techniqueMartinis and Ono (1990) ($`11.2`$ mV) and in a single-electron transistor (SET) fabricated by a multilayer techniquePavolotsky et al. (1999) (1.35 mV). However, it would be difficult to fabricate multijunction devices using the sloped-edge technique, and the charging energy of the SET fabricated by the multilayer technique was $`0.15`$ K, which might be acceptable for flux qubits but would be too small for charge qubits and for most single-electron-tunneling devices.

We have developed a process to fabricate high-quality small-capacitance Nb/(Al–)Al2O3/Nb JJs which includes focused-ion-beam (FIB) etching from two different directions. By employing the process, we have fabricated SETs with a three-dimensional structure, shown in Fig. 1.

Details of the process are as follows. A trilayer of Nb (thickness: 0.3 $`\mu `$m), Al (0.01 $`\mu `$m)–Al2O3 ($``$1 nm), and Nb (0.3 $`\mu `$m) was sputtered onto a SiO2/Si substrate in a single vacuum cycle. The Al2O3 layer was formed by oxidizing the surface of Al. The film was then patterned with standard photolithography and Ar+ milling, where the smallest feature size in this step was $``$$`\mu `$m. The dimensions of the trilayer were decreased further in a Ga+ FIB system, which has three functions: deposition, etching, and observation. From the direction perpendicular to the substrate, we first deposited C (thickness: $``$0.1 $`\mu `$m) by decomposing phenanthrene (C14H10) gas with a Ga+ beam current of 48 pA. The role of the C layer is to minimize damage of the trilayer by Ga+ during the etching process described below. Then, we etched into the pattern shown in Fig. 1(a) except for the fine structure on the narrow track between the source and the drain. For etching, we used two beam currents, 1.3 nA (rough etching) and 9 pA, not only for efficiency but also to minimize the amount of material redeposited on the sides of the narrow track. After tilting the substrate $``$88, three holes were made with a 9 pA beam current on the side of the track as shown in Fig. 1(b). A double-junction structure formed between the holes. An example of the SET fabricated in the above process is shown in Fig. 1(c), which is a secondary-electron image taken in the same FIB system. The final step in our fabrication process is anodization, and there are two reasons for this. One is to eliminate completely short circuit of the junction due to conducting materials (mostly Nb) redeposited during FIB etching, by changing the conducting materials into insulating oxides. The other is to suppress the contribution of the sample surface, which might have been damaged by Ga+, to electrical conduction. Because anodization also reduces the effective area of the junctions, we designed the initial junction size delineated by the FIB to be $`>0.2\times 0.2`$ $`\mu `$m2, although we could fabricate much smaller junctions by this technique. We kept the anodization current constant at a value in the range 1–10 $`\mu `$A/cm2 and, for the sample discussed below, we anodized up to 60 V, which reduced the effective junction area to $`<0.1\times 0.1`$ $`\mu `$m2. The FIB etching was quite reproducible, although the resistance of the SET after anodization varied from sample to sample.

We measured the samples in a 4He continuous-flow cryostat with a Si-diode thermometer mounted next to the sample holder. The current-voltage ($`I`$-$`V`$) characteristics of one of the samples at $`T=3.2`$ K is shown in Fig. 2.

The $`I`$-$`V`$ curve exhibits a sharp superconducting gap, whose value corresponds to two Nb JJs. We obtained $`\mathrm{\Delta }_{\mathrm{Nb}}/e=1.4`$ mV from the distance between the peaks in the differential conductance $`dI/dV`$ vs $`V`$ curve (see the lower data set of Fig. 2). Here the charging energy of the SET is negligibly small for the estimate of $`\mathrm{\Delta }_{\mathrm{Nb}}`$ as will be discussed later. This value of 1.4 mV is the largest for $`0.1\times 0.1`$ $`\mu `$m2 junctions. Moreover, it agrees with the values at liquid 4He temperatures for 1–10 $`\mu `$m scale JJs fabricated by a standard photolithographic technology for integrated JJ circuits. We also measured the temperature dependence of the $`I`$-$`V`$ curve in order to determine $`T_c`$. The $`I`$-$`V`$ curve looks almost linear at $`T>8`$ K. However, when we look at the differential resistance $`dV/dI`$, it is easy to find a qualitative difference between the $`dV/dI`$ vs $`V`$ curves at $`T9.0`$ K and those at $`T9.2`$ K. The inset of Fig. 3 shows the $`dV/dI`$ vs $`V`$ curve at $`T=9.0`$ and 9.2 K.

The curve at 9.0 K has a dip in the middle, while that at 9.2 K does not. Based on this qualitative difference, we determine that $`T_c=9.1\pm 0.2`$ K, which agrees with the bulk value of 9.2 K within the error. Here $`\pm `$0.2 K includes all errors in thermometry and $`T_c`$ determination. We confirm the adequacy of the above $`T_c`$ determination in Fig. 3 by plotting versus $`T`$ the zero-bias resistance for $`T9.2`$ K (normal state) and the subgap resistance for $`T9.0`$ K (superconducting state), where we have defined the subgap resistance as the maximum value of $`V/I(V)`$. From the values of $`\mathrm{\Delta }`$ and $`T_c`$, we conclude that the quality of our JJs is high and it has not deteriorated in the FIB etching process.

Let us look at the properties of the sample as an SET. The gate electrode is located $`>`$$`\mu `$m from the double-junction system (Fig. 1), and the voltage $`V_g`$ applied to the gate modulates the current even up to 5 K (data not shown). Thus, the sample indeed works as an SET. In order to characterize the sample further, we cooled it down to 0.04 K in a 3He–4He dilution refrigerator. At $`T1`$ K, the subgap resistance becomes $`>`$102 M$`\mathrm{\Omega }`$, which is significantly larger than the normal-state resistance, 0.3 M$`\mathrm{\Omega }`$, and again, suggests high quality junctions.

The parameters of a SET, such as the capacitances $`C_1`$ and $`C_2`$ of the tunnel junctions and the capacitance $`C_g`$ between the island electrode and the gate electrode, can be determined by measuring the blockade of single-electron tunneling in the normal state at a low enough temperature.Ingold and Nazarov (1992) We drove the Nb SET into the normal state by applying a magnetic field of 2.7 T perpendicular to the substrate, and measured $`I`$-$`V_g`$ curves for different values of $`V`$ at $`T=0.05`$ K. An example of the $`I`$-$`V_g`$ curve, which is for $`V=0.026`$ mV, is shown in the lower data set of Fig. 4.

One period in the $`V_g`$ axis corresponds to $`e/C_g`$, so that $`C_g=7`$ aF, which is consistent with the geometry. In the upper data set of Fig. 4, we estimate the zero-current region at $`T=0`$ on the $`V`$-$`V_g`$ plane (Coulomb diamond). The horizontal dotted line corresponds to $`2E_C/e`$ of the SET. Thus, from $`2E_C/e=0.18`$ mV in Fig. 4, we obtain $`E_C/k_B=1.1`$ K, which is consistent with the effective junction area of $`<0.1\times 0.1`$ $`\mu `$m2 and sufficiently larger than the base temperature of a typical dilution refrigerator, $`<`$0.1 K. When $`C_gC_1,`$ $`C_2,`$ which is the case in our Nb SET, $`C_1`$ and $`C_2`$ are estimated from the slopes of the solid lines in Fig. 4. We find that $`C_1/C_2=0.75`$, although we designed it so that $`C_1=C_2`$.

The reason for asymmetry in the junction capacitance is that we had to anodize the SET up to a large voltage in order to eliminate the short circuit. The accuracy of the capacitances could be improved by introducing a step of reactive-gas-assisted etching at the end of the FIB process. We have already confirmed that XeF2 gas enhances the etching rate of Nb $``$102 times and “cleans up” redeposited Nb. By introducing the step, we would be able to reduce the anodization voltage considerably. Reducing the anodizaiton voltage would also improve the yield of the fabrication process.

It should be noted that our process is also applicable to more complex circuits. Moreover, it is much more flexible in terms of circuit pattern than the conventional technique based on $`e`$-beam lithography and shadow evaporation, or the multilayer technique in Ref. Pavolotsky et al., 1999. Very recently, a similar FIB-etching technique was independently developed, and nanoscale single junctions with a variety of materials were fabricated. Bell et al. (2003) Thus, the process is not limited to Nb/(Al–)Al2O3/Nb junctions.

In the superconducting state, one expects that the supercurrent flowing through the SET depends on $`V_g`$ periodicallyJoyez (1995) and that the period is $`2e/C_g`$ ($`2e`$ periodic). In many experiments, however, a period of $`e/C_g`$ ($`e`$ periodic), which suggests the existence of subgap quasiparticle states,Joyez (1995) has also been seen. In small-capacitance Al/Al2O3/Nb or Nb/(Al–)Al2O3/Nb systems, only $`e`$ periodicity has been reported. In our Nb SET, the measured supercurrent at $`T=0.04`$ K and $`B=0`$ is also $`e`$ periodic, and its magnitude is on the order of 10 pA, which is $``$10-3 the theoretical maximum, $`I_0/2`$, where $`I_0\pi \mathrm{\Delta }/2eR_n`$ is the Ambegaokar-Baratoff critical current and $`R_n`$ is the normal-state resistance of the junction. Further investigation of the periodicity by measuring high-quality Nb SETs with different parameters would clarify whether the $`e`$ periodicity is intrinsic to Nb or not.

In summary, we have fabricated a high-quality Nb SET with $`E_C/k_B>1`$ K by developing a fabrication process for nanoscale tunnel junctions. The process is much more flexible than conventional ones based on $`e`$-beam lithography.

The authors are grateful to Y. Kitagawa for preparation of the Nb/(Al–)Al2O3/Nb trilayer, to M. Ishida for assistance with the FIB system, to H. Akaike for valuable comments on anodization, and to Yu. Pashkin, T. Yamamoto, and O. Astafiev for fruitful discussions. This work was supported in part by the Special Postdoctoral Researchers Program of RIKEN and by MEXT.KAKENHI (Grant No. 15740190).


text:
# Light and Heavy Dark Matter Particles

## I Introduction

Dark Matter particles are expected to produce indirect annihilation signatures via production of high energy $`e^+,\overline{p},\gamma ,\nu `$ which are potentially detectable . An especially intriguing signal has recently been observed with the SPI spectrometer on the INTEGRAL (INTErnational Gamma-Ray Astrophysics Laboratory) satellite, which detected a bright 511 keV $`\gamma `$-ray line from the galactic bulge . This detection turns out to be in good agreement with previous measurements  and fits reasonably well to a gaussian with full-width half-maximum of $`9^{}`$, with a 2$`\sigma `$ confidence interval of $`6^{}18^{}`$.

Positrons in the galactic bulge can be emitted by such astrophysical objects as for example hypernovae , neutron stars or black holes  etc. However, whether these sources can be at the origin of this 511 keV line is still under debate , especially in the bulge which consists exclusively of old, low mass stars. As an alternative explanation, annihilations of light Dark Matter particles  into $`e^+e^{}`$ seem to provide a surprisingly good explanation for the Integral measurement, if the dark matter halo profile is approximated in the centre of the galaxy as $`\rho (r)r^\gamma `$ with $`\gamma [0.40.8]`$ . Such a value corresponds to a mild cusp intermediate between the NFW distribution  and the profile extrapolated from observations of dwarf and low surface brightness galaxies . This empirically-derived profile is consistent with a density profile that incorporates information on the dark matter contribution constrained by gravitational microlensing measurements of the Milky Way bulge . Moreover the inferred annihilation cross section is concordant with that inferred from the relic density, without any need to assume a boost factor, due for example to dark matter clumpiness.

One way to determine whether the 511 keV line is due to astrophysical sources or Dark Matter (DM) annihilations is to seek a similar signature from low surface brightness dwarf galaxies . Those objects are everywhere dark matter–dominated. Therefore, if the emission line detected in our galaxy is due to DM annihilations, then one should also detect a 511 keV line from nearby dwarf spheroidals. The flux from the Sagittarius Dwarf Galaxy, in particular, appears to fall within the expected sensitivity of Integral/SPI. It is therefore quite likely that, in the near future, Integral/SPI will be able to rule out or confirm the existence of Light Dark Matter (LDM) particles.

Other experiments seek to detect very high energy gamma rays (i.e. of order a hundred GeV) . Clearly, if the assumption of LDM were to be confirmed, it would seem difficult to accomodate low as well as possible high energy gamma ray signatures (unless, of course, astrophysical sources appear to be at the origin of the high energy gamma rays). In fact, the latter may have a natural explanation within the framework of the (Minimal) Supersymmetric Standard Model. Annihilations of the lightest neutralino are indeed expected to give rise to gamma rays of up to several hundred GeV. Whether they would appear a satisfactory solution or not then depends on the angular and energy distribution of the observed gamma rays.

However, it is difficult to accomodate very large neutralino masses (as required if very energetic gamma rays are detected), with the large annihilation cross sections needed to comply with the sensitivity of present/next generation gamma ray experiments, unless one invokes the supermassive black hole at the centre of our galaxy to steepen the innermost dark matter cusp within the zone of influence of the black hole and thereby enhance the predicted annihilation rate . The amplitude of this effect is of course dependent on uncertain astrophysical assumptions about the formation of the supermassive black hole.

Also, it is worth mentioning that the adoption of very large neutralino masses requires one to invoke the coannihilation mechanism in order to produce a low enough relic density. However, the conditions for coannihilation to operate are not always satisfied, as shown in . In particular, the decay process of the particle with which the neutralino is supposed to coannihilate (namely the Next-to-Lightest Supersymmetric Particle, NLSP) may be more efficient than the processes that can regenerate this particle. As a result, one is left with the annihilation region (allowing for neutralino masses up to $``$ 100 GeV), plus potentially only a small region of the supersymmetric parameter space corresponding to the largest neutralino masses (i.e. for which the neutralino is almost mass degenerate with the NLSP and where the efficiency of the NLSP decay process is significantly reduced).

In this Letter, we point out that it is possible to reconcile the low and high energy signatures, even if both of them turn out to be due to Dark Matter annihilations. Indeed, existing theories predict two different kinds of Dark Matter particles . For example, in a $`N=1`$ supersymmetric framework with mirror fermions (somehow reminiscent of $`N=2`$ extended Supersymmetry) , one expects the existence of two neutral and stable particles. One would be a heavy fermion for example, like the lightest neutralino, and the other one a possibly light spin-0 particle. Both of them would be neutral and also stable as a result of two discrete symmetries (say $`R`$ and $`M`$-parities).

The light particle would be crucial in order to achieve the correct relic density, while heavy stable particles would be at the origin of very high energy gamma rays. Light particles would be difficult to detect in direct Dark Matter detection experiments because of their mass but heavy particles would also be hard to find due to their small residual number densities.

## II Relic density, 511 keV line and high energy gamma rays

Let us assume that there indeed exists two DM particles. The first question to deal with is: how do they share the observed relic density ?

If the 511 keV emission line is indeed due to LDM annihilations in the centre of the Milky Way, then one needs a cross section (times the relative velocity, normalised to $`c`$) of about $`10^4`$ to $`10^5`$ pb to fit the observed flux ($`9.9\times 10^4\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1`$) and the angular distribution of $`\gamma `$-rays indicated by Integral/SPI.

On the other hand, it is also necessary to avoid an overproduction of relic LDM for the observed low energy gamma ray flux or an overproduction of gamma rays for a given relic LDM abundance. Since low DM masses actually increase the amount of annihilations in our galaxy, and consequently the flux of low energy gamma rays that has been well measured by several experiments, one has to invoke a $`v^2`$ suppressed annihilation cross section. One way to achieve this is to assume a scalar particle coupled to a new light gauge boson *** This gauge boson is similar to the one introduced in . However, in the case we discuss here, it would be coupled to Dark Matter and, in fact, would mainly decay into two DM particles..

The LDM cross section would then be much larger at early times than is expected in our galaxy . For example, one expects their annihilation cross section to be about $`10^5`$ pb in the Milky Way, to explain the 511 keV line, and about 1 pb in the primordial Universe (since the DM velocity is $`\mathrm{3\; 10}^3`$ the velocity in the primordial Universe).

Given the relationship between the relic density of annihilating DM particles and the annihilation cross section in the primordial Universe, the cosmological parameter associated with these light particles is then expected to be about

$$\mathrm{\Omega }_{\mathrm{dm}}^{th}h^20.024\times \frac{x_F}{\sqrt{g_{}}}\times \left(\frac{\sigma v_p}{1\mathrm{pb}}\right)^1$$
(1)
where $`x_F`$=$`m_{\mathrm{dm}}/T_F1219`$ for particles in the MeV-GeV range, $`g`$ and $`g_{}`$ are the numbers of internal and relativistic degrees of freedom, respectively and $`v`$ is expressed in units of $`c`$. Note that $`\sigma v_p`$ stands for the primordial annihilation cross section of LDM.

For $``$ MeV particles, the value of $`\frac{x_F}{\sqrt{g_{}}}`$ is $``$ 4. Thus one obtains $`\mathrm{\Omega }_{\mathrm{dm}}^{th}h^20.095\left(\sigma v_p/(1\mathrm{pb})\right)^1`$. This has to be compared to $`\mathrm{\Omega }_{\mathrm{dm}}h^2`$ 0.1 as found by WMAP.

LDM particles are therefore likely to fit the observed relic density, leaving little room for another “weakly annihilating” species. Of course, this depends on the profile used and the exact cross section needed to fit the angular resolution of the SPI detection. If the annihilation cross section of heavy DM particles ($`\sigma v^{}`$) is larger than $`1`$ pb at both early epochs and todayOne would then require constant annihilation cross sections for heavy particles., then these particles are likely to be subdominant but, on the other hand, they would yield significant gamma ray fluxes! This picture would then lead to a situation where one can fit the observed relic density, explain the nature of the 511 keV line if it is confirmed and potentially predict, at the same time, a non-negligible flux for high energy gamma rays, even if the dark halo profile is not as cuspy as the NFW profile.

Assuming that LDM fulfills the relic density condition, one can study the effect of a large annihilation cross section associated with heavy particles. Of course, the larger this cross section gets, the less credible it becomes (as one expects $`\sigma 1/m_\chi ^2`$ or $`1/m_S^2`$, where $`\chi `$ and $`S`$ denote the heavy DM particle or the exchanged particle respectively). However, it is worth remembering that resonances are possible, in which case they would significantly enhance the cross section.

The 511 keV emission line seems to favor a quite flat profile (for example $`\gamma 0.4`$). If so, then the annihilation cross section needed may appear extremely large. In fact, this would be probably too large compared to what could be realistically expected in a particle physics model unless, perhaps, new kinds of interactions are invoked. On the other hand, even a quite mild dark halo profile (i.e. $`\gamma 0.8`$) would yield a significant amount of high energy gamma rays if the annihilation cross section associated with heavy DM particles turns out to be large but still “realistic” (i.e. enhanced by resonance effects for example). Such a scenario would then be in agreement with the Dark Matter cosmological parameters measured by WMAP as well as with low and high energy annihilation signatures (if they were both confirmed).

## III A theory with two stable Dark Matter particles

A very important requirement regarding any Dark Matter candidate is its stability. In a Standard Supersymmetric framework, the neutralino (absolute) stability is obtained through a $`Z_2`$ discrete symmetry, called $`R`$-parity . This symmetry implies the decay of a supersymmetric particle into another one, plus standard model particles. This is what prevents the lightest supersymmetric particle (LSP) from decaying.

Since we are now looking for two stable particles, it is particularly interesting to start working within the framework of $`N=2`$ extended supersymmetric theories, as they may naturally allow for two $`Z_2`$ discrete symmetries, and may ultimately lead to two new distinct DM candidates. We shall denote these two discrete symmetries as $`R`$ and $`M`$ parities.

Note that the particles which matter for our purpose (say mirror leptons/quarks and new spin-0 states that couple these mirror particles to ordinary ones ) may be kept even after one abandons the full $`N=2`$ supersymmetry in favor of a $`N=1`$ theory.

### A Two discrete symmetries

As mentioned before, $`N=1`$ supersymmetric theory may admit a discrete $`Z_2`$ symmetry. The latter can be seen, in fact, as a remnant of a continuous $`R`$-symmetry ($`R_p=(1)^R`$).

Thus, quite similarly to the case $`N=1`$, an extended $`N=2`$ theory, for which there are now two supersymmetry generators

$$Q=(\begin{array}{c}Q_L^1\\ Q_L^2\end{array}),$$
(2)
may admit a $`SU(2)\times U(1)`$ global $`R`$-symmetry, that acts on the doublet $`Q`$.

Such a global $`R`$-symmetry allows one to act independently on the two supersymmetry generators $`Q_1`$ and $`Q_2`$. One can therefore define two distinct $`R`$-parity symmetry operators, $`R_{1p}`$ and $`R_{2p}`$; the usual $`R`$-parity symmetry corresponding to the product

$$R_p=R_{1p}R_{2p}.$$
(3)
By abandoning the full $`N=2`$ supersymmetry, one singles out one of the two supersymmetry generators (say $`Q_1`$) and redefines one of the two discrete symmetries (say $`R_{2p}`$) as the $`M`$-parity symmetry.

We are then left with two discrete symmetries: the $`M`$-parity that can be seen as a $`Z_2`$ remnant of a continuous $`U(1)`$ symmetry ($`M`$) acting globally on $`N=1`$ superfields, and the “standard” $`R`$-parity. Both of them are kept intact after spontaneous breaking of the electroweak symmetry.

### B Particle spectrum from extended supersymmetry

We now have to specify what is the expected spectrum from extended supersymmetry and, more precisely, which particle could be the LMP.

When $`N=2`$ extended supersymmetric theories are formulated in terms of $`N=1`$ superfields, ordinary gauge superfields are accompanied by additional ($`N=1`$) “chiral gauge superfields”. Together, they describe $`N=2`$ massless gauge multiplets containing, for example, two gluino octets, two photinos, etc., as well as two color-octets of spin-0 gluons, two spin-0 photons, etc. . (These new “scalar gauge fields” may also appear as originating from the fifth and sixth components of higher-dimensional gauge fields $`V^{\widehat{\mu }}`$, in a six-dimensional spacetime .)

Matter (and also Higgs) chiral superfields systematically occur in pairs so as to describe $`N=2`$ “hypermultiplets”. Since quark and lepton fields belong to left-handed electroweak doublets and right-handed singlets, one has to introduce  mirror  particles which belong to right-handed  doublets  and left-handed  singlets . This also necessitates however the introduction of appropriate symmetry breaking mechanisms, allowing one to reduce the full $`N=2`$ supersymmetry down to $`N=1`$.

Mirror quarks and leptons ($`q_M`$ and $`l_M`$) may then acquire large masses through large Yukawa couplings with electroweak Higgs doublets (namely $`H_1`$ and $`H_2`$ in the standard supersymmetry framework). They would have evaded past accelerator searches if heavier than a few hundred GeV’s, but may still show up at LHC.

The breaking of the ($`N=2`$) supersymmetry may be elegantly obtained by demanding periodic and antiperiodic boundary conditions for ordinary $`R`$-even particles and their $`R`$-odd superpartners, respectively – in which case the masses of the (lowest-lying) gravitinos, gluinos and photinos, which fix the energy scale at which supersymmetric particles should start to show up, would be given, in the simplest case and up to radiative correction effects, by $`m_{3/2}=m_{1/2}=\frac{\pi \mathrm{}}{Lc}=\frac{\mathrm{}}{2Rc}.`$ $`L`$ is the size of the extra dimension responsible for supersymmetry breaking. This led us to consider the possibility of relatively “large” extra dimensions, associated with a compactification scale that could then be as “low” as $``$ TeV scale .

It is quite conceivable that the new spin-0 states, as well as mirror lepton and quark fields, may only manifest themselves at the compactification scale. However, whether the extra spin-0 components of gauge fields actually show up or not in the low-energy theory depends on the details of the mechanism that should be responsible for the breaking of the extended supersymmetry (and/or the compactification of the extra space dimensions). Since this mechanism is unknown, one can still discuss the possibility that they appear in the low-energy theory, much below the compactification scale.

The new spin-0 bosons – originally appearing as extra degrees of freedom for spin-1 gauge fields – will then have Yukawa couplings relating quarks and leptons to their mirror partners. They should lead to new decay modes ($`M`$-parity invariant) like

$$q_M(\text{or}l_M)q(\text{or}l)+\text{neutral spin-0 particle}.$$
(4)
Note that $`M`$-parity is equal to $`+1`$ for ordinary and Supersymmetric Standard Model particles (including those described by the doublet Higgs superfields $`H_1`$ and $`H_2`$). It is equal to $`1`$, on the other hand, for mirrors, their superpartners, the new spin-0 particles, the new inos (second octet of gluinos, additional charginos and neutralinos) and finally also for the extra Higgs bosons which appear as the trace of the underlying extended supersymmetry.

The lightest of the new $`M`$-odd particles, that remains at the end of such decay chains, is then expected to be stable.

### C A Light Dark Matter candidate

In the scenario described above, one ends up with a new neutral and stable ($`M`$-odd) particle that can play the rôle of LDM candidate.

This particle could be a scalar (that couples ordinary to mirror particles). In which case, it is expected to annihilate in pairs through mirror fermion exchanges, as discussed in . However, such a scenario is disfavored by the 511 keV line observation (which requires $`v^2`$ suppressed annihilation cross section), unless the associated cross section turns out to be extremely small in the primordial universe (namely $`10^{31}\mathrm{cm}^3\mathrm{s}^1`$, which requires heavy mirror fermions and probably chiral couplings).

The LDM particle may also annihilate through possible couplings to a spin-1 $`U`$ boson, if the extra $`U(1)`$ symmetry generator includes a contribution involving the continuous $`U(1)`$ $`M`$-symmetry generator Such a gauging of an extra $`U(1)`$ not commuting with the $`N=2`$ supersymmetry may however generate anomalies (at least superficially).. In this case, one would obtain an attractive LDM candidate that could explain the low energy signature without being in conflict with a signal at higher energies.

## IV Summary

Supersymmetric theories originating from $`N=2`$ have the pleasant feature of providing two discrete symmetries, that we denote $`R`$ and $`M`$-parities. When they are conserved, one expects that both the LSP (Lightest Particle invariant under $`R`$-parity) and LMP (Lightest Particle of the spectrum invariant under $`M`$-parity) become absolutely stable. Such a scenario would then offer two DM candidates: one would be the $`R`$-odd LSP and the other one, the $`M`$-odd LMP.

Among the new $`M`$-odd particles, one finds spin-0 states. The case where one of them turns out to be the lightest of the $`M`$-odd spectrum appears particularly appealing as it could play the rôle of Light Dark Matter candidate (with a mass of a few MeV-100 MeV). Provided one introduces a new light gauge boson $`U`$, one could evade low energy gamma ray constraints  and also explain the 511 keV emission line (if it is seen to be due to DM annihilations) §§§There is no need to introduce a gauge boson for $``$ 100 MeV but this case appears borderline regarding the 511 keV line, the relic density criterion and the nucleosynthesis constraint..

To avoid an unwanted production in $`Z`$ decays, this spin-0 particle should have no or negligible direct couplings to the $`Z`$. This would restrict potential LMP candidates to e.g. a spin-0 photon, or companion of the weak hypercharge gauge field $`B^\mu `$.

In such a scenario, one would for example end up with the “standard” (lightest) neutralino as the $`R`$-odd LSP and a new neutral spin-0 particle, as $`M`$-odd LMP.

The new $`M`$-odd spin-0 states are expected to couple quarks and leptons to  mirror  partners. However, if the latter are very heavy (a few hundred GeV for example), one does not expect a significant contribution to low energy gamma rays. (The associated annihilation cross section would indeed be dominated by a constant term but this one is not expected to yield any significant contribution as it would be lower than $`10^5`$ pb. ) On the other hand, one expects high energy gamma rays from the annihilations of heavy DM particles for which $`\sigma v^{}>1`$ pb.

This scenario finally evades several astrophysical and particle physics constraints including too large a contribution to the muon and electron $`g2`$. On the other hand, a gauge boson is expected to bring an extra contribution to the $`\nu e`$ elastic scattering cross section at low energy. The latter has been measured by two experiments (LAMPF and LSND). No significant deviations have been found . Therefore, it is necessary to either impose a dissymmetry between the coupling of the new gauge boson to DM and its coupling to standard model particles or, alternatively, to suppress its coupling to neutrinos (e.g. by having the $`q`$ and $`l`$ contribution to the $`U`$ current proportional to the electromagnetic current, as discussed in ).

Of course, there may exist other ways out to achieve a $`v^2`$ suppressed cross section for LDM. We do not discuss them in this letter.

By having two kinds of DM particles, we then fulfill the correct relic density, potentially explain the 511 keV emission line detected by Integral/SPI and simultaneously predict high energy gamma rays.

This would turn out to be particularly interesting either if the smoking gun proposed in  reveals to be in favor of LDM and if high energy diffuse gamma rays were identified as coming from Dark Matter annihilations, or simply if the high energy gamma rays turned out to be associated with a Dark Matter annihilation cross section that is much larger compared to what is needed to achieve the correct relic density.