text:
# 1 Introduction

## 1 Introduction

The confirmation of the zenith–angle–dependent atmospheric neutrino deficit by the Superkamiokande experiment generally has been understood as first significant hint for neutrino masses and thus particle physics beyond the standard model . The other long-standing puzzle of particle physics is the deficit of solar neutrinos . Altogether they constitute the most important milestone in the search for phenomena beyond the Standard Model (SM), indicating the need for oscillations involving all three active neutrino species . The mounting experimental activity in this field promises a bright future for neutrino physics which may prove to be a most valuable source of information on the structure of a more complete theory underlying the standard model of particle physics.

An elegant way to introduce neutrino masses is via the spontaneous breaking of an ungauged lepton number symmetry through a non-zero $`SU(2)U(1)`$ singlet vacuum expectation value (VEV) of a scalar field. This may be implemented in conventional  as well as supersymmetric models . The couplings of the corresponding Goldstone boson, generically called Majoron and denoted by $`J`$, are rather model-dependent . Here we consider the simplest class of Majoron models, where the Majoron-neutrino coupling matrix $`g_{ij}^Mm_{ij}`$ is proportional to the neutrino mass matrix , so that in the mass eigenstate basis diag $`(m_1,m_2,m_3)`$ the Majoron neutrino couplings $`g_i`$ are diagonal, to lowest order approximation 111This proportionality may be avoided in more complex models, such as those in ref. ,

$$g_{ij}^M\delta _{ij}g_i.$$
(1)
This covers a variety of possibilities including both seesaw-type  as well as radiative models .

Limits on this quantity obtained from laboratory experiments searching for Majoron-emitting pion or kaon decays are rather weak, with the exception of double beta decay . On the other hand Majoron emitting neutrino decays affect the expected neutrino luminosity and spectra which are constrained by the observed signal from SN1987A, providing stringent restrictions . While the limits of laboratory experiments on rare decays are given in the weak basis, bounds from processes in supernovae occur in a dense medium and are expressed in the medium eigenstates (see below). In the present work we discuss the correlations of the different limits and their translation into the mass basis, extending the earlier paper ref. . In the next section we derive the expressions for medium and weak eigenstates, following . In section 3 we review the bounds obtained from the supernova SN1987A using various considerations. In contrast to ref.  here we include the study of the effects associated to the Majorana CP violating phases present in theories of massive neutrinos . Moreover we investigate (Section 4) the recent bounds from neutrinoless double beta decay as well as those that could be attained in future experiments such as GENIUS . The resulting exclusion plots are discussed in section 5 in the mass basis.

## 2 Neutrino mixing in three bases

For neutrinos propagating through a medium one has to deal with three kinds of eigenstates: Flavor eigenstates $`\nu _\alpha `$, mass eigenstates $`\nu _i^{(h_i)}`$ with masses $`m_i`$, and, depending on the environment, medium eigenstates $`\stackrel{~}{\nu }_i^{(h_i)}`$. The flavor eigenstates are defined as

$$\nu _\alpha =\underset{i}{}U_{\alpha i}\nu _i^{(h_i)}$$
(2)
and the medium eigenstates are

$$\stackrel{~}{\nu }_i^{(h_i)}=\underset{j}{}\stackrel{~}{U}_{ij}\nu _j^{(h_j)}.$$
(3)
Here the superscript $`(h_i)=\pm 1`$ refers to the helicity of the state. In the general case it is impossible to diagonalize simultaneously the mass and potential terms. Thus one has to solve the field equations in detail. In a two-component field formalism, where a left-handed four-component field $`\nu `$ expressed in the chiral representation of the $`\gamma `$ matrices  is related to the corresponding two-component field $`\varphi `$ by $`\nu _L^T=(\varphi ^T,0)`$  222The notation here coincides to the one of  up to a factor of $`i`$. , the Lagrangian can be written in the mass basis as

$`_{\mathrm{tot}}`$ $`=`$ $`_0+_{\mathrm{med}}+_{\mathrm{int}}`$ (4)
$`=`$ $`{\displaystyle \underset{i}{}}\varphi _i^{}(i_0i\stackrel{}{\sigma }\times \stackrel{}{})\varphi _i{\displaystyle \frac{m_i}{2}}(\varphi _i^Ti\sigma _2\varphi _i\varphi _i^{}i\sigma _2\varphi ^{})+{\displaystyle \underset{i,j}{}}\varphi _i^{}V_{ij}\varphi _j`$
$`J{\displaystyle \underset{ij}{}}g_{ij}^M(\varphi _i^T\sigma _2\varphi _j+\varphi _i^{}\sigma _2\varphi _j^{})`$
where $`\sigma _2`$ and $`\stackrel{}{\sigma }`$ denote Pauli matrices. Here the free Langrangian $`_0`$ and $`_{\mathrm{med}}`$ describe the propagation in vacuo and the effects of matter described by the potential matrix $`V_{ij}`$, respectively, whereas $`_{\mathrm{int}}`$ takes into account the presence of neutrino-Majoron interactions which may lead to decays. One has now to consider the decays $`\stackrel{~}{\nu }_i^{h_i}(p_i)\stackrel{~}{\nu }_j^{h_j}(p_j)+J(q)`$, where $`\stackrel{~}{\nu }_i^{h_i}(p_i)`$ and $`\stackrel{~}{\nu }_j^{h_j}(p_j)`$ are energy–eigenstate Majorana neutrinos that propagate in matter with four-momenta $`p_i=(E_i^{h_i},\stackrel{}{p}_i)`$ and $`p_j=(E_j^{h_j},\stackrel{}{p}_j)`$, and helicity $`h_i`$ and $`h_f`$ respectively. In order to obtain these energy–eigenstates one has to take $`_0+_{\mathrm{med}}`$ and calculate the resulting field equations.

$$(i_0i\stackrel{}{\sigma }\times \stackrel{}{})\varphi _i(x)+m_ii\sigma _2\varphi _i^{}\underset{j=1}{\overset{N}{}}V_{ij}\varphi _j(x)=0.$$
(5)
One solves these field equations by expanding the fields $`\varphi _i(x)`$ as superpositions of plane-wave spinors with definite helicity, ,

$`\varphi _i(x)={\displaystyle \frac{d\stackrel{}{p}}{\sqrt{(2\pi )^3}}e^{i\stackrel{}{p}\times \stackrel{}{x}}\left\{[P_\alpha ^i(\stackrel{}{p},t)+N_\alpha ^i(\stackrel{}{p},t)]\alpha (\stackrel{}{p})+[P_\beta ^i(\stackrel{}{p},t)+N_\beta ^i(\stackrel{}{p},t)]\beta (\stackrel{}{p})\right\}}`$ (6)
where $`\alpha (\stackrel{}{p})`$ and $`\beta (\stackrel{}{k})`$ are helicity eigenstates and $`P^i`$ and $`N^i`$ denote positive and negative frequency components of the field under consideration. One should now substitute this expression in the equations (5), whose diagonalization would give rise to the desired eigenstates. It can be shown, though, that for relativistic neutrinos the positive-frequency components decouple from the negative-frequency ones and the energy eigenstates obtained in this way result to be the same as those obtained from the diagonalization of the usual MSW equation , which can be stated as

$$i_t\nu _i^{(h)}=(H_{ij}^{\mathrm{rel}}+U_{i\alpha }V_{\alpha \beta }U_{\beta j}^{})\nu _j^{(h)}.$$
(7)
Here $`H_{ij}^{\mathrm{rel}}(p+m_i^2/(2p))\delta _{ij}`$ and $`V_{\alpha \beta }`$ is the potential matrix in the weak basis,

$$V_{\alpha \beta }=\left(\begin{array}{ccc}V_C+V_N& 0& 0\\ 0& V_N& 0\\ 0& 0& V_N\end{array}\right).$$
(8)
The potentials induced by the charged and neutral currents are $`V_C=\sqrt{2}hG_Fn_B(Y_e+Y_{\nu _e})`$ and $`V_N=\sqrt{2}hG_Fn_B\left(\frac{1}{2}Y_N+Y_{\nu _e}\right)`$, where $`Y_i=(n_in_{\overline{i}})/n_B`$ and $`n_B`$ is the baryon density. Diagonalizing $`H^{\mathrm{rel}}+UVU^{}`$ yields the medium eigenstates $`\stackrel{~}{\nu }_i^{(h)}=\stackrel{~}{U}_{ij}^{(h)}\nu _j^{(h)}`$.

In the three-flavor neutrino case the mixing matrix $`U`$ can be parametrized as $`U=U_{23}U_{13}U_{12}U_0`$, where the matrices $`U_{ij}=U_{ij}(\theta _{ij})`$ perform the rotation in the $`ij`$-plane by the angle $`\theta _{ij}`$ and $`U_0`$ includes possible CP-violation effects . In the following we will assumme $`\theta _{13}=0`$, motivated both by detailed fits of the present solar and atmospheric neutrino anomalies  as well as by the reactor results of the Chooz experiment . This simplifies the mixing matrix to $`\nu _\alpha =U_{\alpha i}\nu _i=U_{23}U_{12}U_0\nu _i`$  and allows us to set $`\theta _{12}=\theta _{}`$ and $`\theta _{23}=\theta _{\mathrm{atm}}`$. Now for light neutrinos near the neutrinospheres in supernovae the condition $`|V_{\alpha \alpha }|m_i^2/(2p)`$ holds and, since in the weak basis the potential is diagonal, the medium states can be identified with the weak ones up to an arbitrary rotation in the $`\nu _\mu \nu _\tau `$ subspace. In order to simplify the expressions we exploit this freedom by choosing this arbitrary rotation angle to coincide with $`\theta _{23}`$, see Table 1. This allows us to identify the coupling matrix in the medium basis with the one in the weak basis up to the rotation

$`\stackrel{~}{g}_{ij}=U(\theta _{23})g_{\alpha \beta }^WU^T(\theta _{23})g_{\alpha ^{}\beta ^{}}.`$ (9)
Taking now into account that $`g^M=U^Tg^WU`$ and substituting the explicit expressions for the $`U_{\alpha i}`$ matrices relating mass and weak eigenstates, one gets the following expression, $`\stackrel{~}{g}_{ij}=U_{12}U_0^{}g_{ij}^MU_0^{}U_{12}^T`$, or explicitly

$`\stackrel{~}{g}=\left(\begin{array}{ccc}g_{ee}& g_{e\mu ^{}}& g_{e\tau ^{}}\\ g_{\mu ^{}e}& g_{\mu ^{}\mu ^{}}& g_{\mu ^{}\tau ^{}}\\ g_{\tau ^{}e}& g_{\tau ^{}\mu ^{}}& g_{\tau ^{}\tau ^{}}\end{array}\right)=\left(\begin{array}{ccc}g_1\mathrm{cos}^2\theta _{}+g_2\mathrm{sin}^2\theta _{}e^{2i\delta }& \frac{1}{2}(g_1+g_2e^{2i\delta })\mathrm{sin}2\theta _{}& 0\\ \frac{1}{2}(g_1+g_2e^{2i\delta })\mathrm{sin}2\theta _{}& g_1\mathrm{sin}^2\theta _{}+g_2\mathrm{cos}^2\theta _{}e^{2i\delta }& 0\\ 0& 0& g_3\end{array}\right).`$ (16)
This choice of the rotation in the $`\nu _\mu \nu _\tau `$ subspace leads to a relation between medium and mass eigenstates characterized only by the solar angle $`\theta _{}`$ and by the Majorana CP violating phase $`\delta `$ . Using the definitions $`\mathrm{\Delta }m_{12}^2=\mathrm{\Delta }m_{}^2`$ and $`\mathrm{\Delta }m_{23}^2=\mathrm{\Delta }m_{atm}^2`$ together with the assumptions in 1 one can easily translate the bounds obtained in the weak or medium basis into the mass basis and, in addition, express them in terms of only two independent parameters, for instance $`(m_1,g_1)`$ via

$`g_2=g_1\sqrt{1+{\displaystyle \frac{\mathrm{\Delta }m_{}^2}{m_1^2}}},g_3=g_1\sqrt{1+{\displaystyle \frac{\mathrm{\Delta }m_{}^2+\mathrm{\Delta }m_{atm}^2}{m_1^2}}}.`$ (17)
## 3 Supernova bounds

There is a variety of different arguments based on supernova physics which lead to restrictions on neutrino properties. Processes involving Majoron-neutrino couplings may prevent a successful explosion as well as substantially affect the observed neutrino spectra. A crucial feature to notice is that the effective mass induced by the interactions of neutrinos with background matter breaks the proportionality between the neutrino mass matrix and the neutrino-Majoron coupling matrix $`g_{ij}^M`$. This follows from the fact that the thermal background in the supernova environment consists only of particles of the 1st generation, thus distinguishing the electron flavour from the others. We now describe three different arguments used  in order to restrict the relevant parameters.

### 3.1 Constraints from Neutrino spectra

The idea behind this bound is that Majoron–induced transitions between the neutrino flavors could change the energy spectra of the single flavors. At the typical temperatures of the SN core $`\nu _{\mu ,\tau }`$ only interact with the medium via neutral currents giving rise to a smaller cross section than that corresponding to the electron (anti)neutrinos, which feel both neutral and charged currents. Since the opacity of the heavier $`\nu _{\mu ,\tau }`$ flavors is smaller than for the $`\nu _e`$, their energy-exchanging reactions freeze out in the denser region of the protoneutron star, leading to lower spectral temperatures of $`\nu _e`$ compared to $`\nu _{\mu ,\tau }`$. This expected spectrum can be distorted due to the decays $`\stackrel{~}{\nu }_i^\pm \stackrel{~}{\nu }_j^{}+J`$. Besides the effects of such decays one has to keep in mind the possible oscillation which neutrinos could undergo along their journey to the Earth. In order to consider both aspects we have defined the effective survival probability as

$`N=N_{decay}\times N_{osc},`$ (18)
where $`N_{decay}`$ stands for the survival probability of a $`\stackrel{~}{\nu }_i^\pm `$ emitted from its energy sphere and can be computed as

$`N_{decay}(\stackrel{~}{\nu }_i^\pm )\mathrm{exp}\left\{{\displaystyle _{R_{E,\stackrel{~}{\nu }_i^\pm }}^{\mathrm{}}}𝑑r{\displaystyle \underset{j}{}}\mathrm{\Gamma }(\stackrel{~}{\nu }_i^\pm \stackrel{~}{\nu }_j^{}+J)\right\}.`$ (19)
Within our relativistic approximation, the helicity-flipping neutrino decays rate are given by

$`\mathrm{\Gamma }(\stackrel{~}{\nu }_i^\pm \stackrel{~}{\nu }_j^{}+J)={\displaystyle \frac{|\stackrel{~}{g}_{ij}|^2}{16\pi }}(V_iV_f).`$ (20)
Coming to the oscillation term, the corresponding neutrino survival probability $`N_{osc}`$, will depend on the neutrino mixing angle and squared mass difference. We will analyse separately the different solutions of the solar neutrino problem namely small-angle MSW (SMA-MSW), large-angle MSW (LMA-MSW), LOW-MSW and the just-so case. Details about the present status and required parameters of the various solutions can be found in the global analysis of neutrino data presented in ref. . Such a study favours a rather small value for the angle $`\theta _{13}`$, mainly because of data from reactors . In the first three cases neutrinos will propagate through the supernova environment adiabatically. Therefore they will emerge as energy eigenstates, which in vacuum coincide with the mass eigenstates, without any oscillation occuring on the way from the SN to Earth. If one takes into account that neutrinos have to traverse a distance, $`d`$, of matter in the Earth to reach the detectors, Kamiokande and IMB, one gets the following expression for their survival probability ,

$`N_{osc}=1\{\mathrm{sin}^2\theta _{}\mathrm{sin}2\theta _m\mathrm{sin}(2\theta _{}2\theta _m)\mathrm{sin}^2\left({\displaystyle \frac{\pi d}{l_m}}\right)\},`$ (21)
where $`l_m`$ and $`\theta _m`$ denote the oscillation length and the mixing angle in matter, respectively. As has been previously noted for the simplest choice $`\theta _{13}=0`$ one has that, besides the fact that $`\nu _\mu `$ and $`\nu _\tau `$ behave the same way in the supernova, the conversion $`\nu _e\nu _\mu ^{}`$ will be the only oscillation involving electron (anti)neutrinos, allowing us to set the angle which characterizes their mixing to $`\theta _{}`$.

In the vacuum solution case the neutrinos emerge essentially as flavor eigenstates which then oscillate on their way to Earth. Therefore one has

$`N_{osc}=1{\displaystyle \frac{1}{2}}\mathrm{sin}^22\theta _{}.`$ (22)
In order to get information on the coupling constants we will conservatively require that at least half of the initial electron antineutrinos emerging from the SN1987A survive, $`N>0.5`$, accounting for the rough agreement between the expected and the detected SN1987A signals. In order to analyse the implications of this restriction one must generalize the simplest argument used in  since neutrinos may loose energy as a result of majoron decays.

This allows us to get some limits on the coupling parameter of the order of $`g_1(g_2)\mathrm{few}\times 10^4`$ from the first three solutions to the solar neutrino problem. For the case of vacuum oscillations, though, the solution is already disfavored by the SN1987A data even in the absence of neutrino decays . Though they may narrow it down considerably, the above arguments do not totally close the allowed window of neutrino-Majoron couplings, neither for the SMA, LMA nor LOW solutions, even for a supernova in our milky way.

### 3.2 Constraints from Majoron luminosity

This bound is based on the observation that neutrino decays into Majorons could supress the energy release contained in the neutrino signal. Under the assumption of small $`\nu _e\nu _x`$ mixing the neutrino signal observed in SN 1987A is in good agreement with numerical computations of the total binding energy released in a supernova explosion. An analysis of decay and scattering processes involved yields the exclusion region

$$3\times 10^7<|\stackrel{~}{g}_{ij}|<2\times 10^5.$$
(23)
For $`|\stackrel{~}{g}_{ij}|`$ values smaller than $`3\times 10^7`$ the Majoron neutrino coupling becomes too small to induce any effect. On the other hand for $`|\stackrel{~}{g}_{ij}|>2\times 10^5`$ Majorons get trapped in the core and do not contribute to the energy release.

Another point to observe is that CP violating phases affect these limits. This follows from the appearance of the phase $`\delta `$ in explicit form of the Majoron neutrino coupling constants given in eq. (16). In order to eliminate such an explict CP phase dependence when translating the limit on $`|\stackrel{~}{g}_{ij}|`$ into the mass basis we have analyzed for each term of eq. (16) the excluded region for different values of $`\delta `$ and subsequently considered the intersection of the resulting excluded regions. This conservative procedure allows us to rule out part of the parameter space irrespective of the value of the CP phase.

As an example we illustrate in Figs. 1 and 2 the regions excluded for the LMA MSW solution to the solar neutrino problem. The luminosity bound can be described in two steps. In the first one we take one $`\stackrel{~}{g}_{ij}`$ from eq. (16) and by means of eq. (17) we write it in terms of $`g_1m_1`$. This way we obtain an expression for the energy loss which depends explicitly upon the CP phase $`\delta `$. Now, by varying that CP phase the bound given in eq. (23) is translated into different ruled out regions. We show in Fig. 1 the resulting bound on $`|g_{ee}|`$ assuming two extreme cases, $`\delta =0`$ (solid lines), and $`\delta =\pi /2`$ (dotted lines). Notice that for the latter case the bound disappears because of a cancellation between the two terms in $`|g_{ee}|`$. In order to remove the phase we therefore consider the intersection as the most conservative choice.

Now turn to the implications of the luminosity bound to the other components of the Majoron neutrino coupling matrix elements. Once we have obtained those intersecting regions for each $`\stackrel{~}{g}_{ij}`$ we simply take the union of them, giving rise to a final highly non-trivial exclusion region, as can be seen in Fig. 2.

It is important to notice that the shape of such regions is characterized by the values of the square root of $`\mathrm{\Delta }m_{}^2`$ and $`\mathrm{\Delta }m_{atm}^2`$. Let us first consider $`\stackrel{~}{g}_{ij}`$ with $`i,j3`$. In this case only $`\mathrm{\Delta }m_{}^2`$ appears in eq. (17), so that for $`m_1\mathrm{\Delta }m_{}^2`$ one has $`g_2g_1`$ giving rise to a vertical line with no dependence on $`m_1`$, as noted in the figures. In contrast, for $`m_1\mathrm{\Delta }m_{}^2`$ one has $`g_2g_1\frac{\sqrt{\mathrm{\Delta }m_{}^2}}{m_1}`$ with a explicit dependence on $`m_1`$ which strengthens the bound for lower $`m_1`$ values. Let us now consider the limit coming from $`g_3`$. In this case the characteristic mass scale is always given by $`\mathrm{\Delta }m_{atm}^2`$, eq. (17), irrespective of the particular solutions to solar neutrino problem that one may wish to consider. As a result, for the LOW and (quasi)–vacuum cases the difference between $`\mathrm{\Delta }m_{atm}^2`$ and the solar mass scale $`\mathrm{\Delta }m_{}^2`$ is so large that two branches appear. This explains the two branches observed in figures 5 and 6 corresponding to the LOW and (quasi)–vacuum solutions, respectively (section 4).

Concerning the SMA solution in fig. 4 the main changes arise from the bounds on $`|g_{ee}|`$ and $`|g_{e\mu ^{}}|`$. In the expression of $`g_{ee}`$ in eq. (10) the contributions of $`g_1`$ and $`g_2`$ may cancel for a phase $`\delta \pi /2`$, see fig. 1, as long as these contributions are of comparable magnitude. To fulfill this requirement, a smaller admixture of $`g_2`$ in $`g_{ee}`$, as happens for the SMA solution, has to be compensated by a larger ratio of $`g_2/g_1`$, corresponding to smaller masses $`m_1`$. This is producing a small hollow in the bound at $`m_15\times 10^6`$.

The bound on $`|g_{e\mu ^{}}|`$ is responsible for the sharp peak at the right edge of the excluded region, as can be seen in fig. 2, where all Majoron luminosity bounds are shown explicitly. Here the conservative upper bound on $`|g_{e\mu ^{}}|`$ is obtained, when $`\delta =0`$, corresponding to cancellation of the $`g_1`$ and $`g_2`$ contributions in the large $`m_1^2`$ asymptotics, see eq. (10). Correspondingly the right border of the excluded region is obtained for $`\delta =\pi /2`$, where $`|g_{e\mu ^{}}|`$ becomes constant for large $`m_1`$. The intersection gives rise to the peak. For smaller values of the mixing, as obtained for the SMA solution, the expression for $`|g_{e\mu ^{}}|`$ in eq. (10) is fulfilled for corresponding larger values of $`g_1`$. This shifts the exclusion region to the right, making it more visible in fig. 4.

Before concluding this section we mention the constraints on Majoron-neutrino coupling parameters which arise from the collapsing phase. The idea behind this bound is that a change in the trapped electron fraction could prevent a successful explosion process. At the end of their life massive stars become unstable and, when the iron core reaches the Chandrasekar limit, they implode. Once the nuclear density is reached, a shock wave forms at the core and propagates outwards, turning the implosion into an explosion. The strength and propagation of this shock is sensitive to the trapped electron fraction $`Y_{L_e}=Y_e+Y_{\nu _e}`$, which can be erased by neutrino decays $`\nu _e\overline{\nu }_e+J`$. Requiring $`Y_L(t_{bounce})>0.375`$ leads to a limit of

$$|g_{ee}|<2\times 10^6.$$
(24)
However to the extent that current supernova models do not fully account for the explosion mechanism this limit should be taken only as indicative for the moment.

## 4 Neutrinoless double beta decay

The only laboratory experiment, which is competitive with the supernova bounds, is neutrinoless double beta decay. This decay corresponds to two single beta decays occuring in one nucleus and converts a nucleus (Z,A) into a nucleus (Z+2,A). Limits on the Majoron emitting mode

$${}_{Z}{}^{A}X_{Z+2}^AX+2e^{}+\varphi $$
(25)
are given by two types of experiments. In geochemical experiments the half–life limit is derived from relative abundances of nuclear isotopes found in the earth :

$$|g_{ee}|<3\times 10^5.$$
(26)
However, half–life determinations vary by more than a factor of three.

The best direct laboratory limit (less stringent but more reliable) from the Heidelberg-Moscow experiment  is based on a likelyhood fit to the continuous electron spectrum:

$$|g_{ee}|<8\times 10^5.$$
(27)
Future projects such as GENIUS  and EXO  aim at considerable improvements in the sensitivity. A very rough estimation of the sensitivity of GENIUS 1t is based on the background simulation in , where a background improvement in the interesting energy range of a factor $`1000`$ has been obtained. Since in the Heidelberg-Moscow experiment the Majoron-neutrino coupling bounds are dominated by the systematical error of the background simulation, a considerable reduction of the background by a factor of $`B`$ will reduce the limit on the Majoron-emitting double beta decay half life by $`\sqrt{B}`$ and the coupling constant limit accordingly by $`^4\sqrt{B}`$. This implies a reach of sensitivity down to $`|g_{ee}|10^5`$, which could bridge most of the gap existing between the more reliable limits derived from supernovae.

## 5 Discussion and conclusions

In figures 3 to 6 we present the limits on Majoron-neutrino couplings in terms of $`m_1g_1`$ corresponding to the various solutions of the solar neutrino problem. In Fig. 7 we display the results for the LMA solution also in terms of the equivalent $`m_2g_2`$ variables. This representation has been selected for convenience and generality. By further specifying the underlying model for lepton number violation one can re-express our results in terms of the lepton number breaking VEV, which will provide also useful information for model-builders.

Regions which are excluded by supernova arguments are denoted by the rhombical pattern (obtained from Majoron luminosity) and by the vertical lines (obtained from the neutrino spectra). Also shown are the regions excluded from neutrinoless double beta decay (horizontal lines). The excluded region from Majoron luminosity is a superposition of the bounds on $`\stackrel{~}{g}_{11}`$, $`\stackrel{~}{g}_{12}`$, $`\stackrel{~}{g}_{22}`$, and $`\stackrel{~}{g}_{33}`$, where always the most conservative limits for various CP Majorana phases have been used. Due to the expressions for the helicity-flipping neutrino decays the bound obtained from the neutrino spectra turns out to be independent of the CP phase. For neutrinoless double beta decay one can have a cancellation of the coupling constants $`g_1`$ and $`g_2`$. The expected sensitivity of the GENIUS experiment is shown as a dashed line. It is easy to see that GENIUS could be able to bridge almost the whole gap between the different supernova constraints. Also an upper bound $`m_1<2.3`$ eV from tritium beta decay is displayed.

The limits obtained in this paper apply to the simplest class of models where neutrino masses arise from the spontaneous violation of lepton number. Such Majoron models which cover a wide and attractive class including both models where the smallness of neutrino masses follow from a seesaw scheme, as well as those where it arises from the radiative corrections.

Both neutrinoless double beta decay as well as supernova physics arguments provide stringent limits on Majoron-neutrino interactions. In the present work we have discussed these limits and their translation into the mass basis. Generalizing previous papers  we have now taken into account the effect of CP violating phases, which play a crucial role in the neutrinoless double beta decay limits. Depending on the solution of the solar neutrino problem and the absolute mass scale in the neutrino sector the constraint from the supernova energy release (Majoron luminosity argument) exludes Majoron-neutrino couplings in the wide range of $`10^710^5`$. Upper bounds have been obtained from neutrinoless double beta decay and the SN87A neutrino spectra. An estimate of the potential of the future double beta projects such as GENIUS suggests the possibility to bridge almost the whole gap separating the excluded areas and either to establish Majorons with couplings around a few $`10^5`$ or to restrict neutrino-Majoron couplings down to $`10^7`$.

Last, but not least, let us mention that the propagation of neutrinos produced in the solar interior follows essentially the MSW picture, while any possible effect of decays would happen in vacuo through a non-diagonal neutrino-majoron coupling which is absent in the simplest models considered here . Even in more complex models  where such non-diagonal neutrino-majoron couplings exist in vacuo, one can see that for such small values of the neutrino-majoron coupling strengths indicated by supernova and neutrinoless double beta decay, it is rather unlikely that they can play any role whatsoever in the solar neutrino problem .

## Acknowledgement

We thank M. Hirsch for useful discussions. This work was supported by Spanish DGICYT under grant PB98-0693, by the European Commission TMR contract HPRN-CT-2000-00148 and by the European Science Foundation network grant N. 86. H. P. was supported by TMR contract ERBFMRX-CT96-0090 and DOE grant no. DE-FG05-85ER40226 and R.T. by a grant from the Generalitat Valenciana.


format: text

text:
# Anharmonic parametric excitation in optical lattices

## I Introduction

The phenomenon of parametric excitation of the motion of cold trapped atoms has recently been the subject of several theoretical and experimental investigations . The excitation caused by resonant amplitude noise has been proposed as one of the major sources of heating in far–off–resonance optical traps (FORTs), where the heating due to spontaneous scattering forces is strongly reduced. In particular, the effect of resonant excitation is expected to be particularly important in optical lattices, which usually provide a very strong confinement to the atoms, resulting in a large vibrational frequency and in a correspondingly large transfer of energy from the noise field to the atoms.

Nevertheless, parametric excitation is not only a source of heating, but it also represents a very useful tool to characterize the spring constant of a FORT or in general of a trap for cold particles, and to study the dynamics of the trapped gas. Indeed, the trap frequencies can be measured by intentionally exciting the trap vibrational modes with a small modulation of the amplitude of the trapping potential, which results in heating  or losses for the trapped atoms when the modulation frequency is tuned to twice the oscillation frequency. This procedure usually yields frequencies that satisfactorily agree with calculated values, and are indeed expected to be accurate for the atoms at the bottom of the trapping potential. From the measured trap frequencies is then possible to estimate quantities such as the trap depth and the number and phase space densities of trapped atoms. We note that this kind of measurement is particularly important in optical lattices, since the spatial resolution of standard imaging techniques is usually not enough to estimate the atomic density from a measurement of the volume of a single lattice site.

Recently, 1D lattices have proved to be the proper environment to study collisional processes in large and dense samples of cold atoms, using a trapping potential independent for the magnetic state of the atoms. In this systems, the parametric excitation of the energetic vibrational mode along the lattice provides an efficient way to investigate the cross-dimensional rethermalization dynamics mediated by elastic collisions .

Most theoretical studies of parametric excitation rely on a classical  or quantum  harmonic approximation of the confining potential. Under certain circumstances these expressions show quite good agreement with experimental results . However, general features of the optical lattice could be lost in these approaches. For example, a sinusoidal potential exhibits an energy band structure and a spread of transition energies, while harmonic oscillators have just a discrete equidistant spectrum. Thus, we might expect that the excitation process may happen at several frequencies, and with a non-negligible bandwidth. Such anharmonic effects can be important whenever the atoms are occupying a relatively large fraction of the lattice energy levels. The purpose of this paper is to give a simple description of parametric excitation in a sinusoidal 1D lattice. In the next Section, we briefly discuss general features of the stationary states on such a lattice. Then, we summarize the harmonic description given in Ref. and extend it to the anharmonic case. By a numerical evaluation of transition rates, we make a temporal description of parametric excitation which is compared with experimental results. We discuss about the relevance of broadening of the spectral lines to understand the excitation process in this kind of systems. Some conclusions are given in the last Section.

## II Stationary states of a sinusoidal optical lattice.

The Hamiltonian for an atom in a red detuned FORT is

$$H=\frac{P^2}{2M}+V_{\mathrm{eff}}(\stackrel{}{x}),$$
(1)
with

$$V_{\mathrm{eff}}(\stackrel{}{x})=\frac{1}{4}\alpha |(\stackrel{}{x})|^2,$$
(2)
where $`\alpha `$ is the effective atomic polarizations and $`(x)`$ is the radiation field amplitude. For the axial motion in a sinusoidal 1D lattice we can take

$`H_{ax}`$ $`=`$ $`{\displaystyle \frac{P_z^2}{2M}}+V_0\mathrm{cos}^2(kz)`$ (3)
$`=`$ $`{\displaystyle \frac{P_z^2}{2M}}+{\displaystyle \frac{V_0}{2}}(1+\mathrm{cos}(2kz))).`$ (4)
The corresponding stationary Schrödinger equation

$$\frac{\mathrm{}^2}{2M}\frac{d^2\mathrm{\Phi }}{dz^2}+\frac{V_0}{2}\left(1+\mathrm{cos}(2kz)\right)\mathrm{\Phi }=E\mathrm{\Phi }$$
(5)
can be written in canonical Mathieu’s form

$$\frac{d^2\mathrm{\Phi }}{du^2}+(a2q\mathrm{cos}2u)\mathrm{\Phi }=0$$
(6)
with

$$a=\left(E\frac{V_0}{2}\right)\left(\frac{2M}{\mathrm{}^2k^2}\right)2q=\frac{V_0}{2}\left(\frac{2M}{\mathrm{}^2k^2}\right).$$
(7)
It is well known that there exists countably infinite sets of characteristic values $`\{a_r\}`$ and $`\{b_r\}`$ which respectively yield even and odd periodic solutions of Mathieu equation. These values also separate regions of stability. In particular, for $`q0`$ the band structure of the sinusoidal lattice corresponds to energy eigenvalues between $`a_r`$ and $`b_{r+1}`$ . The unstable regions are between $`b_r`$ and $`a_r`$. For $`q>>1`$, there is an analytical expression for the band width :

$$b_{r+1}a_r2^{4r+5}\sqrt{2/\pi }q^{\frac{1}{2}r+\frac{3}{4}}e^{4\sqrt{q}}/r!.$$
(8)
The quantities defined above can be expressed in terms of a frequency $`\omega _0`$ defined in the harmonic approximation of the potential

$$\frac{1}{2}M\omega _0^2=\frac{V_0}{2}\frac{(2k)^2}{2!},$$
(9)
thus obtaining

$$a=\left(E\frac{V_0}{2}\right)\left(\frac{4V_0}{\mathrm{}^2\omega _0^2}\right);q=\left(\frac{V_0}{\mathrm{}\omega _0}\right)^2.$$
(10)
Thus, the width of the r-band can be estimated using Eq. (8) whenever the condition $`(V_0/\mathrm{}\omega _0)^2>>1`$ is satisfied. In the experiment we shall be working with a 1D optical lattice having $`V_010.5\mathrm{}\omega _0`$. While the lowest band $`r=0`$ has a negligible width $`10^{18}\mathrm{}\omega _0`$, the band widths for highest lying levels $`r=10,11,12`$,and $`13`$ would respectively be $`0.0065`$ ,$`0.1036`$, $`1.52`$, and $`20.56`$ in units of $`\mathrm{}\omega _0`$.

In order to determine the energy spectrum, a variational calculation can be performed. We considered a harmonic oscillator basis set centered in a given site of the lattice, and with frequency $`\omega _0`$. The diagonalization of the Hamiltonian matrix associated to (4) using 40 basis functions gives the eigenvalues $`E_n<V_0`$ shown in Table I for $`V_0=10.5\mathrm{}\omega _0`$. According to the results of last paragraph, the eigenvalues $`12`$ and $`13`$ belong to the same band while the band width for lower levels is smaller than $`0.11\mathrm{}\omega _0`$

## III Parametric excitation

As already mentioned, parametric excitation of the trapped atoms consists in applying a small modulation to the intensity of the trapping light,

$$H=\frac{P^2}{2M}+V_{\mathrm{eff}}[1+ϵ(t)].$$
(11)
Within first order perturbation theory, this additional field induces transitions between the stationary states $`n`$ and $`m`$ with an averaged rate

$`R_{mn}`$ $`=`$ $`{\displaystyle \frac{1}{T}}\left|{\displaystyle \frac{i}{\mathrm{}}}{\displaystyle _0^T}𝑑tT(m,n)ϵ(t)e^{i\omega _{mn}t}\right|^2`$ (12)
$`=`$ $`{\displaystyle \frac{\pi }{2\mathrm{}^2}}|T(m,n)|^2S(\omega _{mn});\omega _{mn}={\displaystyle \frac{E_mE_n}{\mathrm{}}}`$ (13)
where

$`T(m,n)`$ $`=`$ $`m|V_{\mathrm{eff}}|n`$ (14)
$`=`$ $`E_n\delta _{nm}{\displaystyle \frac{1}{2M}}m|\widehat{P}^2|n`$ (15)
is the matrix element of the space part of the perturbation and

$$S(\omega )=\frac{2}{\pi }_0^T𝑑\tau \mathrm{cos}\omega \tau ϵ(t)ϵ(t+\tau )$$
(16)
is the one-sided power spectrum of the two-time correlation function associated to the excitation field amplitude.

If the confining potential is approximated by a harmonic well, the transition rates different from zero are

$`R_{nn}`$ $`=`$ $`{\displaystyle \frac{\pi \omega _0^2}{16}}S(0)(2n+1)`$ (17)
$`R_{n\pm 2n}`$ $`=`$ $`{\displaystyle \frac{\pi \omega _0^2}{16}}S(2\omega _0)(n+1\pm 1)(n\pm 1)`$ (18)
The latter equation was used in  to obtain a simple expression for the heating rate,

$$\dot{E}=\frac{\pi }{2}\omega _0^2S(2\omega _0)E,$$
(19)
showing its exponential character. The dependence on $`2\omega _0`$ is characteristic of the parametric nature of the excitation process. The fact that $`\mathrm{}`$ is not present is consistent with the applicability of Eq. (19) in the classical regime.

Classically, parametric harmonic oscillators exhibit resonances not just at $`2\omega _0`$ but also at $`2\omega _0/n`$ with $`n`$ any natural number . In fact, the resonances corresponding to $`n`$=2, i.e. at an excitation frequency $`\omega _0`$, have been observed in optical lattices . A quantum description of parametric harmonic excitation also predicts resonances at the same frequencies via $`n`$-th order perturbation theory . In particular, the presence of the resonance at $`\omega _0`$ can be justified with the following argument. According to the standard procedure, the second order correction to the transition amplitude between states $`|n`$ and $`|m`$ is given by

$`R_{mn}^{(2)}=n|𝒰^{(2)}(t_0,t)|m={\displaystyle \underset{k}{}}\left({\displaystyle \frac{i}{\mathrm{}}}\right)^2T(n,k)T(k,m)`$ (20)
$`{\displaystyle _{t_0}^t}𝑑t^{}e^{i\omega _{nk}t^{}}ϵ(t^{}){\displaystyle _{t_0}^t}𝑑t^{\prime \prime }e^{i\omega _{km}t^{\prime \prime }}ϵ(t^{\prime \prime })`$ (21)
with $`𝒰^{(2)}(t_0,t)`$ the second order correction to the evolution operator $`𝒰`$. Therefore, the transition may be described as a two step procedure $`|m|k|n`$. For harmonic parametric excitation the matrix element of the space part of the perturbation differs from zero just for transitions $`|n|n`$ and $`|n\pm 2|n`$. Consider a transition in Eq. (21) involving a ”first” step in which the state does not change $`|n|n`$ and a ”second” step for which $`|n\pm 2|n`$. Then resonance phenomena occur when the total energy of the two excitations, 2$`\mathrm{}\mathrm{\Omega }`$, coincides with that of the second step transition, $`i.`$ $`e.`$ for an excitation frequency $`\mathrm{\Omega }=\omega _o`$.

These ideas can be directly extended to anharmonic potentials: the corresponding transition probability rates $`R(n,m)`$ would be determined by the transition matrix $`T(n,m)`$, by the transition frequencies $`\omega _{nm}`$ and by the time dependence of the excitation $`ϵ(t)`$. In general, anharmonic transition matrix elements $`T(n,m)`$ will be different from zero for a wider set of pairs $`(n,m)`$. Besides, the transition energies will not be unique so that the excitation process is $`not`$ determined by the excitation power spectrum at a single given frequency $`2\omega _0`$ and its subharmonics $`2\omega _o/n`$. As an example the transition energies for the specific potential considered in this work are reported in Table I. Therefore, within the model Hamiltonian of Eq. (11), resonance effects can occur for several frequencies that may alter the shape of the population distribution within the trap. However, in general these resonant excitations will not be associated with the escape of trapped atoms.

Here we are interested in a 1D lattice; the direct extension of the formalism mentioned above requires the evaluation of the matrix elements $`T(n,m)`$ among the different Mathieu states that conform a band. This involves integrals which, to our knowledge, lack an analytical expression and require numerical evaluation. As an alternative, we consider functions which variationally approximate the Mathieu functions. They are the eigenstates of the Hamiltonian (4) in a harmonic basis set of frequency $`\omega _0`$:

$$|n=\underset{i=1}{\overset{i_{max}}{}}c_{ni}|i_{\omega _0},$$
(22)
These states are ordered according to their energy: $`E_nE_{n+1}`$ as exemplified in Table I. Within this scheme one obtains a very simple expression for $`T(n,m)`$

$$T(n,m)=E_n\delta _{nm}\underset{i,j=1}{\overset{i_{max}}{}}c_{ni}c_{nj}\frac{1}{2M}i|\widehat{P}^2|j.$$
(23)
It is recognized that any discrete basis set approximation to a system with a band spectrum will lack features of the original problem which have to be carefully analyzed. Anyway, alternatives to a discrete basis approach may be cumbersome and not necessarily yield a better approach to understand general properties of experimental data. While the discrete basis approach is exact for transitions between the lowest levels, which have a negligible width, eigenstates belonging to a band of measurably width should be treated with special care. Thus, we shall assume that matrix elements $`T(\nu ,\mu )`$ involving states with energies $`_\nu `$ and $`_\mu `$, so that $`E_n(E_nE_{n1})/2_\nu E_n+(E_{n+1}E_n)/2`$ with an analogous expressions for $`_\mu `$, are well approximated by $`T(n,m)`$.

Within this scheme the equations which describe the probability $`P(n)`$ of finding an atom in level $`n`$, given the transition rates $`R_{mn}`$ are

$$\dot{P}(n)=\underset{m}{}R_{mn}^{(1)}(P(m)P(n))$$
(24)
in the first order perturbation theory scheme, and the finite difference equations

$`P_n(t)`$ $`=`$ $`P_n(t_0)+`$ (27)
$`{\displaystyle \underset{m}{}}R_{mn}^{(1)}(P_m(t_0)P_n(t_0))(tt_0)+`$
$`{\displaystyle \underset{m}{}}R_{mn}^{(2)}(P_m(t_0)P_n(t_0))(tt_0)^2,`$
valid up to second order time-dependent perturbation theory whenever $`tt_0`$. Both sets of equations are subjected to the condition

$$\underset{n}{}P(n)=1.$$
(28)
Now, according to Eqs. (13) and Eq (21), the evaluation of $`R_{nm}^{(r)}`$ also requires the specification of the spectral density $`S(\omega )`$. In the problem under consideration, the discrete labels $`m,n`$ are used to calculate interband transitions which are actually spectrally broad. This broadening might arise not only from the band structure of the energy spectra associated with the Hamiltonian Eq. (4), but also from other sources, which we will discuss below. Broad spectral lines can be introduced in our formalism by defining an effective spectral density $`S_{\mathrm{eff}}(\omega )`$ which should incorporate essential features of this broadening without simulating specific features. Having this in mind, an effective Gaussian density of states $`𝒮_n(\omega )`$ is associated to each level $`|n`$ of energy $`E_n`$

$$𝒮_n(\omega )=\frac{1}{\sqrt{2\pi }\sigma _n}e^{{\displaystyle \frac{(\mathrm{}\omega E_n)^2}{2(\mathrm{}\sigma _n)^2}}}.$$
(29)
The spectral effective density $`S_{\mathrm{eff}}(\omega _{nm})`$ associated to the transition $`mn`$ is obtained by the convolution of $`𝒮_n(\omega )`$ with $`𝒮_m(\omega )`$ and with the excitation source spectral density $`S(\omega )`$. For a monochromatic source the latter is also taken as a Gaussian centered at the modulation frequency that once integrated over all frequencies yields the square of the intensity of the modulation source. The net result is that $`S_{\mathrm{eff}}(\omega _{nm})`$ has the form

$$S_{\mathrm{eff}}(\omega )=S_0e^{{\displaystyle \frac{(\omega \omega _{\mathrm{eff}})^2}{2\sigma _{\mathrm{eff}}^2}}}$$
(30)
with $`\omega _{\mathrm{eff}}`$ determined by the modulation frequency $`\mathrm{\Omega }`$ and the energies $`E_n`$ and $`E_m`$. The effective width $`\sigma _{\mathrm{eff}}`$ contains information about the frequency widths of the excitation source and those of each level.

## IV Comparison with experimental results.

We have tested the procedure described in last section to model parametric excitation in a specific experiment conducted at LENS. In this experiment 40K fermionic atoms are trapped in a 1D lattice, realized retroreflecting linearly polarized light obtained from a single–mode Ti:Sa laser at $`\lambda `$=787 nm, detuned on the red of both the D1 and D2 transition of potassium, respectively at 769.9 nm, and 766.7 nm. The laser radiation propagates along the vertical direction, to provide a strong confinement against gravity. The laser beam is weakly focused within a two-lens telescope to a waist size $`w_090\mu `$m, with a Rayleigh length $`z_R`$=3 cm; the effective running power at the waist position is $`P`$=350 mW.

The trap is loaded from a magneto-optical trap (MOT), thanks to a compression procedure already described in , with about $`5\times 10^5`$ atoms at a density around $`10^{11}`$cm-3. The typical vertical extension of the trapped atomic cloud, as detected with a CCD camera (see Fig. 1), is $`500\mu `$m, corresponding to about 1200 occupied lattice sites with an average of 400 atoms in each site. Since the axial extension of the atomic cloud is much smaller than $`z_R`$, we can approximate the trap potential to

$$V(r,z)=V_0e^{{\displaystyle \frac{2r^2}{w_0^2}}}\mathrm{cos}^2(kz);k=2\pi /\lambda ,$$
(31)
thus neglecting a 5% variation of $`V_0`$ along the lattice. The atomic temperature in the lattice direction is measured with a time–of–flight technique and it is about 50 $`\mu `$K.

In order to parametrically excite the atoms we modulate the intensity of the confining laser with a fast AOM for a time interval T$``$100 ms, with a sine of amplitude $`ϵ`$=3% and frequency $`\mathrm{\Omega }`$. The variation of the number of trapped atoms is measured by illuminating the atoms with the MOT beams and collecting the resulting fluorescence on a photomultiplier. In Fig.2 the fraction of atoms left in the trap after the parametric excitation is reported vs the modulation frequency $`\mathrm{\Omega }/2\pi `$. Three resonances in the trap losses are clearly seen at modulation frequencies 340 kHz, 670 kHz and 1280 kHz. By identifying the first two resonances with the lattice vibrational frequency and its first harmonic, respectively, we get as first estimate $`\omega _02\pi \times 340`$kHz. As we will show in the following, these resonance are actually on the red of $`\omega _0`$ and 2$`\omega _0`$, respectively, and therefore a better estimate is $`\omega _02\pi \times 360`$kHz. Therefore the effective trap depth is, from Eq. (9), $`V_0185\mu `$K $`10.5\mathrm{}\omega _0`$. Since the atomic temperature is about $`V_0`$/3.5, we expect that most of the energy levels of the lattice have a nonnegligible population and therefore the anharmonicity of the potential could play an important role in the dynamics of parametric excitation. Note that the third resonance at high frequency, close to 4$`\omega _0`$, is not predicted from the harmonic theory. It is possible to observe also a much weaker resonance in the trap losses around 1.5 kHz, which we interpret to be twice the oscillation frequency in the loosely confined radial direction. Anyway, in the following we will focus our attention just on the axial resonances.

As discussed in the previous Section, the overall width of the excitation assumed for our model system could play an important role in reproducing essential features of experimental data. Since the source used in the experiment has a negligible line width, it is necessary to model just the broadening of the atomic resonances. The spread of the transition energies due to the axial anharmonicity is reported in Table I, while the broadening of each energy level, due to the periodic character of the sine potential, is estimated using Eq. (8). We now note that the 1D motion assumed in Section II is not completely valid in our case, since the atoms move radially along a Gaussian potential. Since the period of the radial motion is about 500 times longer than the axial period, the atoms see an effective axial frequency which varies with their radial position, resulting in a broadening of the transition frequency. Other sources of broadening are fluctuations of the laser intensity and pointing, and inhomogeneities along the lattice. We note that also elastic collisions within the trapped sample, which tend to keep a thermal distribution of the trap levels population, can contribute to an overall broadening of the loss resonances. Since it is not easy to build a model which involve all these sources, we introduce semiempirically an effective broadening for the $`r`$-th level (see Eq. (29)). Recognizing that the width could be energy dependent we considered the simple expression

$$\sigma _r^2=\lambda _1\left(\frac{E_r}{V_0}\right)^p+\lambda _0$$
(32)
for several values of the constants $`\lambda _1,\lambda _0`$ and $`p`$. When $`p=0`$, i.e. for a constant value of the band width we were not able to reproduce the general experimental behavior reported in Fig.2. The best agreement between the simulation and the experimental observations is obtained for $`\lambda _0=0.0002`$, $`\lambda _1=0.0135`$, in units of $`\omega _0^2`$, and $`p=3`$. Similar results are obtained also for slightly higher (lower) values of $`\lambda _{0,1}`$ together with slightly higher (lower) values of the power $`p`$. In Table I the resulting widths are shown for the lower twelve levels. Note that we have intentionally excluded levels 11, 12 and 13 from the calculation, since their intrinsic width is so large that the atoms can tunnel out of the trap along the lattice in much less than 100 ms . Anyway, the inclusion of these levels proved not to change substantially the result of the simulation.

The comparison of experimental and theoretical results is made in Fig. 3; the abscissa for the experimental data has been normalized by identifying $`\omega _0`$ with $`2\pi \times 360`$kHz. As already anticipated, the principal resonance in trap losses appears at $`\mathrm{\Omega }1.85\omega _0`$. This result follows from the fact that the excitation of the lowest trap levels is not resulting in a loss of atoms, as it would happen for a harmonic potential. On the contrary, the most energetic atoms, which have a vibrational frequency smaller than the harmonic one, are easily excited out of the trap. The asymmetry of the resonances, which has been observed also in , is well reproduced in the calculations and it is a further evidence of the spread of the vibrational frequencies. The first interesting result obtained by our study of parametric excitation is therefore the correction necessary to extract the actual harmonic frequency from the loss spectrum. For the specific conditions of the present experiment, we find indeed that the principal resonance in the trap losses appears at $`\mathrm{\Omega }1.85\omega _0`$. Anyway, the calculation shows that the resonance is nearby this position for all the explored values of $`\lambda _{0,1}`$ and $`p`$ also for deeper traps, up to $`V_0`$=25 $`\mathrm{}\omega _0`$, and therefore it appears to be an invariant characteristic of the sinusoidal potential.

The result of the numerical integration of Eqs. (27) reported in Fig.3 reproduces relatively well the subharmonic resonance, which in the harmonic case would be expected at $`\omega _0`$. On the contrary, both experiment and calculation show that the actual position of the resonance is $`\mathrm{\Omega }0.9\omega _0`$. It must be mentioned that the accuracy of these results is restricted by the finite difference character of Eqs. (27) and by the fact that some noise sources which have not been included could be resonant at a nearby frequency. In particular, a possible modulation of the laser pointing associated to the intensity modulation is expected to be resonant at $`\mathrm{\Omega }=\omega _0`$ in the harmonic problem , and it could play an analogous role in our sinusoidal lattice.

The higher order resonance around 3.5 $`\omega _0`$ observed in the experiment is also well reproduced by the calculations based on first order perturbation theory. Note that a simpler approximation to the confining potential by a quartic potential $`V_Q(z)=k_2z^2+k_4z^4`$ would yield a resonance around $`4\omega _0`$ and not $`3.5\omega _0`$. Anyway, it is possible to understand qualitatively one of the features of this resonance considering a quartic perturbation of the form $`ϵ(t)V_Q`$ to a harmonic potential. In this case the ratio of the transition rates at the 2$`\omega _0`$ and 4$`\omega _0`$ resonances is set by Eqs.13 to

$$|T(n\pm 2,n)|^2/|T(n\pm 4,n)|^2\frac{V_0^2}{\omega _0^2}.$$
(33)
This result can qualitatively explain the absence of the corresponding high order resonance in the radial excitation spectrum (see Fig.2): since the radial trap frequency is a factor 500 smaller than the axial one, the relative strength of such radial anharmonic resonance is expected to be suppressed by a factor (500)2. In conclusion, high harmonics resonances, which certainly depend on the actual shape of the anharmonic potential, are expected to appear only if the spring constant of the trap is large.

In Fig.4 theoretical and experimental results for the evolution of the total population of trapped atoms at the resonant exciting frequency $`\mathrm{\Omega }=2\omega _0`$ are shown. Although there is a satisfactory agreement between the model and the experiment, we notice that experimental data exhibit a different rate for the loss of atoms before and after 100 ms. This change is probably due to a variation of the collision rate as the number of trapped atoms is modified, which cannot be easily included in the model. The comparison of the experimental evolution of the trap population with and without modulation shows the effectiveness of the excitation process in emptying the trap on a short time-scale.

We have also simulated the energy growth of the trapped atoms due to the parametric excitation, which is reported in Fig. 5. Our calculations show nonexponential energy increase in contrast with what expected in the harmonic approximation, Eq. (19). The fast energy growth at short times is related to the depopulation of the lowest levels, which are resonant with the 2$`\omega _0`$ parametric source. The saturation effect observed for longer times is due to the fact that the resonance condition is not satisfied for the upper levels so that they do not depopulate easily.

## V Conclusions.

We have studied both theoretically and experimentally the time evolution of the population of atoms trapped in a 1D sinusoidal optical lattice, following a parametric excitation of the lattice vibrational mode. In detail, we have presented a theoretical model for the excitation in an anharmonic potential, which represents an extension of the previous harmonic models, and we have applied it to the actual sinusoidal potential used to trap cold potassium atoms. The simulation seems to reproduce relatively well the main features of both the spectrum of trap losses, including the appearance of resonances beyond 2$`\omega _0`$, and the time evolution of the total number of trapped atoms.

By comparing the theoretical predictions and the experimental observations the usefulness of a parametric excitation procedure to characterize the spring constant of the trap has been verified. Although the loss resonances are red-shifted and wider than what expected in the harmonic case, the lattice harmonic frequency can be easily extracted from the experimental spectra, to estimate useful quantities such as the trap depth and spring constant.

We have also made emphasis on the need of modeling the broadening of bands with a negligible natural width in order to reproduce the observed loss spectrum. In a harmonic model this broadening is not necessary since the equidistant energy spectrum guarantees that a single transition energy characterizes the excitation process. We think that most of the broadening in our specific experiment is due to the fact that the actual trapping potential is not one-dimensional, and also to possible fluctuations and inhomogeneities of the lattice.

To conclude, we note that the dynamical analysis we have made can be easily extended to lattices with larger dimensionality, and also to other potentials, such as Gaussian potentials, which are also commonly used for optical trapping.

We acknowledge illuminating discussions with R. Brecha. This work was supported by the European Community Council (ECC) under the Contracts HPRI-CT-1999-00111 and HPRN-CT-2000-00125, and by MURST under the PRIN 1999 and PRIN 2000 Programs.