format: text

text:
# Predicting Neutron Production from Cosmic-ray Muons

## I introduction

Fast neutrons from cosmic-ray muons represent an important background for low-energy underground experiments such as searches for proton decay and dark matter, and low-energy neutrino oscillation. Unlike charged hadrons which can be tagged by a veto detector system, neutrons usually cannot be identified until they are finally captured, mimicking the signal. The occasional neutron scattering with the sensitive material of the detector and the long lifetime of neutrons in the detector and surrounding materials further complicate the situation. For example, the Palo Verde reactor neutrino oscillation experiment found such neutrons to be their dominant background , and a similar situation is expected at the ultra-long baseline detector KamLAND . Low-energy solar neutrino experiments such SNO and Borexino also have to estimate such backgrounds and the understanding of neutron backgrounds may be relevant in resolving the controversy between the CDMS  and DAMA  results on dark matter searches. Finally, low energy accelerator experiments at shallow depths, such as Karmen, LSND and OrLaND have also similar problems.

Although the total neutron yield from cosmic muon spallation has been measured by several experiments, contradictory results are given in the literature . Theoretical calculations  are also not consistent with each other and with data. The fact that primary neutrons, pions and protons can all produce secondary neutrons through hadronic interactions makes analytical calculations very difficult. A simple cascade model  suggests that the number of nuclear cascade products such as neutrons, pions and protons, increases with the average muon energy approximately like $`\mathrm{E}_\mu ^{0.7}`$. This formula agrees with most measurements  except the recent one by the LVD collaboration .

In addition to these problems with neutron yields, the few measurements of the neutron energy spectrum  are not well reproduced by theoretical calculations . The interpretation of experimental data is complicated by the fact that the neutron energy spectrum depends upon the muon spectrum that, in turn, is a function of the depth at which the measurement is carried-on and the geometrical configuration of the underground site. There is a broad range of results reported in the literature. Barton  suggests that the spectrum of neutrons from hadronic cascade follows $`E^{1/2}`$ between 10-50 MeV, while the spectrum of neutrons from $`\pi ^{}`$ capture follows a flat spectrum up to 100 MeV. Perkins  suggests that the neutron spectrum from muon spallation follows $`E^{1.6}`$. The combination of $`(9.7E^{1/2}+6.0e^{E/10})`$ has been used in a measurement  at a shallow site. It has also been suggested  to use proton and neutron spectra following $`E^{1.86}`$ as measured at accelerators for photo-nuclear interactions . Experimentally the Karmen experiment reported a visible energy spectrum following $`e^{E/39}`$ for spallation neutrons  and the LVD experiment reported a visible energy spectrum following $`E^1`$.

In order to put some order in this area, we have studied the neutron yield, the neutron energy spectrum, the multiplicity and the angular distribution using the 1999 version of the FLUKA  Monte Carlo program which is expected to give an accurate description of all the processes involved. It is our intention to provide an unified approach and obtain a reliable estimate of neutron background for experiments at all depths. While we perform the simulation in a particular detector configuration, results should be applicable to a range of substances, including rock and other shielding materials.

## II Muon spallation models

Fast neutrons from cosmic-ray muons are produced in the following processes:

* Muon interactions with nuclei via a virtual photon producing a nuclear disintegration. This process is usually referred to as “muon spallation” and is the main source of theoretical uncertainty.
* Muon elastic scattering with neutrons bound in nuclei.
* Photo-nuclear reactions associated with electromagnetic showers generated by muons.
* Secondary neutron production following any of the above processes.

Processes b) and c) are reasonably well understood while a) and d) are the root of the difficulties described in previous calculations. Neutrons can be also produced from muons which stop and are captured, resulting in highly excited isotopes emitting one or more neutrons. This process is reasonably well understood and its contribution to total neutron yield can be calculated. All the experimental results referred to in this paper do not include these neutrons since they can be easily identified and eliminated. Neutrons produced from neutrinos are negligibly small at the depths considered, and thus are not discussed in this paper.

The muon spallation process is schematically illustrated in Figure 1. The desired $`\mu N`$ cross section is then calculated as

$$\sigma _{\mu N}=\frac{N_\gamma (\nu )\sigma _{\gamma N}^{virt}(\nu )}{\nu }𝑑\nu $$
(1)
where $`\nu =EE^{}`$, $`\mathrm{E}`$ and $`\mathrm{E}^{}`$ are energies of initial and final muons, and $`\mathrm{N}_\gamma (\nu )`$ the virtual photon energy spectrum. Theoretical calculations often treat the virtual photons according to the Weizsäcker-Williams approximation , in which the passage of a charged particle in a slab of material produces the same effects as of a beam of quasi-real photons. A general expression of the Weizsäcker-Williams formula is given in ref. :

$`N_\gamma (\nu )`$ $`=`$ $`{\displaystyle \frac{\alpha }{\pi }}[{\displaystyle \frac{E^2+E^2}{p^2}}\mathrm{ln}{\displaystyle \frac{EE^{}+PP^{}m^2}{m\nu }}`$ (2)
$``$ $`{\displaystyle \frac{(E+E^{})^2}{2P^2}}\mathrm{ln}{\displaystyle \frac{(P+P^{})^2}{(E+E^{})\nu }}{\displaystyle \frac{P^{}}{P}}]`$ (3)
where $`\mathrm{m}`$ is the muon mass, and $`\mathrm{P}`$ and $`\mathrm{P}^{}`$ are momentum of initial and final muons.

Since in the above approximation, it is assumed that the $`\gamma N`$ cross section is the same for real and virtual photons, the measured $`\gamma N`$ cross section can be used to calculate the $`\mu N`$ cross section in Eq.(1). At low muon energy the situation is more complicated. Here, the virtuality of the photon becomes comparable to its energy and cannot be neglected. It follows that the Weizsäcker-Williams approximation can no longer be used. In addition, the interaction of the virtual photon with the nucleus is a collective excitation of the nucleus (Giant Dipole Resonance, GDR) rather than a single photon-nucleon interaction. This implies that the GDR model would have to be applied to virtual photons introducing further theoretical and technical complications. However, it might be reasonable to assume that neutron production by low-energy muon interactions is small as compared to neutron photoproduction by low-energy bremsstrahlung photons and adds therefore only a minor contribution to the total neutron yield.

In addtional to these assumptions, there are a number of problems associated with analytical calculations: first, they cannot reliably calculate all daughter products for every nucleus if the $`\gamma N`$ interaction is very violent so that the nucleus becomes highly excited; second, they cannot properly take into account secondary neutron production. Hence, while these calculations provide useful guidance and, at shallow depths, where hadronic shower effects are small, they may even give quantitatively sound predictions , in general they cannot be considered particularly reliable.

Monte Carlo approaches are commonly used to properly model hadronic cascades. Currently the most complete code to describe both hadronic and electromagnetic interactions up to 20 TeV is FLUKA . In this program, different physical models, or event generators, are responsible for the various aspects of particle production at different energies . High-energy hadronic interactions are described based on the Dual Parton Model followed by a pre-equilibrium-cascade model. In addition, models for nuclear evaporation, break-up of excited fragments and $`\gamma `$-deexcitation treat the disintegration of excited nuclei. Hadronic interactions of photons are simulated in detail from threshold (Giant Dipole Resonance interactions) up to TeV-energies (Vector Meson Dominance Model). For nuclei up to copper, measured photo-nuclear cross sections in the low energy region are used . Hadronic interactions of muons are based on the Weizsäcker-Williams approximation as formulated by Bezrukov and Bugaev . A spectrum of virtual photons is generated which interact with the nuclei similar to real photons. Due to the above theoretical and technical complications in the description of hadronic interactions of virtual photons at very low energies the simulation is restricted to photon energies above the delta resonance threshold. The implementation of hadronic interactions of muons has been shown to give reliable predictions for the MACRO experiment .

In the following FLUKA is used to obtain a consistent and complete estimate of neutron production from cosmic muons. We model a simple cubic detector filled with liquid scintillator $`C_nH_{2n+2}`$, where $`n`$ is taken to be 10. Muons with monochromatic energy are tracked in the detector and secondary neutrons and other hadrons are analyzed.

## III Neutron yield

The total neutron yield is probably the most measured quantity in our problem. There are many experimental results from different depths which can be compared with the model.

All neutrons, either primary or secondary, are included. Double counting due to neutron scattering or neutron spallation is carefully avoided. Due to the limited size of the detection volume in the simulation, some neutrons can escape, resulting in fewer secondary neutrons. The total number of neutrons thus depends on the size of the detector. The detector size has to be limited so that the muon energy-loss is small and the initial muon energy can be used as a constant. To correct this problem, we run the simulation with different detector sizes at each muon energy, and fit the neutron yield as a function of percentage of neutrons which escape the detector. An exponential behavior is found and the total number of neutrons can be extrapolated. The corrected neutron yield as a function of the muon energy is shown in Fig. 2. The neutron yield per muon can be fit as

$$\mathrm{N}_n=4.14\times E_\mu ^{0.74}\times 10^6neutron/(\mu gcm^2)$$
(4)
where E is in GeV. This relation is consistent with $`\mathrm{E}_\mu ^{0.7}`$ law suggested in ref. .

While many experiments report their results as a function of the detector’s depth underground, clearly the proper physical parameter is the mean muon energy at the detector. The conversion between the two quantities is not entirely trivial as the average energy depends upon the geometry of the overburden, particularly in the case of deep sites. Here for consistency we report all experimental results as a function of the average muon energy $`\overline{}E_\mu `$. For measurements performed at deep laboratories such as Gran Sasso or Mont Blanc the conversion is given in the original papers , while for shallow sites we estimate the average muon energy-loss using the simple relation  and assuming a flat geometry:

$$\frac{dE_\mu }{dX}=\alpha E_\mu /\xi $$
(5)
where $`\mathrm{E}_\mu `$ is in GeV, $`\alpha =(1.9+0.08log(E_\mu /m_\mu ))\times 10^3GeV/(gcm^2)`$ and $`\xi =2.5\times 10^5g/cm^2`$ in rock. All known experimental neutron-production rates are reported in Figure 2 along with the FLUKA predictions. The agreement is substantially better than obtained with previous calculations .

Given this agreement we can proceed to extract from the Monte Carlo simulation information on the origin of the neutrons. This is of great interest as it can provide hints for a better theoretical understanding of the processes involved. In Figure 3 we show the fractional composition of the neutrons by origin. As expected the fraction of neutrons from primary processes such as muon spallation and real photo-nuclear interactions decreases with energy, while secondary neutrons from neutron and proton spallation, from pion absorption and from other minor processes such as $`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }`$ decays, increase in relative importance with energy.

The $`\pi ^+`$ yield in muon spallation can be estimated from analytical calculations more reliably than that of neutrons. This is because 1) secondary processes for neutrons are more important; 2) there is an uncertainty in the direct neutron production from virtual or real pions. Table I shows the comparison of our simulation with an analytical calculation  which does not include secondary pions from hadronic showers. In the experimental data , $`\pi ^+`$’s are identified through their $`\pi ^+\mu ^+\mathrm{e}^+`$ decays, and hadronic shower effects are not corrected for. In the pion case, good agreement between data, FLUKA and the analytical calculation can be seen at shallow depth, where secondary hadrons from showers are not important. We obtained the pion yield per muon as:

$`\mathrm{N}_{\pi ^+}=4.45\times E_\mu ^{0.80}\times 10^7pion/(\mu gcm^2)`$ (6)
where $`\mathrm{E}_{mu}`$ is in GeV. The $`\mathrm{E}_\mu ^{0.7}`$ law appears to be universal as suggested by Ryazhskaya  and consistent with results from ref. .

## IV Neutron energy spectrum

The neutron energy spectrum is particularly controversial, with a wide range of results reported in theoretical calculations and in the few experimental measurements. In Figure 4 we show some of the energy spectra obtained with FLUKA. Each histogram is fitted to the universal empirical function

$$\frac{dN}{dE_n}=A\left(\frac{e^{7E_n}}{E_n}+B(E_\mu )e^{2E_n}\right)$$
(7)
where $`A`$ is a normalization factor, and

$$B(E_\mu )=0.520.58e^{0.0099E_\mu }.$$
(8)
This simple function reproduces fairly well the FLUKA distributions with $`\chi ^2`$ per degree of freedom of 3.9, 4.5, 9.6 and 3.7 for 11, 20, 90 and 270 GeV respectively.

We are aware of two direct measurements that can be compared to our calculations, as shown in Fig. 5. The Karmen collaboration measured neutron energy spectrum up to 50 MeV. Within this modest range, their parameterization $`\mathrm{N}_{soft}e^{E_n/2.1}+N_{hard}e^{E_n/39}`$, is in a reasonable agreement with our result. They attribute all the soft component to muon capture but it seems that muon spallation also produces some soft neutrons. The LVD experiment reported a $`E_n^1`$ spectrum up to 400 MeV, also in a reasonable agreement with our result.

Some of the theoretical estimates of neutron energy spectrum are comparable to our results at low energies while our calculation gives results over a much wider range of energies. For example, the power laws $`E_n^{1.6}`$  and $`E_n^{1.86}`$  agree with our results up to 400 MeV. Other functions, such as $`(9.7E_n^{1/2}+6.0e^{E_n/10})`$, or that suggested by Barton , are significantly different.

## V Neutron multiplicity and angular distribution

The neutron multiplicity is probably the least known quantity in the neutron production problem. In most cases muon spallation only happens once and produces only a few primary hadrons. But these hadrons can shower and generate secondary hadrons, including neutrons. Using our simulation we have found that in some cases the number of secondary neutrons exceeds 50. The average number of neutrons is about 3 for a 11 GeV muon, and it increases to about 7 at a muon energy of 385 GeV. Figure 6 shows the neutron multiplicity distributions at different muon energies from FLUKA, together with the universal empirical parameterization:

$$\frac{dN}{dM}=A(e^{A(E_\mu )M}+B(E_\mu )e^{C(E_\mu )M})$$
(9)
where $`M`$ is the multiplicity,

$`A(E_\mu )`$ $`=`$ $`0.085+0.54e^{0.075E_\mu }`$ (10)
$`B(E_\mu )`$ $`=`$ $`{\displaystyle \frac{27.2}{1+7.2e^{0.076E_\mu }}}`$ (11)
$`C(E_\mu )`$ $`=`$ $`0.67+1.4e^{0.12E_\mu }`$ (12)
The $`\chi ^2`$’s per degree of freedom are 0.6, 2.0,1.5 and 1.5 for 11, 20, 90 and 270 GeV respectively.

There are only few experimental results  on multiplicity. The Palo Verde experiment at the shallow depth corresponding to a mean muon energy of about 16.5 GeV, observed a two-neutrons to one-neutron ratio between 5 to 10, depending on the assumption on three neutron yield. Similar results are also reported by Bezrukov et al. . These numbers appear to be substantially larger than the FLUKA prediction of 2. While more data would be desirable to better understand this discrepancy it is possible that the data is affected by the incomplete efficiency matrix for four or more neutrons, whose contribution, previously assumed small, seems significant from our simulation.

There are no experimental results on the neutron angular distribution with respect to muons. Experimental results  on $`\gamma +^{12}Cp+X`$ have been used  as a first approximation of the neutron angular distribution. The distribution is expected to be forward-peaked, smoothed somewhat by the contribution of secondary neutrons. Our simulation is well parameterized by the angular distribution

$$\frac{\mathrm{dN}}{\mathrm{dcos}\theta }=\frac{A}{(1cos\theta )^{0.6}+B(E_\mu )}$$
(13)
where $`B(E_\mu )=0.699E_\mu ^{0.136}`$. Figure 7 shows the $`\mathrm{cos}\theta `$ distribution at the usual set of energies together with the above parameterization. The $`\chi ^2`$ per degree of freedom are 0.51, 0.51, 0.72 and 0.72 for 11, 20, 90 and 270 GeV respectively.

## VI Summary

We obtained for the first time a complete description of neutron production from cosmic-ray muons using the FLUKA Monte Carlo program. Results have been compared with existing data and some analytical calculations reported in literature. With a few exceptions, our results agree well with data, and our predictions cover an energy range much broader than what has been discussed before.

Analytical parameterizations have been obtained for neutron yield, energy spectrum, multiplicity and angular distribution. These formulae can be used as a starting point to estimate neutron backgrounds from cosmic-ray muons when the muon energy spectrum is known. For detailed calculation of muon-induced neutron backgrounds in low-energy underground experiments, a complete Monte Carlo simulation is needed to account for the correlation among different variables and can be setup using FLUKA. Intermediate processes such as neutron scattering before neutron spallation are also automatically taken into account in this method that has great promise to improve the quality of predictions.

## VII Acknowledgments

We would like to thank P. Vogel for many useful discussions. This work was supported in part by DoE grant DE-FG03-96ER40986 and DE-AC03-76SF00515.


text:
# Stability of exact solutions of the defocusing nonlinear Schrödinger equation with periodic potential in two dimensions

## 1 Introduction

The cubic nonlinear Schrödinger equation with attractive or repulsive nonlinearity and a potential is used as a mean-field model for the dynamics of a dilute-gas Bose Einstein condensate (BEC) . In this case, the equation is often referred to as the Gross-Pitaevskii equation. These BECs are of interest in both the theoretical and experimental physics community: they are examples of macroscopic quantum phenomena which display phase coherence  and which have potential applications in such diverse areas as quantum logic , matter-wave diffraction , and matter-wave transport .

After rescaling of the physical variables, the governing equation is

$$i\frac{}{t}\psi (\stackrel{}{x},t)=\frac{1}{2}\mathrm{\Delta }\psi (\stackrel{}{x},t)+\alpha \left|\psi (\stackrel{}{x},t)\right|^2\psi (\stackrel{}{x},t)+V(\stackrel{}{x})\psi (\stackrel{}{x},t),$$
(1)
where $`\mathrm{\Delta }`$ denotes the Laplacian operator, $`\psi (\stackrel{}{x},t)`$ is the macroscopic wave function of the condensate in one, two or three dimensions, with $`\stackrel{}{x}`$ defined as $`x`$, $`(x,y)`$ or $`(x,y,z)`$ respectively. The real-valued function $`V(\stackrel{}{x})`$ is an experimentally generated macroscopic potential. The parameter $`\alpha `$ determines whether (1) is repulsive ($`\alpha =1`$, defocusing nonlinearity), or attractive, ($`\alpha =1`$, focusing nonlinearity). Note that for BEC applications, both signs of $`\alpha `$ are relevant.

Many BEC experiments use harmonic confinement, but recently there has been interest in confinement of repulsive BECs using standing light waves, resulting in a sinusoidal potential in the Nonlinear Schrödinger equation. So far, most interest has focussed on the quasi-one-dimensional regime, in which the longitudinal dimension of the condensate is far greater than the transverse dimensions, which are themselves of the order of the healing length of the condensate. In this case, the governing mean-field equation is (1), with $`\alpha =1`$, in one dimension, and the experimentally generated potential $`V(x)=V_0\mathrm{sin}^2(x)`$. In , a number of exact solutions of this equation was constructed and their stability was investigated. In fact, a potential more general than a sinusoidal potential was considered:

$$V(x)=V_0\mathrm{sn}^2(x,k),$$
(2)
where $`\mathrm{sn}(x,k)`$ denotes the Jacobian elliptic sine function, with elliptic modulus $`k`$. In the limit $`k=0`$, this potential reduces to a sinusoidal potential. It is important to realize that for values of $`k`$ up to $`k=0.9`$, the shape of the potential is virtually indistinguishable from a sinusoidal one. On the other hand, considering values of $`k`$ close to $`1`$, $`i.e.,`$ $`k>0.999`$ gives periodic potentials with well-separated peaks, allowing the study of BEC dynamics in an entirely new regime.

Currently, no experiments are being performed where a BEC is trapped in a higher-than-one-dimensional periodic potential. However, the interest in the applications mentioned above strongly suggests that these experiments may eventually take place. Already, there are theoretical investigations of BECs in multidimensional lattice potentials , suggesting the realization of such experiments. Although motivated by the developments in BECs, in this paper we consider (1) with repulsive nonlinearity in two dimensions in its own right. Thus we consider

$$i\frac{\psi }{t}=\frac{1}{2}\left(\frac{^2\psi }{x^2}+\frac{^2\psi }{y^2}\right)+\left|\psi \right|^2\psi +V(x,y)\psi ,$$
(1)
with $`\psi =\psi (x,y,t)`$. As in one dimension , making the proper choices for the potential allows the construction of a large class of exact solutions. The family of potentials considered is

$`V(x,y)`$ $`=`$ $`\left(A_1\mathrm{sn}^2(m_1x,k_1)+B_1\right)(A_2\mathrm{sn}^2(m_2y,k_2)+B_2)+`$ (3)
$`+m_1^2k_1^2\mathrm{sn}^2(m_1x,k_1)+m_2^2k_2^2\mathrm{sn}^2(m_2y,k_2).`$
Here $`A_1,A_2,B_1,B_2,m_1`$ and $`m_2`$ are real constants. The two elliptic moduli $`k_1`$ and $`k_2`$ are in the interval $`[0,1]`$. The first term is a straightforward generalization of the one-dimensional potential (2), with an additive constant. The other terms are incorporated to facilitate the construction of exact solutions. Although this exact expression for the potential is necessary to allow the construction of exact solutions, it is the qualitative features, $`i.e.,`$ its periodicity and amplitude, that are most important. Numerical and analytical results throughout this paper demonstrate that the behavior of a solution in a lattice potential is largely independent of the quantitative features of the potential. Figure 1 displays the potential (3) for various values of its parameters.

In the next section, a family of exact solutions to (1) is given and their linear stability is discussed. Numerical results for various representative solutions are given in Section 3. The last section concludes the paper with a brief summary of the relevant mathematical and physical results.

## 2 Analytical Results

A judicious choice of the potential allows for the cancellation of the nonlinear term in (1) so that exact solutions can be constructed. Of course, one can always find a suitable potential by solving (1) for $`V(x,y)`$, given a certain $`\psi (x,y,t)`$. However, this generically results in time-dependent potentials and is hence not of interest. The kind of potentials we seek are two-dimensional generalizations of the one-dimensional potential (2). This one-dimensional potential is periodic. Its two-dimensional generalizations should be periodic in both its spatial variables. Furthermore, it is desirable that for a given potential more than one exact solution exists.

A derivation of the exact solutions is not presented in this paper. Rather we state the solution formulas and discuss them. For the potential (3), a class of exact solutions is given by

$$\psi (x,y,t)=r_1(x)r_2(y)e^{i\theta _1(x)+i\theta _2(y)i\omega t},$$
(4)
with

$`r_1^2(x)=A_1\mathrm{sn}^2(m_1x,k_1)+B_1,\theta _1(x)=c_1{\displaystyle _0^x}{\displaystyle \frac{dz}{A_1\mathrm{sn}^2(m_1z,k_1)+B_1}},`$ (5a)
$`r_2^2(y)=A_2\mathrm{sn}^2(m_2y,k_2)+B_2,\theta _2(y)=c_2{\displaystyle _0^y}{\displaystyle \frac{dz}{A_2\mathrm{sn}^2(m_2z,k_2)+B_2}},`$ (5b)
and

$`\omega `$ $`=`$ $`{\displaystyle \frac{1}{2}}m_1^2\left(1+k_1^2+k_1^2{\displaystyle \frac{B_1}{A_1}}\right)+{\displaystyle \frac{1}{2}}m_2^2\left(1+k_2^2+k_2^2{\displaystyle \frac{B_2}{A_2}}\right),`$ (6a)
$`c_1^2`$ $`=`$ $`m_1^2{\displaystyle \frac{B_1}{A_1}}\left(A_1+B_1\right)\left(A_1+k_1^2B_1\right),`$ (6b)
$`c_2^2`$ $`=`$ $`m_2^2{\displaystyle \frac{B_2}{A_2}}\left(A_2+B_2\right)\left(A_2+k_2^2B_2\right).`$ (6c)
Choosing the parameters $`A_1A_2`$, $`B_1/A_1`$, $`B_2/A_2`$, $`m_1`$, $`m_2`$, $`k_1`$ and $`k_2`$ determines the potential. Thus the solution class given by (5a-b) is a one-parameter family of solutions, with free parameter $`A_1/A_2`$. Existence of these solutions requires $`r_1^20,c_1^20`$ and $`r_2^20,c_2^20`$. This imposes conditions on the parameters: $`A_i0,B_i0`$ or $`B_i0,A_iB_iA_i/k_i^2`$, where $`i=1,2`$.

From (5a-b), it follows that $`r_1(x)`$ is periodic in $`x`$ with period $`2K(k_1)/m_1`$, where $`K(k)=_0^{\pi /2}𝑑z/\sqrt{1k_1^2\mathrm{sin}^2(z)}`$, and $`r_2(y)`$ is periodic in $`y`$ with period $`2K(k_2)/m_2`$. In general, $`\theta _1(x+2K(k_1)/m_1)\theta _1(x)+2n\pi `$, for some integer $`n`$. Thus, the solution (4) is usually not periodic in $`x`$ or in $`y`$. Imposing this periodicity requires a quantization of the phase $`\theta _1(x)`$ and $`\theta _2(y)`$ . It is unclear if such a quantization is possible, since only one free parameter is available to satisfy two conditions. There are two special cases in which phase quantization is not a concern. The first case results in trivial-phase solutions, for which $`c_1=0`$, $`c_2=0`$. The second case is the trigonometric limit, in which $`k_1=0`$, $`k_2=0`$.

The solution has trivial phase in the $`i`$-direction if $`c_i`$ is zero, where $`i=1(2)`$ corresponds to the $`x(y)`$-direction. There are three possibilities:

$`B_i`$ $`=`$ $`0:r_i=\sqrt{A_i}\mathrm{sn}(m_ix,k_i),`$ (7a)
$`B_i`$ $`=`$ $`A_i:r_i=\sqrt{A_i}\mathrm{cn}(m_ix,k_i),`$ (7b)
$`B_i`$ $`=`$ $`A_i/k_i^2:r_i=\sqrt{A_i/k_i^2}\mathrm{dn}(m_ix,k_i),`$ (7c)
where $`\mathrm{cn}(m_ix_i,k_i)`$ is the Jacobian elliptic cosine function, and $`\mathrm{dn}(m_ix_i,k_i)`$ denotes the third Jacobian elliptic function.

The solution is trigonometric in the $`i`$-direction if $`k_i`$ is zero. Then

$$r_i^2(x)=A_i\mathrm{sin}^2(x_i)+B_i,\mathrm{tan}(\theta _i)=\sqrt{1+\frac{A_i}{B_i}}\mathrm{tan}(m_ix_i),$$
(8)
and phase quantization is satisfied. Notice that it is possible for the solution to have trivial phase in one direction and be trigonometric in the other.

The main stability results from  can be generalized to two dimensions in a straightforward manner. The details of this will be presented elsewhere. Here a brief summary is given. Examining the linear stability of exact trivial-phase stationary solutions $`\psi (x,y,t)=r(x,y)\mathrm{exp}(i\omega t)`$ of equation (1) gives rise to the eigenvalue problem

$$\left(\begin{array}{cc}0& L_{}\\ L_+& 0\end{array}\right)\left(\begin{array}{c}u_1\\ u_2\end{array}\right)=\lambda \left(\begin{array}{c}u_1\\ u_2,\end{array}\right)$$
(9)
where the linear Schrödinger operators $`L_+`$ and $`L_{}`$ are

$`L_+`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(_x^2+_y^2\right)+3r^2(x,y)+V(x,y)\omega ,`$ (10a)
$`L_{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(_x^2+_y^2\right)+r^2(x,y)+V(x,y)\omega .`$ (10b)
Reasoning similar to the arguments in  leads to the following conclusions:

* If $`r(x,y)>0`$, then $`r(x,y)`$ is the ground state of $`L_{}`$ and $`\psi (x,y,t)`$ is linearly stable.
* If $`r(x,y)`$ has any zeros, it is not the ground state of $`L_{}`$. If in addition $`L_+`$ is a positive operator, then $`\psi (x,y,t)`$ is unstable.
* In all other cases, no conclusion is obtained from this stability analysis.

Primarily, this analysis immediately establishes the linear stability of a solution $`\psi (x,y,t)=r_1(x)r_2(y)\mathrm{exp}(i\omega t)`$, where both $`r_1(x)`$ and $`r_2(y)`$ are strictly positive. Thus, the solution with both the $`x`$\- and $`y`$-dependence specified by (7c) is linearly stable.

## 3 Numerical Results

In this section, the result of numerically solving (1) with initial conditions chosen from the exact solutions given in the previous section are discussed. The numerical procedure uses a fourth-order Runge-Kutta method in time and a filtered pseudospectral method in space. For each experiment, a small amount of white noise was added to the initial conditions as a perturbation.

Consider the solution

$$\psi (x,y,t)=\frac{\sqrt{A_1}\sqrt{A_2}}{k_1k_2}\mathrm{dn}(m_1x,k_1)\mathrm{dn}(m_2y,k_2)e^{i\omega t}.$$
(11)
where $`A_1<0`$ and $`A_2<0`$.

As stated in the previous section, this solution is linearly stable. This is confirmed by the numerics presented in Fig. 2. The three different columns give from left to right the dynamics of the solution, the same using contour plots, and the arctan of the Fourier power spectrum of the solution. In the first two columns, the bottom figure shows the potential $`V(x,y)`$. The stability of this solution is reminiscent of the stability of a plane wave solution of the two-dimensional defocusing nonlinear Schrödinger equation .

Next, consider the trivial-phase solution

$$\psi (x,y,t)=\sqrt{A_1}\sqrt{A_2}\mathrm{cn}(m_1x,k_1)\mathrm{cn}(m_2y,k_2)e^{i\omega t}.$$
(12)
The numerics of Fig. 3 shows that this solution is unstable. For the equation (1), the $`L^2`$-norm of the solution is conserved, hence so is the $`L^2`$-norm of its Fourier spectrum. This is not reflected in Fig. 3, due to the arctan transformation, which is used to diminish the range of the power spectrum. The onset of instability occurs between $`t=15`$ and $`t=30`$. A detailed look at this onset is given in Fig. 4.

The middle column of this figure shows $`||\psi (x,y,t)|^2|\psi (x,y,t=0)|^2|`$. The instability begins to develop around $`t=18`$. Before the instability occurs, $`|\psi (x,y,t)|^2`$ is equal to its initial condition, up to effects due to the added noise and numerical round-off, as is expected, since $`\psi (x,y,t)`$ is an exact stationary solution of (1). The solution

$$\psi (x,y,t)=\sqrt{A_1A_2}\mathrm{sn}(m_1x,k_1)\mathrm{sn}(m_2y,k_2)e^{i\omega t}.$$
(13)
is observed numerically to be unstable as well. Its dynamics is very similar to that of the solution (12). These results are consistent with the one-dimensional results obtained previously .

Analytical results for the stability of the nontrivial phase solutions are difficult to obtain, and so far numerical investigations are the only means of examining these solutions. Since phase quantization is a serious complication for the numerics, it is convenient to work with solutions (8) that are trigonometric in both directions since phase quantization is automatically satisfied. Thus, we consider solutions of the form

$$\psi (x,y,t)=\sqrt{\left(A_1\mathrm{sin}^2(m_1x)+B_1\right)\left(A_2\mathrm{sin}^2(m_2y)+B_2\right)}e^{i\theta _1(x)+i\theta _2(y)i\omega t},$$
(14)
with

$$\mathrm{tan}(\theta _1(x))=\sqrt{1+\frac{A_1}{B_1}}\mathrm{tan}(m_1x),\mathrm{tan}(\theta _2(y))=\sqrt{1+\frac{A_2}{B_2}}\mathrm{tan}(m_2y).$$
(15)
These solutions are stable or unstable depending on the offset parameters $`B_1`$ and $`B_2`$. As shown in Fig. 2 the off-set $`\mathrm{dn}(m_ix_i,k_i)`$-type solution is stable whereas in Fig. 3 the unstable $`\mathrm{sn}(m_ix_i,k_i)`$ and $`\mathrm{cn}(m_ix_i,k_i)`$ solutions without offset are unstable. These trivial phase solutions suggest that offset is essential for stability. In Fig. 5, $`A_1=A_2=1`$ and $`B_1=B_2=0.5`$ and the onset of instability of the nontrivial phase solution is observed. Here the chosen offset ($`B_1=B_2=0.5`$) is insufficient to stabilize the dynamics. In contrast, Fig. 6 has parameter values $`A_1=A_2=1`$ and $`B_1=B_2=1`$ and is stable. Thus with a higher offset, the nontrivial phase solution is stabilized. The boundary between the stable and unstable regions is difficult to calculate analytically and we rely on numerical simulations to determine the amount of offset required for stability. This behavior is again consistent with the one-dimensional case .

## 4 Summary and Conclusions

We considered the repulsive nonlinear Schrödinger equation with an elliptic function potential in two dimensions. Periodic solutions of this equation were found and their stability was investigated both analytically and numerically. Using analytical results for trivial-phase solutions, we showed that solutions with sufficient offset are linearly stable. This is confirmed with extensive numerical simulations on the governing nonlinear equation. Likewise, nontrivial-phase solutions are stable if they are sufficiently off-set.

Since our equations is a model for a Bose-Einstein condensate trapped in a lattice potential, our results imply that a large number of condensed atoms is sufficient to form a stable, periodic condensate. Physically, this implies stability of states near the Thomas–Fermi limit. The results are consistent with the one-dimensional trapping of a BEC condensate in a standing light wave. To quantify this phenomena, we consider the $`k=0`$ limit and note that in the trigonometric limit $`k0`$ the number of particles $`n`$ per potential well is given by $`n=|\psi (x,y,t)|^2=r_1^2(x)r_2^2(y)=(A_1/2+B_1)(A_2/2+B_2)`$, where $``$ denotes the spatial average. In the context of the BEC, and for a fixed atomic coupling strength, this means a large number of condensed atoms $`n`$ per potential well is sufficient to provide an offset on the order of the potential strength. This ensures stabilization of the condensate. Alternatively, a condensate with a large enough number of atoms can be interpreted as a developed condensate for which the nonlinearity acts as a stabilizing mechanism.

Acknowledgments: The authors acknowledge useful conversations with J. C. Bronski, L. D. Carr, R. Carretero-González, K. Promislow, and W. Reinhardt.