text:
# Stability of strange stars (SS) derived from a realistic equation of state.

## 1 Introduction

Recently there has been some excitement about the possibility that some compact stars are made from unusual forms of matter . Further, from an analysis of over 1 million seismic data reports sent to the U.S. Geological Survey in the years 1990-93, which were not associated with traditional epicentral sources, Anderson et al. (2002) claim to have found two events that ‘have the properties predicted for the passage of a strange quark nugget’ through the earth . Approximately 8000 separate seismic stations around the world are included in the database. If confirmed this would be a discovery similar to the detection of gamma ray bursts .

The best observational evidence for the existence of quark stars seems to come from some compact objects like the X-Ray burst sources SAX J1808.4$``$3658 (the SAX in short) and 4U 1728$``$34, the X-ray pulsar Her X-1 and the super burster 4U 1820$``$30. Among these the first is the most stable pulsating X-ray source known to man as of now. This star is claimed to be an ReSS with mass  $`1.3M_{}`$ and a radius of about 7 km. The mass of 4U 1728$``$34 is claimed to be less than 1.1 $`M_{}`$ in Li et al. , which places it much lower in the M-R plot (Fig.1). So it could be still gaining mass and shift to another stable point on the M-R graph. Thus in the model proposed in  there is a possible answer to the question posed by Franco : why are the pulsations of SAX not attenuated, as they are in the 4U 1728$``$34 ?

Fig.1 presents M-R relations for neutron as well as strange stars. The current phenomenology of compact objects could be interpreted to indicate that the mass of a star increases due to accretion so that the radius of a star changes from one point of the stable M-R curve to another. For neutron stars, exemplified by the EOS BBB 111this is one of the set calculated by Baldo, Bombaci and Burgeo  using a realistic nuclear equation of state. \- the curve on the right in Fig.1 \- a smaller mass would imply a larger radius for the star. If the mass of the star increases it should contract. Therefore, expansion due to an increase of mass, subsequent to accretion, is an unstable process for a neutron star. The expected behaviour for SS is directly opposite to that of neutron stars as Fig.1 shows and therefore may support the claim that some compact stars are ReSS.

Coupled to our claim are various other evidences for the existence of ReSS, such as explaining the compact M-R relations of the two candidates given in  (namely, the Her X-1 and the 4U 1820$``$30), as also the possible explanation of two kHz quasiperiodic oscillations in the 4U 1728 - 34 .

Recently we found that the matter of an ReSS may have diquarks on the surface which could account for the delayed emission of huge amounts of energy, after a thermonuclear catastrophe. This could be a possible scenario for superbursts .

Some of the stars like the SAX J1808.8$``$3658 or the PSR 1937$`+`$21 are fast rotors. ReSS have the possibility of withstanding high rotations which neutron stars or even bag SS cannot sustain. The maximum frequencies for the two EOS of D98 are 2.6 and 2.8 kHz respectively when they are on the mass shed limit (supramassive model) and 1.8 kHz and 2 kHz when they are in the normal evolutionary sequence as shown in Gondek-Rosińska et al. . In the present paper we further show that the ReSS are not only stable under fast rotation but also against radial oscillations.

The strange matter hypothesis has been used to postulate a scenario whereby a neutron star collapses to a quark star, simultaneously with a gamma ray burst - the so called phenomenon of a quark nova .

In the cosmic separation of phase scenario of Witten  SS are created along with baryons in a hot environment during the expansion of the early universe. The problem was investigated for ReSS  and it was found that they are formed at a temperature $`T70`$ $`MeV`$. Since it is self-sustaining, the system expands as the Universe cools, the ReSS expands like the Universe itself with cooling . The calculation also suggests that the shift in entropy due the change from normal to strange matter is not very large, indicating that the phase transition is relatively smooth.

We note that in a series of early papers van Paradijs () had noted that (i) the X-rays for bursters originate from stars with radii around 7 kms, assuming a canonical mass of 1.4 $`M_{}`$ for them and (ii) if one assumes a lower mass the estimated radii also becomes lower, which fits the M-R relation for the ReSS (Fig. 1)

Because of the extremely strong electric field stretching outside the stellar surface within $`10^{10}`$ cm  Usov  suggested that at a finite surface temperature ($`5\times 10^8K`$) of an SS one expects to have the creation of $`e^+e^{}`$ pairs on its surface. Such an effect could be an additional observational signature of SS with nearly bare quark surfaces. The strange matter equilibrium will not allow the pair production to quench, contrary to the comment in  as pointed out in . Usov claimed that the chemical potential of the electrons at the surface is large compared to their mass, being around 20 MeV and the mean velocity of the electrons is very high, so that the electric field will be restored very quickly. In our calculations we specifically find this to be $``$29 MeV, thus strongly supporting the conclusions of Usov.

Usov has further claimed that from the nature of the two step process, viz. $`e^+e^{}`$ pair production and subsequent $`\gamma `$ emission, the soft $`\gamma `$ repeaters are indeed very young SS  and their genesis may be due to the impact of comet like objects with these stars .

We must also add that radio pulsars may be SS as was suggested recently by Xu et al  for the PSR 0943+10 and all other drifting pulsars. Further, Kapoor and Shukre  used a remarkably precise observational relation for pulsar core component widths of radio pulsars to get stringent limits on pulsar radii, strongly indicating that some pulsars are strange stars. This is achieved by including general relativistic effects due to the pulsar mass on the size of the emission region needed to explain the observed pulse widths. A recent paper supporting their ideas is Xu, Xu and Wu  for PSR 1937$`+`$21,the fastest known radio pulsar.

The calculations for cold ReSS by Dey et al.  enables us to draw conclusions about chiral symmetry restoration (CSR) in QCD when the EOS is used to get ReSS fitting definite M-R relations . The empirical M-R relations were derived from astrophysical observations like luminosity variation and some properties of quasiperiodic oscillations from compact stars. The density dependence of the strong coupling constant can be deduced from the CSR described above .

Recently Glendenning  has argued that the SAX could be explained as a neutron star rather than a bare SS, not with any of the existing known EOS, but with a hypothetical one, satisfying however, the well-accepted restrictions based on general physical principles and having a core density about 26 $`\rho _0`$. Of course, such high density cores imply hybrid strange stars, subject to Glendenning’s assumption that such stars can exist with matching EOS for two phases. There is the further constraint that if the most compact hybrid star has a given mass, all lighter stars must be larger. It was found in Li et al.  that the star 4U 1728$``$34 may have a mass less than that of the SAX and yet have a smaller radius. Another serious difference is that the EOS of D98, using the formalism of large $`N_\mathrm{c}`$ approximation, indeed shows a bound state in the sense of having a minimum at about $`4.8n_0`$, whereas in Glendenning  one of the assumptions is that strange matter has no bound state.

## 2  ReSS model - the EOS

The model solves the relativistic mean field equation with a realistic qq interaction, the quark masses decreasing with density and restoring to their current mass values at high enough density. Moreover, the bare, confining qq interaction is expected to get screened in the medium. The inverse Debye screening length ($`D^1`$) is an increasing function of density. With these two conditions, namely the tendencies towards deconfinement and chiral symmetry restoration with density, we set out to obtain the required EOS for beta-equilibrated, charge neutral strange quark matter (SQM) containing uds quarks.

The plausibility condition for SQM to be preferred to ordinary matter is,

$$E_{ud}/A>E_{Fe^{56}}/A>E_{uds}/A$$
(1)
In Fig.2 we have shown that the energy per baryon for our equation of state (EOS1 of ) has a minimum at E/A = 888.8 $`MeV`$ compared to 930.4 of $`Fe^{56}`$. The pressure at this point is zero and this marks the surface of the star as can be seen when the well known TOV equation is solved. The curve clearly shows that the system can fluctuate about this minimum - so that the zero pressure point can vary. The EOS is parameterized to a linear form in  as

$$p=a(ϵϵ_0),$$
(2)
where $`a`$ and $`ϵ_0`$ are two parameters, $`a=0.463`$ and $`ϵ_0=1.15\times 10^{15}gmcm^3`$.
For EOS3 they are $`a=0.455`$, $`ϵ_0=1.33\times 10^{15}gmcm^3`$. Sharma et al  showed that the model possesses scaling properties similar to the bag model.

The bag model EOS, which satisfies the plausibility condition Eq.(1), must have a bag constant which cannot be higher than $`75MeV/fm^3`$ when the strange quark mass is taken to be $`150MeV`$. We have plotted this EOS as a dashed curve in Fig.2. The minimum in E/A is seen to occur at a low density, about $`2n_0`$ where $`n_0`$ is the normal nuclear matter density of $`0.17/fm^3`$. This kind of density may have been already reached in present day heavy ion collision experiments and yet no clear signature of stable strange matter has been observed. This problem does not occur with ReSS where the surface density is about $`5n_0`$. In addition to this, the bag does not fit the M-R data for the strange star candidates as has been repeatedly stated before . The minimum of energy of the bag model (dashed curve) is 921 $`MeV`$. The gain in energy in strange matter, using the bag model, is only a few $`MeV`$ over $`Fe^{56}`$, of the order of thermonuclear energy release. For the realistic EOS it is much larger $`40MeV`$.

The unusually hard X-ray burster GRO J1744$``$28  or the soft $`\gamma `$ repeaters  require the strange star models, since these stars require an energy release which is large compared to thermonuclear energies. With the bag model that would fit the requirement Eq.(1), it is difficult to explain the hard X-ray bursters or soft $`\gamma `$ repeaters since the energy gain over $`Fe^{56}`$ is only $`3`$ $`MeV`$.

## 3 Mass - Radius

Quasiperiodic oscillations in the 4U 1728$``$34, led to the idea that it is an ReSS  with a mass 1.1$`M_{}`$. The star is 4.3 kpc away.

We elaborate on the point mentioned above, namely that if an SS gains mass due to accretion its radius might not necessarily change. We demonstrate this in Fig.3 by the horizontal tangent line drawn on the ReSS M-R relation. It is possible that the X-rays or $`\gamma `$ rays emitted by a compact object is due to the conversion of the accreting normal matter to strange matter . The energy gain by a baryon on conversion is nearly $`40`$ $`MeV`$ or so, making it a very favourable event. Thus a star of small mass may become heavier with a consequent increase of the radius. However, when it reaches the horizontal line $`\frac{dR}{dM}`$ is zero and its radius does not increase for some time. We have fitted the M-R curve to a polynomial,

$$R(M)=\underset{i=1}{\overset{7}{}}a_i\left(\frac{M}{M_{}}\right)^i$$
(3)
The parameters ai are given in Table 1 for the convenience of possible users. The parameters obtained from the two EOS are given in Table 2.

It is found that some X-ray bursters like the 4U 1728$``$34 do not show any radial expansion whereas others like the KS 1731$``$260 are usually observed only in bursts that exhibit photospheric expansion . We conjecture that the persistently bright X-ray transient KS 1731$``$260 was a low mass star gaining mass and thus a radius expansion was allowed by the TOV stability criterion. After 11.5 years of activity this source was detected in January 2001 with Chandra in quiescence with a luminosity that is comparable to other normal X-ray emitting compact stars . This could also be explained if the star, in these decade of expansion, has reached the same stage as the 4U 1728$``$34.

## 4  Radial oscillations of a relativistic star

Thirty five years ago Chandrasekhar  investigated these radial modes. They give information about the stability of a stellar model. Recently such studies have been carried out for many existing neutron star EOS . Sharma et al  made a stability analysis of Eqn. (2) and found that the EOS is stable against radial oscillations. However, a more detailed analysis is needed to restrict the M-R relation. For completeness an outline of the scheme is given below.

The spherically symmetric metric is given by the line element

$$ds^2=e^{2\nu }dt^2+e^{2\mu }dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2).$$
(4)
Together with the energy-momentum tensor for a perfect fluid, Einstein’s field equations yield the TOV equation which can be solved if we have an EOS, $`p(n_B)`$ and $`ϵ(n_B)`$. Given the central density $`ϵ_c`$, we can arrive at an $`MR`$ curve by solving the TOV. Without disturbing the spherical symmetry of the background we define $`\delta r(r,t)`$, a time dependent radial displacement of a fluid element located at the position $`r`$ in the unperturbed model which assumes a harmonic time dependence, as

$$\delta r(r,t)=u_n(r)e^{i\omega _nt}.$$
(5)
The dynamical equation governing the stellar pulsation in its $`n`$th normal mode ($`n=0`$, is the fundamental mode) has the Sturm-Liouville’s form ( for details, see  )

$$P(r)\frac{d^2u_n(r)}{dr^2}+\frac{dP}{dr}\frac{du_n}{dr}+\left[Q(r)+\omega _n^2W(r)\right]u_n(r)=0,$$
(6)
where $`u_n(r)`$ and $`\omega _n`$ are the amplitude and frequency of the $`n`$th normal mode, respectively. The functions $`P(r)`$, $`Q(r)`$ and $`W(r)`$ are expressed in terms of the equilibrium configuration of the star and are given by

$$P(r)=\frac{\mathrm{\Gamma }p}{r^2}e^{\mu +3\nu }$$
(7)
$$Q(r)=e^{\mu +3\nu }\left[\frac{(p^{})^2}{r^2(ϵ+p)}\frac{4p^{}}{r^3}\frac{8\pi }{r^2}(ϵ+p)pe^{2\mu }\right]$$
(8)
$$W=\frac{(ϵ+p)}{r^2}e^{3\mu +\nu },$$
(9)
where the varying adiabatic index $`\mathrm{\Gamma }`$ is given by

$$\mathrm{\Gamma }=\frac{(ϵ+p)}{p}\frac{dp}{dϵ},$$
(10)
$`ϵ`$ and $`p`$ being the energy density and pressure of the unperturbed model, respectively. Eigenfrequencies can be obtained with the boundary conditions,

1. at the centre $`r=0`$, $`\delta r=0`$ and
2. at the surface $`\delta p=0`$ leading to $`\mathrm{\Gamma }pu(r)^{}=0`$,

Since $`\omega `$ is real for $`\omega ^2>0`$, the solution is oscillatory. However for $`\omega ^2<0`$, the angular frequency $`\omega `$ is imaginary, which corresponds to an exponentially growing solution. This means that for negative values of $`\omega ^2`$ the radial oscillations are unstable. For a compact star the fundamental mode $`\omega _0`$ becomes imaginary at some central density $`ϵ_c`$ less than the critical density $`ϵ_{critical}`$ for which the total mass $`M`$ is a maximum. At $`ϵ_c=ϵ_c^0`$, $`\omega _0`$ vanishes. All higher modes are zero at even higher central densities. Therefore, the star is unstable for central densities greater than $`ϵ_c^0`$. To illustrate, we plot the eigenfrequencies $`\omega _n`$ against $`ϵ_c`$, the central density in Fig. 5. The fundamental frequency $`\omega _0`$ does vanish at some $`ϵ_c^0`$ while the higher modes remain nonzero.

Numerical values of masses, radii, central densities and the corresponding eigenfrequencies $`\omega _0`$, $`\omega _1`$ and $`\omega _2`$ are given in Tables 3 and 4 for our EOS1 (SS1) and EOS3(SS2) respectively. The corresponding linear fits described by Eqn.(2) for SS1 and SS2 are given in Table 5 and 6. Tables 7 and 8 are for the bag model EOS with different parameters.

## 5  Discussions and summary

In summary, evidence for the existence of strange stars have been accumulating. In the present paper we review possible candidates and suggest that the properties peculiar to some of the compact stars can be explained if they are SS. In particular we point out that mass accumulation due to accretion does not lead to an increase in the radius for stars like the 4U 1728$``$34, claimed to be low mass SS from accretion data . If some compact objects are proven to be ReSS then parameterizations of QCD chiral symmetry restoration at high densities for quarks, the smallest particles known, can be achieved with the help of data from some of the heaviest objects in the Universe.

ReSS are stable against radial oscillations close to the maximum attainable mass. For example, the EOS of SS1 sustains gravitationally, $`M_{max}1.4M_{}`$, R=7 km with a central number density $`n_c16n_0`$. However, the fundamental frequency of radial oscillations becomes zero at around $`n_c9.5n_0`$, destabilizing the star after M=1.36 $`M_{}`$ with R= 7.24 km (Table 3). It is still on the $`\frac{dM}{dR}>0`$ region. Thus the maximum mass star which is stable against radial oscillations has a number density $`9.5n_0`$ at the centre and $`4.7n_0`$ at the surface. Macroscopically, upto this density small vibrations may be sustained.

A corresponding linear fit has the stable values M=1.34 $`M_{}`$ and R=7.35 km(Table 5). Radial oscillation is rather sensitive and the fit is better than 3$`\%`$(Tables 3 and 5). The same pattern is seen for SS2 and the corresponding linear fit (Tables 4 and 6). Thus we see that almost all of the M-R region is stable.

There has been a recent controversy about the star RXJ1856.$``$3754 $``$ whether it can be inferred to be a strange star . Analysis of the observational data of this “no pulsar”, 120 pc away, cool star is not conclusive. As noted by  the stable portion of the $`MR`$ region shown in Fig. 4 can accommodate this star very easily. This could be a possible candidate of our SS model provided its mass and radius are established beyond controversy.

Acknowledgements: The sad demise of our collaborator Dr. Ranjan Ray was painful for us. We dedicate this paper to his memory. MD and MS thank the RRI for hospitality and JD, MS and SB thank the IUCAA.


format: text

text:
# Three-wave mixing of Bogoliubov quasi-particles in a Bose condensate

\[

## Abstract

A dressed basis is used to calculate the dynamics of three-wave mixing between Bogoliubov quasi-particles in a Bose condensate. Due to the observed oscillations between different momenta modes, an energy splitting, analogous to the optical Mollow triplet, appears in the Beliaev damping spectrum of the excitations from the oscillating modes.

\]

Since the experimental realization of Bose-Einstein condensation in trapped atomic gases, which serves as a mono-energetic and dense atomic source, experiments in non-linear atom optics have become feasible. Atomic four-wave mixing (4WM) , Superradiance , and matter-wave amplification ,  are all examples of non-linear atom optics, which involve the mixing of several atomic fields or the mixing of atomic and electromagnetic fields.

In the case of atomic 4WM, bosonic amplification directs one of the products of a collision between an excitation atom and a condensate atom into an, initially, largely populated mode ,,. Thus, a new mode is macroscopically populated through the second collision product. When the energy, $`\epsilon _k`$, of an excitation mode $`𝐤`$ is low compared to the condensate chemical potential, the excitation can no longer be described as a free atom moving with momentum $`𝐤`$, but rather as a collective phonon excitation, which involves a large number of atoms. Atomic 4WM is therefore inadequate for the description of phonon decay. Phonon excitations are described in the framework of Bogoliubov theory . As Bogoliubov excitations are bosons, bosonic amplification will also direct one of the damping products of a phonon into an, initially, largely populated mode, leading to three wave-mixing (3WM) of Bogoliubov quasi-particles.

In this letter we introduce a basis of states which are dressed by the interaction between three, largely populated, modes of Bogoliubov quasi-particles in a Bose-Einstein condensate at zero temperature. Using this basis of quantum states we calculate the time evolution of the system, which exhibits non-linear oscillations between the different momentum modes. The use of the wave-mixing eigen-basis enables the calculation of the system dynamics efficiently and with no approximations. In contrast to other theoretical approaches to wave-mixing, which lead to non-linear differential equations , the dressed state approach turns the problem into a linear one, for which propagation in time is trivial. We show that, due to relative number squeezing, the variance in the number difference between the two low momentum modes remains constant.

In analogy to the treatment of spontaneous photon scattering in the atom-laser dressed system , we treat damping into the quasi-continuum of empty modes, during the 3WM dynamics, as transfer between dressed state manifolds. Thus, the effects of 3WM on the damping process, are calculated. The damping energy spectrum is presented. A transition from elastic to inelastic damping is observed, which, in analogy to the optical Mollow triplet, leads to a splitting of the spectrum into a doublet of resonance energies.

Previously, nonlinear mixing of quasi-particles in a Bose condensate was studied, between the two lowest energy excitation modes in the discrete regime, where the energy of the excited mode is of the order of the energy separation between modes , , .

Our model system is a homogenous condensate of finite volume $`V`$, with $`N_0`$ atoms in the ground state. The Hamiltonian for the system in the Bogoliubov basis, taken to the $`\sqrt{N_0}`$ order, is given by

$$H=H_0+H_{int},$$
(1)
where $`H_0=\frac{1}{2}gnN_0_𝐤\epsilon _kv_k^2+_𝐤\epsilon _kb_𝐤^+b_𝐤`$ is the part of $`H`$ which is diagonalized by the Bogoliubov transformation. $`N_0`$ is the number of atoms in the $`𝐤=0`$ mode. $`\epsilon _k=\sqrt{\epsilon _k^0\left(\epsilon _k^0+2gn\right)}`$ is the recently measured Bogoliubov energy  and $`\epsilon _k^0`$ is the free particle energy $`\epsilon _k^0=\mathrm{}^2k^2/2m`$. $`g=4\pi \mathrm{}^2a/m`$ is the coupling constant, $`a`$ is the s-wave scattering length, $`m`$ is the mass of the atoms. $`b_𝐤^+`$ and $`b_𝐤`$ are the creation and annihilation operators respectively, of quasi-particle excitations with momentum $`𝐤`$. $`v_k`$ and $`u_k`$ are the Bogoliubov quasi-particle amplitudes. This part of the Hamiltonian describes excitations in the condensate with momentum $`𝐤`$ and energy $`\epsilon _k`$.

$`H_{int}`$ is the part of the Hamiltonian which is responsible for the interaction between excitations

$$H_{int}=\frac{g}{2V}\sqrt{N_0}\underset{𝐤,𝐪}{}A_{kq}\left(b_𝐤^+b_𝐪b_{𝐤𝐪}+b_𝐪^+b_{𝐤𝐪}^+b_𝐤\right).$$
(2)
The first term in parentheses is referred to as Landau damping and is analogous to photonic up conversion, such as second harmonic generation in optics , whereas the second term is referred to as Beliaev damping, and is analogous to photonic down conversion . $`A_{kq}`$ is the many-body suppression factor

$`A_{kq}=2u_k\left(u_qu_{kq}v_qu_{kq}u_qv_{kq}\right)`$ (3)
$`2v_k\left(v_qv_{kq}u_qv_{kq}v_qu_{kq}\right).`$ (4)
The main damping mechanism from a single, largely populated mode $`𝐤`$, will be elastic Beliaev damping of the excitations into empty modes, which are on a mono-energetic surface in momentum space .

We now consider the case where the condensate is excited with $`N`$ excitations of momentum $`𝐤`$ and $`M`$ excitations of momentum $`𝐪`$, $`|N_𝐤,M_𝐪,0_{𝐤𝐪}`$, such that $`𝐤`$ and $`𝐪`$ fulfill the Bragg condition, i.e. $`\epsilon _𝐤=\epsilon _𝐪+\epsilon _{𝐤𝐪}`$. Two-photon Bragg transitions can be used to excite the condensate, and populate different momentum modes with a variable number of excitations . The amplitude given by $`H_{int}`$ for Beliaev damping of the $`𝐤`$ momentum excitation into two excitations with momenta $`𝐪`$ and $`𝐤𝐪`$, $`|\left(N1\right)_𝐤,\left(M+1\right)_𝐪,1_{𝐤𝐪}`$, will be $`\sqrt{N}`$ or $`\sqrt{M}`$ fold larger than any of the other damping channels. This will result in 3WM into a newly populated $`𝐤𝐪`$ momentum mode. We assume that $`N_0>>N,M>>\mathrm{\Gamma }t_0`$, where $`\mathrm{\Gamma }`$ is the total Beliaev damping rate of the excitations into the quasi-continuum of empty modes and $`t_0`$ is the time of the experiment. Thus, during the experiment time, 3WM dynamics dominates over damping into empty modes. In general the set of $`N+1`$ excitation Fock states $`|\left(Ni\right)_𝐤,\left(M+i\right)_𝐪,i_{𝐤𝐪}`$, where $`i`$ varies between $`0`$ and $`N`$, spans a degenerate subspace of the eigenstates of $`H_0`$. $`H_{int}`$ couples between pairs of states in this subspace which have a difference of $`1`$ in $`i`$, and can be represented by the $`(N+1)\times (M+1)`$ tri-diagonal matrix

$`Ns,M+s,s\left|H_{int}\right|Ni,M+i,i=`$ (5)
$`{\displaystyle \frac{g}{2V}}\sqrt{N_0}A_{kq}(\sqrt{Ni}\sqrt{M+i+1}\sqrt{i+1}\delta _{i+1,s}+`$ (6)
$`\sqrt{Ni+1}\sqrt{M+i}\sqrt{i}\delta _{i1,s}).`$ (7)
When $`H_{int}`$ is diagonalized, we get a new set of $`N+1`$ eigenstates, $`|j`$, where $`j`$ varies between $`1`$ and $`N+1`$, that are dressed by the interaction. Note that since the dressed states are superpositions of degenerate eigenstates of $`H_0`$, they are eigenstates of the complete Hamiltonian in Eq. (1).

The filled circles in Fig. 1 show the calculated energy spectrum of the dressed basis for $`N=M=100`$ in units of $`(g\sqrt{N_0}/2V)A_{kq}`$. The degeneracy between the excitation Fock states is removed by the interaction. The spectrum appears to be roughly linear with an average energy spacing of $`dE2.48\sqrt{N}`$ in the above units . The hollow circles in Fig. 1 show the value of the energy difference between each two dressed states. Due to the non-linearity of the problem, energy differences between each pair of dressed states are slightly different.

Figure 2 shows the absolute value squared of the transfer matrix between the bare excitation Fock state basis and the dressed basis for $`N=M=100`$. Since $`H_{int}`$, after being applied a sufficient number of times, can couple each two bare states, most dressed states are spanned by a large number of bare states .

In the limit of $`M=N`$ $`\mathrm{}`$ we can use the tri-diagonality of $`H_{int}`$ to approximate the dressed state as the solution to the differential equation

$$\frac{a_j^2\left(x\right)}{x^2}=\left(\frac{\lambda _j}{\sqrt{\left(1x^2\right)x}}f_{0}^{}{}_{}{}^{2}\right)a_j\left(x\right),$$
(8)
with the boundary conditions $`da_j\left(x\right)/dx=\left(\lambda _j/\sqrt{N}N\right)a_j\left(x\right)`$ at $`x=0`$, and $`da_j\left(x\right)/dx=\left(N\lambda _j/\sqrt{2N}\right)a_j\left(x\right)`$ at $`x=1`$. Here $`dx=1/N`$, $`x=idx`$, and $`a_j\left(x\right)=Ni,M+i,i|j`$ is the projection of $`|j`$ on the Fock state with $`i`$ excitations removed from the $`𝐤`$ momentum mode. Also, $`\lambda _j=E_j\sqrt{N}`$, $`E_j`$ is the energy eigenvalue of $`|j`$, and $`f_{0}^{}{}_{}{}^{2}=2N^2`$. The main contribution to the dressed state superposition of Fock states comes from the two states which solve $`\left(\frac{\left|\lambda _j\right|}{\sqrt{\left(1x^2\right)x}}f_{0}^{}{}_{}{}^{2}\right)=0`$. The dashed-dotted line in Fig. 2 draws the solution of this equation. In order for this equation to have real solutions the resulting spectrum must satisfy $`\left|E_j\right|<1.24N^{3/2}`$. Since there are $`N+1`$ dressed states we get an average energy difference between dressed states of $`dE2.48\sqrt{N}`$, consistent with our numerical observation.

We choose as a model system a condensate of $`3\times 10^{5\text{ }87}Rb`$ atoms in the $`F=2`$, $`m_f=2`$ ground state. The condensate, which is similar to the experimental parameters of  is homogeneous and has a density of $`3\times 10^{14}`$ $`atoms/cm^3`$. The $`k=0.7`$ $`\mathrm{}/\xi `$ mode and $`q=k/\sqrt{2}`$ mode, where $`\xi `$ is the healing length of the condensate given by $`\xi =\sqrt{8\pi na}^1`$, are populated with a 100 excitations each.

We start from the above excitation Fock state, written as a linear superposition of dressed states, which are obtained by diagonalizing (4). The state of the system is then readily propagated in time by evolving each of the dressed states phases according to its energy. Figure 3 shows $`N_𝐤`$, the expectation value of the number of excitations with momentum $`𝐤`$, as a function of time. Excitations oscillate between the $`𝐤`$ momentum mode and the $`𝐪`$ and $`𝐤𝐪`$ momenta modes. Inset (a) of Fig. 3 shows $`N_k`$ during a much longer time. Since the energy spectrum is not precisely linear, at longer times beating between the different oscillation frequencies gives rise to a slow amplitude modulation of the oscillations.

Even though we start from an excitation Fock state, the system immediately evolves into a superposition of Fock states. Inset (b) of Fig. 3 shows the standard deviation of $`N_𝐤`$ vs. time. After a time scale which is set by the non-linearity of the spectrum, the average value of $`N_𝐤`$ and its standard deviation are of similar size. However, the expectation value of the number difference between mode $`𝐪`$ and mode $`𝐤𝐪`$ equals $`M`$ and is constant in time. The standard deviation of that difference always remains zero, which implies relative number squeezing between the two modes.

Thus far we have discussed the time evolution of the system within the dressed state manifold, which is defined by the initial population of the $`𝐤`$ and $`𝐪`$ modes, and have neglected damping into empty modes. In analogy to the treatment of spontaneous photon scattering as transfer between dressed states manifolds of the dressed atom-laser system , we now consider scattering into empty modes as transfer between dressed state manifolds. Thus, the damping of excitations from mode $`𝐤`$ is no longer elastic, but rather carries the energy difference between the dressed states among which it occurred.

Evaluating $`H_B=\frac{g}{2V}\sqrt{N_0}_𝐪^{}A_{kq^{}}\left(b_𝐪^{}^+b_{𝐤𝐪^{}}^+b_𝐤\right)`$, between every pair of states in the two manifolds reveals the spectral structure of the damping process. We find that, the spectrum of the $`N^{}=N1`$, $`M^{}=N`$ manifold, which has one less energy eigenvalue, is shifted by roughly half the energy difference with respect to the $`N=M`$ energy spectrum. $`H_B`$ significantly couples only dressed states with neighboring energies. This results in a structure of a doublet in the Beliaev damping spectrum .

The finite lifetime of the dressed state results from the fact that there is a quasi-continuum of $`N^{}=N1`$, $`M^{}=N`$ manifolds, all with an identical energy spectrum, to which a dressed state in the $`N=M`$ manifold can couple to via $`H_B`$. To determine the width of each transition between two dressed states, we use the Fermi golden rule. The damping rate between a state $`|j_{N,M}`$ in the $`N`$, $`M`$ manifold and a state $`|i_{N1,M}`$ in the $`N1`$, $`M`$ manifold is then given by

$`\mathrm{\Gamma }={\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{𝐪^{}}{}}{\displaystyle \frac{g^2N_0}{2V^2}}\left|A_{kq^{}}|^2\right|{}_{N1,M}{}^{}i\left|b_𝐤\right|j_{N,M}^{}|^2\times `$ (9)
$`\delta \left(\epsilon _𝐤+\epsilon _0\epsilon _𝐪^{}\epsilon _{𝐤𝐪^{}}\right),`$ (10)
where $`\epsilon _0`$ is the energy difference between $`|j_{N,M}`$ and $`|i_{N1,M}`$. Using momentum conservation the energy conservation $`\delta `$-function becomes a geometrical condition on the angle $`\theta `$ between $`𝐤`$ and $`𝐪^{}`$

$`\mathrm{cos}\left(\theta \right)=`$ (11)
$`{\displaystyle \frac{1}{2kq^{}}}\left[k^2+q^2+1\sqrt{\left(\epsilon _k+\epsilon _0\epsilon _q^{}\right)^2+1}\right],`$ (12)
where momentum is in units of $`\mathrm{}/\xi `$, and energy is in units of $`gn`$. Equation (6) is then simplified to

$`\mathrm{\Gamma }={\displaystyle \frac{g}{2\xi ^3h}}{\displaystyle }\left|A_{kq^{}}|^2\right|{}_{N1,M}{}^{}i\left|b_𝐤\right|j_{N,M}^{}|^2\times `$ (13)
$`{\displaystyle \frac{q^{}}{2k}}{\displaystyle \frac{\left(\epsilon _k+\epsilon _0\epsilon _q^{}\right)}{\sqrt{\left(\epsilon _k+\epsilon _0\epsilon _q^{}\right)^2+1}}}dq^{}.`$ (14)
The spectrum of each transition is taken as a Lorenzian with a width of $`\mathrm{\Gamma }/2\pi N_𝐤`$ and a normalization of $`\mathrm{\Gamma }/N_𝐤`$ around $`ϵ_0`$. $`N_𝐤`$ is the average occupation of mode $`𝐤`$ for the two dressed states involved. Averaging over all of the possible transitions between the manifolds, Fig. 4 shows the damping spectrum between the $`N=M=5\times 10^3`$ and the $`N=5\times 10^31`$, $`M=5\times 10^3`$ manifolds, for the same model system as in Fig. 3 for $`k`$=3.2, 1.6, and 0.7 $`\mathrm{}/\xi `$, and $`q=k/\sqrt{2}`$ . A clear doublet structure is evident. The inset of Fig. 3 shows the energy-conserving surfaces for the two center energies of the $`k`$=0.7 $`\mathrm{}/\xi `$ curve (solid line), and the energy conserving surface for elastic damping from the same mode (dashed line).

The separation between the peaks in Fig. 4, which is equal to the oscillation frequency, decreases with $`k`$ as $`A_{kq}`$. The decrease in the width of each peak as a function of $`k`$ comes from two contributions. Firstly, since damping into empty modes is incoherent, it decreases as $`\left|A_{kq^{}}\right|^2`$. Secondly, as $`k`$ decreases there are less allowed empty modes on the energy-conserving surface. Therefore, the doublet structure is more resolved for the lower $`k`$ values. The difference in width between the positive and negative energy peaks is mainly due to the difference in the number of allowed empty modes on the energy conserving surfaces. For small enough $`k`$ values, with respect to the oscillation frequency, the energy conserving surface for the negative energy peak completely disappears and only the positive energy peak remains. Another contribution to the width of the resonance in Fig. 4 is due to the non-linearity of the dressed state spectrum, which gives a slightly different energy for transitions between different dressed state pairs.

Several other mechanisms contribute to the broadening of the two peaks which are not included in Fig. 4. The fact that only the first scattering event occurs between the $`N=M=5\times 10^3`$ and the $`N=5\times 10^31`$, $`M=5\times 10^3`$ manifolds will further broaden the resonances. Since the energy splitting scales as $`\sqrt{N}`$, in an experiment where one scatters $`dN`$ atoms from mode $`𝐤`$, this will result in a relative broadening of $`\frac{dN}{2N}`$. According to the same scaling, an initial coherent, rather than Fock, state will cause a relative broadening of $`\frac{\sqrt{N}}{2N}`$. The condensate finite size or inhomogeneous density profile, will further contribute to the width of the resonance. We estimate that for the experimental parameters of  a doublet structure in the Beliaev damping spectrum can be resolved. Experimentally, the energy doublet can be observed by computerized tomography analysis of time of flight absorption images of the 3WM system .

Another process which will transfer the system between manifolds is that of a two-photon Bragg transition from the condensate to either the $`𝐤`$ or $`𝐪`$ modes . The line-shape of the Bragg transition into the $`𝐤`$ mode also splits into a doublet structure .

The dressed state formalism can be applied to other processes which involve bosonic amplification. We have applied this method to a strongly Bragg-driven condensate, with large momentum transfer, by diagonalizing $`H_{int}=\frac{h\mathrm{\Omega }_R}{2}\left(a_𝐤^+a_0+a_0^+a_𝐤\right)`$, $`\mathrm{\Omega }_R`$ is the two-photon Rabi frequency, and to atomic 4WM with large relative momentum by diagonalizing $`H_{int}=\frac{g}{2V}\left(a_𝐤^+a_0^+a_{𝐤𝐪}a_𝐪+a_{𝐤𝐪}^+a_𝐪^+a_𝐤a_0\right)`$. Both systems show a doublet structure in the spectrum of collisions from a certain momentum mode, due to rapid oscillations in the population of that mode .

In conclusion, we calculate the wave-mixing dynamics between three, low $`k`$, Bogoliubov quasi-particles in a Bose condensate. The Hamiltonian of this system is diagonalized to the next order in $`\sqrt{N_0}`$. The resulting basis of dressed states allows for the efficient, linear, propagation of the system in time. Non-linear oscillations between the different momentum modes are observed. Relative number squeezing between the $`𝐪`$ and the $`𝐤𝐪`$ momentum modes is shown. Beliaev damping of excitations from these modes is treated as a transfer between dressed states manifolds. The damping process is shown to become inelastic and, similarly to the optical Mollow triplet, exhibits a doublet-like energy spectrum, which is more resolved for the low $`k`$ values.

We thank Pierre Meystre and Eitan Rowen for helpful discussions. This work was supported by the Israeli Science Foundation.


text:
# 1 Introduction

## 1 Introduction

The firmly established disappearance of muon neutrinos of cosmic ray origin  strongly points toward the existence of neutrino oscillations .

The approved first generation long baseline (LBL) experiments — K2K , MINOS , ICARUS  and OPERA  — will search for a conclusive and unambiguous signature of the oscillation mechanism. They will provide the first precise measurements of the parameters governing the main muon disappearance mechanism. In particular, the CERN-NGS beam, specifically optimized for tau appearance, will allow to directly confirm the hints for neutrino flavor oscillation.

In addition to the dominant $`\nu _\mu \nu _\tau `$ oscillation, it is possible that a sub-leading transition involving electron-neutrinos occur as well. In the “standard interpretation” of the 3-neutrino mixing, the $`\nu _\mu \nu _e`$ oscillations at the $`\mathrm{\Delta }m^22.5\times 10^3\mathrm{eV}^2`$ indicated by atmospheric neutrinos is driven by the so-called $`\theta _{13}`$ angle. Indeed, given the flavor eigenstates $`\nu _\alpha (\alpha =e,\mu ,\tau )`$ related to the mass eigenstates $`\nu _i^{}(i=1,2,3)`$ where $`\nu _\alpha =U_{\alpha i}\nu _i^{}`$, the mixing matrix $`U`$ is parameterized as:

$$U(\theta _{12},\theta _{13},\theta _{23},\delta )=\left(\begin{array}{ccc}c_{12}c_{13}& s_{12}c_{13}& s_{13}e^{i\delta }\\ s_{12}c_{23}c_{12}s_{13}s_{23}e^{i\delta }& c_{12}c_{23}s_{12}s_{13}s_{23}e^{i\delta }& c_{13}s_{23}\\ s_{12}s_{23}c_{12}s_{13}c_{23}e^{i\delta }& c_{12}s_{23}s_{12}s_{13}c_{23}e^{i\delta }& c_{13}c_{23}\end{array}\right)$$
(1)
with $`s_{ij}=\mathrm{sin}\theta _{ij}`$ and $`c_{ij}=\mathrm{cos}\theta _{ij}`$.

The best sensitivity for this oscillation is expected for ICARUS at the CERN-NGS. Limited by the CNGS beam statistics at low energy, this search should allow to improve by roughly a factor 5 (see Ref.) the CHOOZ limit on the $`\theta _{13}`$ angle for $`\mathrm{\Delta }m^23\times 10^3\mathrm{eV}^2`$. Beyond this program, new methods will be required in order to improve significantly the sensitivity.

At present, the only well established proposal in this direction is the JHF-Kamioka project. In its first phase, 5 years of operation with the Super-K detector, it aims to a factor 20 improvement over the CHOOZ limit.

In Ref., we have studied an optimization of the CNGS optics that would allow to increase the neutrino flux yield at low energy by a factor 5 compared to the baseline $`\tau `$-optimization of the CNGS beam. This would yield an improvement in the sensitivity by about a factor two, or equivalently an improvement of a factor 10 compared to CHOOZ.

In this paper, FLUKA Monte Carlo simulations are employed to perform a systematic study of particle production and neutrino yields for different beam energy and baselines. Focusing systems adapted to the low energy range are also investigated, to obtain realistic estimates of the achievable neutrino rates.

The whole procedure is assumed to be detector and accelerator independent. Nevertheless, to give a first estimate of the needed beam intensities, the ICARUS T3000 detector size has been assumed as a reference. This is also justified by a recent neutrino detector comparison, following which the 2.35 kton LAr ICARUS fiducial mass is equivalent to a detector of approximately 20 kton of steel or 50 kton of water. For this reason, we concentrate on the intrinsic electron-neutrino background from the beam and do not explicitly calculate other sources of backgrounds that are strongly related to detector performances, such as backgrounds from $`\pi ^0`$ production in neutral current events . We assume that in reality the “best” beam would be complemented by the “best” detector and that these backgrounds will only introduce a small correction to our sensitivity estimates. In addition, we are primarily interested at this stage in the comparison among proton drivers.

We consider so-called “conventional” neutrino superbeams, in which neutrinos are produced by the decay of secondary pions obtained in high-energy collisions of protons on an appropriate target and followed by a magnetic focusing system.

In this kind of beams, the neutrino beam spectrum and its flux are essentially determined by three parameters that can be optimized appropriately:

* the primary proton energy $`E_p`$,
* the number of protons on target $`N_{pot}`$ per year,
* the focusing system, which focus a fraction of the secondary charged pions and kaons (positive or both signs depending on the focusing device).

In order to simplify the problem, we consider three “classes” of proton energies:

* Low energy: protons in a range of a few GeV. We take as reference the CERN-SPL proton driver design with 2.2 GeV kinetic energy and an “upgraded SPL” with similar characteristics but with 4.4 GeV protons;
* Medium energy: we take 20 GeV proton energy, similar to the CERN PS machine and the 50 GeV of the JHF facility;
* High-energy: we take the highest energies, i.e. the 400 GeV protons like in the case of the CERN SPS.

The purpose of this work is to understand the required intensities for the various proton energies, in other words, is it more favorable to employ low-energy, high-intensity proton or high-energy, low-intensity proton machines?

## 2 Observing the oscillation – choosing $`L`$ and $`E`$

In order to maximize the probability of an oscillation, we must choose the energy of the neutrino $`E_{max}`$ and the baseline $`L`$ such that

$$1.27\frac{L(km)}{E_{max}(GeV)}\mathrm{\Delta }m^2(eV^2)\frac{\pi }{2}$$
(2)
However, in order to observe the oscillation, we must at least see the maximum preceded by a minimum, given by

$$1.27\frac{L(km)}{E_{min}(GeV)}\mathrm{\Delta }m^2(eV^2)\pi $$
(3)
The mass difference indicated by Superkamiokande is given by $`\mathrm{\Delta }m^22.5\times 10^3\mathrm{eV}^2`$ and lies within the range $`1\times 10^3<\mathrm{\Delta }m^2<4\times 10^3\mathrm{eV}^2`$.

Since we are considering one single detector at a given location as target, our baseline $`L`$ for a given source is fixed. This implies that the neutrino beam spectrum should not be too narrow, e.g. it should be a relatively wide-band beam in order to cope with the uncertainty on the $`\mathrm{\Delta }m^2`$. Since the knowledge on this parameter will improve with time, in particular with the running of K2K, MINOS and CNGS experiments, the range of energies will eventually be limited by the desire to observe the minimum and the maximum. In the meantime, we consider the full mass range indicated by atmospheric neutrino observation.

We do not wish to consider baselines longer than the current CNGS baseline of $`L=730\mathrm{km}`$ and use therefore a range of baselines between 100 km and 730 km. Table 1 shows the values (in MeV) of $`E_{min}`$ and $`E_{max}`$ defined above as a function for the baselines. The energies range from about 50 MeV to 300 MeV for $`L=100\mathrm{km}`$ and from 300 MeV to 2400 MeV for $`L=730\mathrm{km}`$. Clearly, we wish to optimize for high intensity, low energy neutrino superbeams.

We note that at these low energies, neutrino interactions are clearly identifiable and have generally easily reconstructible final-states. This is an advantage for detector-related background suppression. In a detector with the granularity like ICARUS, neutral pions can be easily suppressed via energy ionization or by direct reconstruction of the two (well separated) decay photons.

## 3 The approximate scaling with proton energy

The understanding of the relationship between primary beam energy and neutrino production is the first step toward any optimization. To isolate this relationship from the other parameters, it is more convenient to work in the “perfect focusing” approximation, where all (positively charged) mesons produced within a given solid angle are supposed to be focused exactly on the detector direction111In order to derive the scaling, we consider only neutrino fluxes. However, in Section 6, when we compute our oscillation sensitivities we shall include both neutrinos and antineutrinos components obtained from the full simulation of the mesons of the relevant charge.. In our case, we set $`P_T=0`$ for all secondary positively charged particles produced within 1 rad.

All the simulations presented in this section have been performed with the same target and decay tunnel geometry. The target is in graphite, density 1.8 g/cm3, one meter long, 2 mm radius; the decay tunnel is 150 m long and 3.5 m in radius. Target optimization will not change the general scaling, and the effect of the tunnel length on muon neutrino production can be accounted for in first approximation assuming that all neutrinos are produced by pions, and that they carry the maximum possible momentum.

The pion production rate at the exit of the target for various incident proton energies $`E_p`$ as estimated with FLUKA are shown in Figure 1 for energies ranging from 2.2 GeV up to the 400 GeV of the SPS for the CNGS beam. In order to compare pion productions at different proton energies, we divide the spectra by the proton energy $`E_p`$.

All normalized spectra have similar shape, with the maximum yield at low energies ($`p_\pi 500`$ MeV). Far from the endpoint, all spectra scale approximately with the incoming proton energy. Departures from the overall scaling consist in a slightly different shape at low $`E_p`$, and harder spectra at high $`E_p`$.

If the detector is far away in the forward direction, the neutrino event rate as a function of neutrino energy can be derived from these spectra by considering that:

* The neutrino carries a momentum that is 0.43 times the parent pion momentum
* The Lorentz boost gives a factor proportional to $`E_\pi ^2E_\nu ^2`$ on the solid angle
* The neutrino cross section grows approximately as $`E_\nu `$ (not true at the lowest energies, where the quasi-elastic grows more rapidly, however, for the optimization of the proton energy, this approximation is adequate.)

For a given Eν, these factors apply independently on the primary proton energy. Thus, we expect that the energy scaling observed on the pion production is reflected in the neutrino production. This is shown in Figure 2, where we plot the simulated muon neutrino event rate scaled with the primary proton energy.

The super-position of the curves at the lowest energies is impressive, except for one at $`E_p=400\mathrm{GeV}`$. Of course, by rising $`E_p`$ the energy-integrated event rate raises dramatically because the spectra extend to higher energies, but the event rate for a given $`E_\nu <0.15\times E_p`$ is simply proportional to $`E_p`$.

This approximate scaling implies that we can define a power factor $``$ for neutrino production as the product of the proton energy times the beam intensity:

$$=E_p\times N_{pot}$$
(4)
Up to now we did not take into account the neutrino oscillation probability. If the aim is that of having a neutrino beam centered on the oscillation maximum, $`E_\nu `$ and the baseline $`L`$ have to be chosen such as

$$\mathrm{sin}^2\left(1.27\mathrm{\Delta }m^2\frac{L}{E_\nu }\right)1,$$
(5)
thus

$$E_{max}L.$$
(6)
How to choose the baseline $`L`$ to maximize the rate of oscillated neutrino events per proton incident on the neutrino target? The neutrino flux grows scales like $`1/L^2`$, like the solid angle. But $`L`$ and $`E_{max}`$ are proportional, thus dividing bin by bin the simulated neutrino spectra with a factor $`E_\nu ^2`$ one gets the shape of the event rate as a function of $`E_{max}`$. This has been done for the perfect focusing approximation in Figure 3 for various proton energies $`E_p`$.

It is evident that the most optimal situation is for $`E_{max}`$ in the range 200-600 MeV and this results holds essentially independently of the proton energy, at least as long as we discard the very high 400 GeV energy case. We can go beyond the $`1/E_\nu ^2`$ approximation and in the following we will use the exact expression for the oscillation probability and will compute exactly the number of oscillated events by folding the expected neutrino spectra.

Nonetheless, the situation gives us a large freedom in the choice of proton energy, provided that the intensity can be re-scaled accordingly, so that the power factor remains essentially constant

$$=E_p\times N_{pot}const.$$
(7)
## 4 The optimal baseline $`L`$ for each proton energy

In order to study the scaling with the proton energy $`E_p`$ and the baseline $`L`$, we compute the number $`N_e`$ of oscillated $`\nu _\mu \nu _e`$ events:

$$N_e_{pot}_0^{\mathrm{}}𝑑E_\nu \sigma (E_\nu )\varphi (E_\nu )P(\nu _\mu \nu _e,E_\nu ,L,\mathrm{\Delta }m^2,\mathrm{sin}^22\theta _{13})$$
(8)
where $`\varphi (E_\nu )`$ is the neutrino flux per proton on target. We normalize this number to a 2.35 kton (argon) detector, and assume 100% electron detection efficiency. For definiteness, we assume unless otherwise noted that $`\mathrm{\Delta }m^2=3\times 10^3\mathrm{eV}^2`$ and $`\mathrm{sin}^22\theta _{13}=1\times 10^3`$.

In order to take into account the effect of focusing (realistic focusing is discussed in section 5), we focus ideally all particles with the acceptance of 1 rad and apply a constant “focusing efficiency” of 20%, i.e. scale down the rates by a factor 5. We will see that this assumption is quite realistic and at least conservative for the low energy neutrinos.

Similarly, one can also estimate the “goodness” of the neutrino flux by computing the number of muon charged current events $`N_{\mu ,CC}^0`$ within the energies $`E_{min}(\mathrm{\Delta }m^2=1\times 10^3\mathrm{eV}^2)`$ and $`E_{max}(\mathrm{\Delta }m^2=4\times 10^3\mathrm{eV}^2)`$ (see Table 1) of the oscillation:

$$N_{\mu ,CC}^0_{E_{min}}^{E_{max}}𝑑E_\nu \sigma (E_\nu )\varphi (E_\nu )$$
(9)
We choose to normalize $`N_{\mu ,CC}^0`$ to $`10^{19}`$ pots and per kton of (argon) target.

With these two parameters, we can easily compare all the possible options. It is fair to note that the number of oscillated events $`N_e`$ assumes that the search is not limited by background, which in practice is not the case, as we will show in section 6. However, we prefer to separate the issue of background (which will be discussed in section 7) and concentrate for the moment on the question of the needed proton intensity. As we will show, the background turns out to be also quite independent of the proton energy and hence it does not enter in the optimization of the proton energy.

The results as a function of the baseline $`L`$ are summarized in Table 2. The configuration corresponding to the minimum $`_{pot}`$ is shown in bold.

We observe that for each baseline there is an optimal proton energy $`E_p^{optimal}`$, which minimizes the required integrated proton intensity $`_{pot}`$ to observe a fixed number of oscillated events. This is also visible in Figure 8, where we plot the integrated beam intensity needed to obtain 5 oscillated $`\nu _e`$ events in a 2.35 kton detector for $`\mathrm{\Delta }m^2=3\times 10^3\mathrm{eV}^2`$ and $`\mathrm{sin}^2(2\theta _{13})=10^3`$ for different beam energies and baselines. The point of minimum pots corresponds approximately to the maximum $`\nu _\mu `$ event rate.

This is because in that point, the secondary pion yield energy spectrum per proton on target, is best matched to the neutrino oscillation probability. Conversely, for each proton energy there is an optimal baseline $`L^{opt}`$, which maximizes the integrated neutrino oscillation probability in the neutrino energy region which corresponds to the largest weighted pion yield at that proton energy.

For a 2.2 GeV proton driver, the optimal baseline $`L^{opt}`$ is approximately $`L^{opt}150\mathrm{km}`$. For a 4.4 GeV proton driver, it is approximately $`L^{opt}200\mathrm{km}`$. For 20 GeV, we find $`L^{opt}450\mathrm{km}`$. For energies above 50 GeV, the optimal baseline is around 700 km.

For actual baselines smaller than the optimal baseline, $`L<L^{opt}`$, the neutrino oscillation maximum occurs at lower energy and the yield of corresponding pions for the given proton energy is lower than in the optimal case. Hence, we need a higher intensity to compensate for this effect.

For actual baselines greater than the optimal baseline, $`L>L^{opt}`$, the neutrino oscillation maximum occurs at higher energy and the yield of corresponding pions for the given proton energy is lower than in the optimal case. Indeed, at some point, the optimal neutrino oscillation energy corresponds to a pion energy, which is kinematically forbidden for the given incident proton energy. Hence, we need again a higher proton intensity to compensate for the kinematical suppression.

At the optimal baselines $`L_{opt}`$ the power factors are the following:

$`_{2.2}2.2\times 1.9\times 10^{24}=4.2\times 10^{24}\mathrm{GeV}\times \mathrm{pot},`$ (10)
$`_{4.4}4.4\times 8.2\times 10^{23}=3.6\times 10^{24}\mathrm{GeV}\times \mathrm{pot},`$ (11)
$`_{20}20\times 1.0\times 10^{23}=2.0\times 10^{24}\mathrm{GeV}\times \mathrm{pot},`$ (12)
$`_{50}50\times 5.2\times 10^{22}=2.6\times 10^{24}\mathrm{GeV}\times \mathrm{pot},`$ (13)
$`_{400}400\times 1.3\times 10^{22}=5.2\times 10^{24}\mathrm{GeV}\times \mathrm{pot}.`$ (14)
This directly confirms the approximate scaling with $``$, apart at the lowest and the highest proton energies. Strictly speaking, a proton energy of 20 GeV appears to be the most economical choice in terms of protons, with a proton economy of about a factor 2 compared to the 2.2 and 400 GeV cases.

Summarizing, the required protons on target for 5 oscillated events at $`\mathrm{\Delta }m^2=3\times 10^3\mathrm{eV}^2`$ and $`\mathrm{sin}^22\theta _{13}=10^3`$ in the 2.35 kton mass located at the optimal baseline assuming a 20% ideal focusing are $`2\times 10^{24}`$ for 2.2 GeV protons, $`10^{24}`$ for 4.4 GeV, $`10^{23}`$ for 20 GeV, $`5\times 10^{22}`$ for 50 GeV and $`10^{22}`$ for 400 GeV.

## 5 Results in real focusing

The standard focusing system for all neutrino beams up to now is based on magnetic horns. However, new solutions should be envisaged for low-energy, intense beams. The angular acceptance of the system has to be large (the average transverse momentum of reaction products is of the order of 200-300 MeV, comparable with the total momentum for the shortest baselines), and the amount of material within the secondary beam cone must be as little as possible, both to preserve the flux and to avoid heating/damage at high intensities.

In principle, the focusing system should cancel the transverse momentum of all secondary particles relative to the direction of flight toward the detector and this independently of their momentum $`p`$. Indeed, this is the definition of the ideal focusing.

In practice, a horn system can be designed to focalize a given signed momentum, i.e. like $`p`$. A FODO system focalizes like $`|p|`$. A series of coils should focalize like $`1`$, i.e. independent of $`p`$. Hence, this last solution seems to be attractive in our situation.

We have found with the help of full simulations and tracking of the particles in the focusing system that (1) for the lowest energy configuration, a series of coils does indeed provide a quite optimal focusing and (2) for the higher energy configurations, where the coil focusing becomes impracticable, the traditional horn focusing can be efficiently used.

### 5.1 Details on coil focusing

One possible solution is to exploit the focusing capabilities of magnetic field gradients, such as the fringe field at the ends of a solenoid, or the decreasing field far from a current loop. The fact that a particle traveling almost parallel to the axis of a solenoid suffers a change of its transverse momentum when traversing the fringe region is not intuitive but can be understood by realizing that any change of the longitudinal component of an axially symmetric field is associated through Maxwell’s divergence equation to a radial component of the field.

The particle motion in axially symmetric fields can be described on the basis of the Busch’s theorem, which in turn follows from the conservation of canonical angular momentum. The derivation and applications of this theorem can be found in textbooks (see for instance ), here we give a summary to better illustrate our focusing system. Suppose we have a particle having electric charge $`q`$ and total energy $`E`$, moving in a static, non-uniform axially symmetric field, at a distance $`r`$ from the field axis ($`z`$axis), with an angular velocity $`\dot{\varphi }`$. Busch’s theorem states that

$$rp_\varphi +\frac{q}{2\pi }\mathrm{\Psi }=const$$
(15)
where $`p_\varphi =\gamma m\dot{\varphi }`$ is the azimuthal component of the particle momentum, and $`\mathrm{\Psi }`$ is the magnetic flux linked by a circle of radius $`r`$ centered on the field axis.

If the field can be considered constant within the circle of radius $`r`$, we have

$$\mathrm{\Psi }=\pi r^2B_z,\mathrm{thus}rp_\varphi +\frac{q}{2}r^2B_z=const$$
(16)
If particles are emitted with zero angular velocity by a source located on the axis of an axially symmetric magnetic field, both the angular momentum and $`\mathrm{\Psi }`$ are null, thus $`const=0`$.

When exiting to a region of zero field, according to Busch’s theorem the particles will have zero angular momentum. This is very appealing remembering that in an uniform (solenoidal) field, all the transverse momentum is azimuthal. However, this is not the case when the field varies, and part of the momentum can also be directed along the radius. The radial motion can also be derived exploiting Busch’s theorem, but the practical use is not straightforward.

However, if the variation of the magnetic field is small on the scale of the revolution time, the theorem of adiabatic invariance states that the magnetic flux encircled by the particle trajectory remains a constant of motion. The trajectory radius can be assumed constant over a revolution and given by the usual expression

$$R=\frac{p_T}{qB_z},$$

and from adiabatic invariance follows

$$B_zR^2\pi =\frac{\pi p_T^2}{q^2B_z}=const.$$

Under this condition, all the rotational motion is converted into the longitudinal one, and the particle transverse momentum varies as

$$\frac{p_T^2}{B_z}=\frac{p_{T0}^2}{B_{z0}}$$
(17)
where $`B_{z0},p_{T0}^2`$ are the initial magnetic field and transverse momentum.

Since the period of the motion is

$$\tau =\frac{2\pi }{\omega _B}=\frac{2\pi E}{qB_z},$$

the condition of adiabatic motion means that $`\mathrm{\Delta }B_z`$ must be small over distances

$$\mathrm{\Delta }z=\frac{2\pi Ev_{}}{qB_z}=\frac{2\pi p_{}}{qB_z}.$$

For a pion having a longitudinal momentum of 1 GeV/c, in B=1 T, $`\mathrm{\Delta }z`$ is about 20 m. Is is thus difficult to focus efficiently high energy particles through fringe fields, while this method can be very effective for low energy ones. It should be stressed that this method applies equally well to positive and negative particles, even though the resulting advantage is not great since the positive component is dominant for low energy beams.

### 5.2 Results on horn and coil focusing

We report in Table 3 our results with full detailed simulations of the focusing systems for 2.2 GeV 4.4 GeV, 20 GeV and 400 GeV proton energies.

The coil method has been applied here to the 2.2 and 4.4 GeV proton beams, where the produced pions have small energies (only positive mesons have been considered in Table 3). The target has been assumed to be a short (30 cm) mercury target, like in the SPL proposal. The non-uniform magnetic field has been obtained with ten circular loops, having a radius of 1 m, positioned from 0 to 14 meters from the target, carrying decreasing currents to give a field from 20 T to zero. Examples of particle orbits and the central magnetic field intensity are shown in Figure 4. The effect on the transverse momentum distribution can be appreciated from Figure 5, where positive pions emitted from the target within 1 radian have been considered. The decay tunnel is 150 m long.

For higher energy neutrino beams of 20 GeV and 400 GeV, the traditional two-horns system has been used. The calculations presented here refer to a first horn to focus 2 GeV/c particles, followed by a reflector to focus 3 GeV/c. The horn is placed around a graphite target, has a length of 4 meters and a current of 300 kA. The reflector starts at 6 m from the target, is 4 m long with a current of 150 kA. Examples of particle trajectories can be seen in Figure 6, and the effect on the transverse momentum in the case of a 20 GeV primary beam is shown in Figure 7. The decay tunnel is 350 m long.

In the total energy range (0-2.5 GeV), the resulting focusing efficiency varies in between 0.5 and 0.2. This motivated our choice in Section 4, when comparing all possible beam/distance options, where a common energy-independent 20% focusing efficiency had been assumed. In reality, the focusing profile is not energy independent, and the real situation can be better than this, as can be derived from the comparison of Tables 2 and 3. This effect is also visible in Figure 8.

Summarizing, the required protons on target for 5 oscillated events at $`\mathrm{\Delta }m^2=3\times 10^3\mathrm{eV}^2`$ and $`\mathrm{sin}^22\theta _{13}=10^3`$ in the 2.35 kton mass located at the optimal baseline assuming optimized realistic focusing are $`1\times 10^{24}`$ for 2.2 GeV protons, $`4\times 10^{23}`$ for 4.4 GeV, $`5\times 10^{22}`$ for 20 GeV and $`5\times 10^{21}`$ for 400 GeV.

The cautious reader can wonder why we have apparently been “conservative” in assuming a 20% efficiency for the ideal focusing. The point here is that in the ideal case we assume a constant efficiency over the whole meson energy range, while in the real focusing case one effectively reaches a situation where a part of the meson energy range is focalized with an efficiency better than 20% while other parts have lower efficiencies. We have of course optimized focusing for the energy relevant to the oscillation. It would however be incorrect to assume that this higher efficiency is constant over the full energy range.

## 6 A superbeam to Gran Sasso (BNGS222BNGS stands for “Better Neutrinos to GS”.) ?

### 6.1 Finding the location for the source

As a working hypothesis, we take for granted the existence of LNGS as an underground laboratory that can host large neutrino detectors. In particular, we have assumed for this study the existence of ICARUS with 2.35 kton fiducial mass.

If we take the location of the detector as fixed, the baseline $`L`$ is determined by the location of the neutrino source. We have investigated various potential locations within Italy where large ENEA444ENEA stands for “Ente per le Nuove Tecnologie, l’Energia e l’Ambiente”. or INFN infrastructures are already existing, as shown in Table 4. In these locations, the required conditions to host a high intensity proton driver could be met and the machine would find other applications in addition to neutrino physics.

In Table 4, the column called “angle” describes the space angle relative to the orientation of the LNGS Halls. Inclination is the incoming neutrino angle relative to the horizontal plane at LNGS. We note that due to the fortunate orientation of the LNGS Halls according to the geographical axis of the Italian peninsula, it is possible to find various laboratories within the orientation of the LNGS Halls.

Indeed, ENEA Aquilone, ENEA Brasimone, ENEA Trisaia, INFN Pavia and ENEA Ispra appear with an angle less than $`15^o`$ with respect to the LNGS Hall direction. This is an advantage for the acceptance of higher energy events given the natural longitudinal orientation of the detectors. On the other hand, in the case of ENEA Casaccia situated near Rome, the angle is $`63^o`$. However, for the shortest baseline, we expect the relevant neutrinos to have an energy similar to that of most atmospheric neutrinos, so that we can argue that the acceptance will not be a problem given the isotropical nature of a detector like ICARUS.

### 6.2 The medium baselines ($`120<L<470\mathrm{km}`$)

We study the cases ENEA Casaccia ($`L=120\mathrm{km}`$), ENEA Brasimone ($`L=270\mathrm{km}`$) and INFN Pavia ($`L=465\mathrm{km}`$). At these distances, the neutrino energy range relevant for $`\nu _\mu \nu _e`$ search is 0.1-1.0 GeV. We assume 2.2, 4.4 and 20 GeV proton energies. Expected neutrino fluxes and rates, obtained assuming ideal (p.f.) and real (coil or horn) focusing systems, are reported in Table 5. The figures are normalized to the baseline $`L=270\mathrm{km}`$ and for $`0E_\nu 2.5\mathrm{GeV}`$. Other baselines can be rescaled accordingly. With real focusing, we find about 100(10) $`\nu _\mu `$($`\overline{\nu }_\mu `$) CC events/kton per $`10^{23}`$ pots for 2.2 GeV protons, about 340(40) $`\nu _\mu `$($`\overline{\nu }_\mu `$) CC events/kton for 4.4 GeV protons and about 26000(250) $`\nu _\mu `$($`\overline{\nu }_\mu `$) CC events/kton for 20 GeV protons. Intrinsic $`\nu _e`$($`\overline{\nu }_e`$) beam contaminations are in the range from 1% to 2% with respect to $`\nu _\mu `$($`\overline{\nu }_\mu `$).

The 2.2 and 4.4 GeV proton energies are the closest to the optimal energies and only these two cases are considered in the following. We will consider 20 GeV in section 6.3 in the context of the CERN-LNGS baseline.

To appreciate the matching between the neutrino beams and the oscillation probability, we show in Figures 9 and 10 the charged current event rate at 120 km as a function of the neutrino energy. The dotted lines correspond to the oscillation probability (arb. norm.) for a $`\mathrm{\Delta }m^2=3\times 10^3\mathrm{eV}^2`$.

Clearly, these neutrino beams offer optimal condition to study $`\nu _\mu `$ disappearance. Indeed, the maximum of the oscillation is very well covered by the neutrino beam and hence it is quite obvious that a very precise determination of the main oscillation parameters will be accomplished. We here do not consider this any further.

We, however, concentrate instead on the $`\nu _e`$ appearance measurement. Since we are in the presence of intrinsic $`\nu _e`$ background from the beam at the level of 1%-2% of the $`\nu _\mu `$ component, we can improve our sensitivity by studying the energy spectrum of the $`\nu _e`$ charged current events. This method is more sensitive than simple event counting.

In order to estimate the sensitivity, we adopt our standard fitting procedure of the various reconstructed event classes (See Ref. ). We assume that the neutrino and antineutrino interactions cannot be distinguished on an event-by-event basis, and hence add the $`\nu _e`$ and $`\overline{\nu }_e`$ contributions from the beam. Similarly, the oscillated spectrum is calculated by summing both $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations, assuming the same oscillation probability for neutrino and antineutrinos.

In the present study, we considered only the energy distribution of electron events and computed the $`\chi ^2`$ as a function of the $`\mathrm{sin}^22\theta _{13}`$ mixing angle, scanning in $`\mathrm{\Delta }m^2`$. The 90% C.L. sensivity region is defined by the condition $`\chi ^2>\chi _{min}^2+4.6`$, defined by the condition that the actually observed events in the experiment coincide with the expected background. In an actual experiment, a simultaneous fit of the muon disappearance and electron appearance spectra will constrain the $`\mathrm{\Delta }m^2`$, $`\mathrm{sin}^2\theta _{23}`$ parameters and in case of negative result will limit the $`\mathrm{sin}^22\theta _{13}`$ within the allowed $`\mathrm{\Delta }m^2`$ region.

The results of the $`\mathrm{\Delta }m^2`$ scans are shown in Figure 11 for 2.2 GeV and Figure 12 for 4.4 GeV proton energy for the three assumed baselines. The curves correspond to 5 years running with $`2\times 10^{23}`$ pots/year equivalent to a continous proton current of 1 mA and a fiducial mass of 2.35 kton. The assumed protons on target is compatible with an accelerator with performances similar to those of the planned CERN-SPL.

The expected sensitivities represent a great improvement to the CHOOZ limit which gives for $`\mathrm{\Delta }m^2=2.5\times 10^3\mathrm{eV}^2`$,

$`(\mathrm{sin}^22\theta _{13})_{CHOOZ}<0.14\mathrm{or}\theta _{13}<11^\mathrm{o}`$ (18)
Indeed, we find for $`2.2\mathrm{GeV}`$ protons:

$`(\mathrm{sin}^22\theta _{13})_{BNGS,120km}`$ $`<`$ $`0.006\mathrm{or}\theta _{13}<2.2^\mathrm{o}`$ (19)
$`(\mathrm{sin}^22\theta _{13})_{BNGS,270km}`$ $`<`$ $`0.015\mathrm{or}\theta _{13}<3.4^\mathrm{o}`$ (20)
$`(\mathrm{sin}^22\theta _{13})_{BNGS,465km}`$ $`<`$ $`0.03\mathrm{or}\theta _{13}<5^\mathrm{o}`$ (21)
and for $`4.4\mathrm{GeV}`$ protons:

$`(\mathrm{sin}^22\theta _{13})_{BNGS,120km}`$ $`<`$ $`0.0035\mathrm{or}\theta _{13}<1.7^\mathrm{o}`$ (22)
$`(\mathrm{sin}^22\theta _{13})_{BNGS,270km}`$ $`<`$ $`0.006\mathrm{or}\theta _{13}<2.2^\mathrm{o}`$ (23)
$`(\mathrm{sin}^22\theta _{13})_{BNGS,465km}`$ $`<`$ $`0.02\mathrm{or}\theta _{13}<4^\mathrm{o}`$ (24)
For comparison, the JHF proposal with the OAB beam gives similar results

$`(\mathrm{sin}^22\theta _{13})_{JHF,OAB}<0.006\mathrm{or}\theta _{13}<2.2^\mathrm{o}.`$ (25)
### 6.3 The CERN-GS baseline

The baseline between CERN and GS is 730 km. At this distance, the neutrino energy range relevant for $`\nu _\mu \nu _e`$ search is 0.3-2.5 GeV.

The present CNGS design is optimized for $`\nu _\tau `$ appearance, thus for a relatively high-energy neutrino beam. The 400 GeV/c SPS beam will deliver 4.5 1019 protons per year on a graphite target, made of spaced rods to reduce the re-interaction rate. The two magnetic horns (horn and reflector) are tuned to focus 35 and 50 GeV/c mesons, with an acceptance of the order of 30 mrad. The decay tunnel length is 1 km. With the standard CNGS parameters, the low-energy neutrino flux is low, as can be seen from the entries flagged by  in Table 6.

In Ref. , we have studied a L.E. optimization of the 400 GeV protons of the CNGS in order to improve the sensitivity to $`\theta _{13}`$. This yielded an improvement of a factor 5 in flux at low energy compared to the $`\tau `$ optimization.

Here, we study the 20 GeV proton energy (we call this the PS++) and compare it to 400 GeV. Expected neutrino fluxes and rates obtained with real focusing systems are reported in Table 6.

The results of the $`\mathrm{\Delta }m^2\mathrm{sin}^22\theta _{13}`$ sensitivity scans are shown in Figure 14. The curves correspond to 5 years running with $`2\times 10^{21}`$ or $`2\times 10^{22}`$ pots/year and a fiducial mass of 2.35 kton.

We find for $`\mathrm{\Delta }m^2=2.5\times 10^3\mathrm{eV}^2`$:

$`(\mathrm{sin}^22\theta _{13})_{PS++,2\times 10^{21}pot/year}<0.016\mathrm{or}\theta _{13}<3.6^\mathrm{o}`$ (26)
$`(\mathrm{sin}^22\theta _{13})_{PS++,2\times 10^{22}pot/year}<0.005\mathrm{or}\theta _{13}<2^\mathrm{o}`$ (27)
Clearly, a strongly intensity-upgraded CERN PS booster would provide very interesting possibilities for the oscillation searches over the CERN-GS baseline.

## 7 The intrinsic $`\nu _e`$ background

As well known, the intrinsic electron (anti)neutrinos in the beam are produced either in the decay of muons coming from kaons or pions via the chain

$$\pi ^+/K^+\mu ^++\nu _\mu ,\mu ^+\nu _e+e^++\overline{\nu }_\mu ,$$

or directly in the three-body $`K_{e3}`$ kaon decays.

When looking for $`\nu _\mu \nu _e`$ oscillations, this contamination will eventually be the limiting factor. It is therefore essential to understand its level and it is also worth understanding if the beam design can be optimized to minimize this background.

Rather than the $`\nu _e/\nu _\mu `$ ratio, we decide to consider the ratio $`\nu _\mu /\sqrt{\nu _e}`$ in order to better estimate the effect of the backrgound on the $`\nu _\mu \nu _e`$ oscillation sensitivity.

Which source of $`\nu _e`$ is relevant to our study?

Due to the large difference between the $`\pi `$ and $`\mu `$ decay lengths, the electron neutrino background depends on the length of the decay tunnel $`l`$. This dependence should be more evident for low beam energies.

Kaon production strongly depends on the proton energy, as shown in Figure 15, where a threshold effect is clearly visible. For the lowest proton energies, 2.2 and 4.4 GeV, the fraction of kaon relative to pions is a few per mille. Above 20 GeV, it is on the order of a little less than 10%.

At low $`\nu _e`$ energies, the production is shared by kaon and muon decays, while kaon decays alone are responsible for the high energy tail. For low proton beam energies, kaon production is much lower and practically all the intrinsic electron contamination comes from muon decays.

A first order estimate of the effect of the decay tunnel length $`l`$ on the muon-induced background can be derived assuming forward decay at each step, and counting for each neutrino energy the fraction of parent particles that decay within a path $`l`$. The fraction $`D_\pi `$ of $`\pi `$ decayed after a length $`l`$ is simply given by

$$D_\pi =(1e^{\frac{l}{\lambda _\pi }})$$
(28)
Muons have to be generated by a pion first, thus

$`D_\mu `$ $`=`$ $`{\displaystyle _0^l}\lambda _\pi ^{}e^{\frac{y}{\lambda _\pi ^{}}}\left(1e^{\frac{ly}{\lambda _\mu }}\right)𝑑y`$ (29)
$`=`$ $`1{\displaystyle \frac{1}{\lambda _\mu \lambda _\pi ^{}}}\left(\lambda _\mu e^{\frac{l}{\lambda _\mu }}\lambda _\pi ^{}e^{\frac{l}{\lambda _\pi ^{}}}\right)`$
The decay lengths $`\lambda _\pi ^{},\lambda _\mu `$ depend on the energy of the meson. We can fix a neutrino energy $`E_\nu `$, the same for $`\nu _\mu `$ and $`\nu _e`$. For $`\nu _\mu `$ production, we can then assume $`p_\pi =E_\nu /0.43`$. For the $`\nu _e`$ produced in the decay of a muon, there is no fixed relation, but we can take the average of the $`\nu _e`$ energy in the muon rest frame, getting $`p_\mu E_\nu /0.6`$. The grand-parent pion had therefore a momentum equal to $`p_\pi ^{}p_\mu E_\nu /0.6`$.

With these approximations, $`D_\pi `$ and $`D_\mu `$ can be expressed as a function of $`l/E_\nu `$. In these units, one $`\pi `$ decay length corresponds to $`l/E_\nu =130`$ m/GeV.

It is a priori obvious that short decay tunnels reduce the relative probability of muon production and decay. However, the pion decay yield is also affected.

To study the background from muons, we consider the statistically significant ratio $`\pi decays/\sqrt{\mu decays}`$ as the correct estimator for $`\nu _\mu /\sqrt{\nu _e}`$ as a function of the decay tunnel length. This ratio is shown in Figure 16, where the assumption of a fixed $`\nu _e`$ energy in $`\mu `$ decay has been relieved: the hatched band corresponds to 20%-80% of the maximum $`\nu _e`$ energy in the $`\mu `$ rest frame.

We find that the ratio does not show dramatic variations between 0 and 4 pion decay lengths. We therefore conclude that not much is to be gained by reducing the length of the decay tunnel.

We have verified these results directly by the full simulation of the neutrino beams for various decay tunnel lengths. Answers are reported in Table 7 for various proton energies and decay tunnel lengths. The 6th column shows the expected $`\nu _e`$ contamination relative to the $`\nu _\mu `$ and the last column lists the statistically relevant ratio $`\nu _\mu /\sqrt{\nu _e}`$. For the 2.2 GeV proton energy, the ratio $`\nu _e/\nu _\mu `$ varies from 0.3% for $`l=20m`$ up to 1.7% for $`l=150m`$, but this happens at a high cost of genuine $`\nu _\mu `$’s. The statistically relevant ratio $`\nu _\mu /\sqrt{\nu _e}`$ varies from 0.67 for $`l=20m`$ down to 0.47 for $`l=150m`$. This is a modest loss. We also stress that the naive $`\sqrt{N}`$ scaling is not adequate for an appearance experiment where we are looking for few events and hence we conclude that genuine $`\nu _\mu `$ rate is more important than a slightly better $`\nu _e/\nu _\mu `$ ratio.

Remains the issue of the proton energy. Naively, one would expect that the higher the proton energy, the higher is the background. We find however that the intrinsic electron neutrino background does not strongly depend on the proton energy.

This was verified directly for various proton energies and baselines. Results of the calculation are shown in Table 8, all normalized to $`10^{19}`$ pots. We observe that (1) for the shortest baselines, the ratio $`\nu _e/\nu _\mu `$ is increasing dramatically with proton energy. Accordingly, (2) the ratio $`\nu _\mu /\sqrt{\nu _e}`$ decreases. However, we must rescale this ratio to take into account the approximate scaling of the number of events with the proton energy. Since we expect

$$\frac{\nu _\mu }{\sqrt{\nu _e}}\frac{E_p}{\sqrt{E_p}}\sqrt{E_p}$$
(30)
we consider the rescaled ratio $`\nu _\mu /\sqrt{\nu _eE_p}`$ in the last columns of the Table 8. These are also plotted in Figure 17 as a function of the baseline $`L`$. Numerically, we find that the rescaled ratios $`\nu _\mu /\sqrt{\nu _eE_p}`$ are almost the identical at the optimal baselines of each proton energy, so not much is too be gained by varying the proton energy.

Summarizing, we find that the electron neutrino background in the relevant region is not dependent on the proton energy and only determined by the decay tunnel length, but that its optimisation is very limited. For the maximum neutrino muon flux, it is at the level of the 1% for any of the considered setups.

## 8 Summary and Conclusions

In this document, we have performed a two-dimensional scan, varying the beam energy and baseline parameters to optimize the conditions for the investigation of $`\theta _{13}`$ driven neutrino oscillations in the whole Superkamiokande allowed $`\mathrm{\Delta }m^2`$ range. We find that:

* The optimal baselines for $`\theta _{13}`$ searches are in the range 100-700 km for proton energies varying from 2.2 to 400 GeV.
* The needed beam intensity scales approximately with the inverse of the beam energy.
* In terms of proton economics, the optimum beam energy is around 20 GeV, but lower beam energies are appealing for the shortest baselines.

Realistic focusing system for low and medium baselines have also been studied. In this case:

* Focusing efficiencies of 30-50% can be achieved in the energy range of interest.

The whole procedure is detector and accelerator independent. Nevertheless, to give a first estimate of the needed beam intensities, the ICARUS detector has been assumed as a reference.

We can draw the following observations:

* a 2.2 GeV or 4.4 GeV high-intensity proton machine (i.e. à la CERN-SPL) is well matched to a baseline in the range 100-300 km. It is not matched to a baseline of 730 km (i.e. CERN-LNGS).
* a 20 GeV machine is best matched to a baseline of 730 km (i.e. CERN-LNGS). However, an integrated intensity in the range of $`10^{23}`$ pots are required, which is about two orders of magnitude higher than the intensity deliverable by the current CERN-PS in a reasonable amount of time.
* a 400 GeV energy is reasonably matched to a baseline of 730 km (i.e. CERN-LNGS). For the 400 GeV, the required intensity is in the range of $`10^{22}`$ pots, which is about one order of magnitude higher than the intensity deliverable by the current CERN-SPS in a reasonable amount of time.

Finally, we stress that the present study is essentially a theoretical one. All the “real” work has still to be accomplished in order for one of these options to become reality.