text:
# IFIC/02-38 UWThPh-2002-25 Constraining Majorana neutrino electromagnetic properties from the LMA-MSW solution of the solar neutrino problem

## 1 Introduction

Present solar Cleveland:nv , Abdurashitov:2002nt , GaGNO , Fukuda:2002pe , SNONC , SNODN  and atmospheric neutrino data atm02 , Fukuda:1998mi  provide the first robust evidence for neutrino flavour conversion and, consequently, the first solid indication for new physics. Neutrino oscillations constitute the most popular explanation for the data (for recent analyses see Refs. Maltoni:2002ni , solfits , 3flavors ) and are a natural outcome of gauge theories of neutrino mass Schechter:1980gr , rev . Non-zero neutrino masses manifest themselves through non-standard neutrino electromagnetic properties. If the lepton sector in the Standard Model (SM) were minimally extended in analogy with the quark sector, neutrinos get Dirac masses ($`m_\nu `$) and their magnetic moments (MMs) are tiny Fujikawa:yx ,

$$\mu _\nu 3\times 10^{19}\mu _B\left(\frac{m_\nu }{1\mathrm{eV}}\right),$$
(1)
where $`\mu _B`$ is the Bohr magneton. Laboratory experiments give 90% C.L. bounds on the neutrino MMs of $`1.8\times 10^{10}\mu _B`$ Derbin:wy , Derbin:ua  and $`1.3\times 10^{10}\mu _B`$ munu  for the electron neutrino, $`6.8\times 10^{10}\mu _B`$ for the muon neutrino LSND  and $`3.9\times 10^7\mu _B`$ for the tau neutrino DONUT  (see also Ref. Hagiwara:pw ). On the other hand, astrophysics and cosmology provide limits of the order of $`10^{12}`$ to $`10^{11}`$ Bohr magnetons Raffelt:gv . Improved sensitivity for the electron neutrino MM from reactor neutrino experiments is expected, while a tritium $`\overline{\nu }_e`$ source experiment H3  aims to reach $`3\times 10^{12}\mu _B`$.

It has for a long time been noticed, on quite general “naturality” grounds, that Majorana neutrinos constitute the typical outcome of gauge theories Schechter:1980gr . On the other hand, precisely such neutrinos also emerge in specific classes of unified theories, in particular, in those employing the seesaw mechanism rev . If neutrinos are indeed Majorana particles the structure of their electromagnetic properties differs crucially from that of Dirac neutrinos Schechter:1981hw , being characterized by a $`3\times 3`$ complex antisymmetric matrix $`\lambda `$, the so-called Majorana transition moment (TM) matrix. It contains MMs as well as electric dipole moments of the neutrinos. The existence of any electromagnetic neutrino moment well above the expectation in Eq. (1) would signal that some special mechanism—which goes beyond the SM—is at work. Therefore, neutrino electromagnetic properties are interesting probes of new physics.

As noted in Ref. Schechter:1981hw  the existence of Majorana transition moments leads to the phenomenon of spin-flavour precession, when they move through a magnetic field, as might happen in the Sun. This was subsequently shown to affect the propagation of solar neutrinos in an important way, due to the effects of matter Akhmedov:uk , and this has recently been shown to provide excellent fits to current data Miranda .

In the present paper we will take a different attitude. Motivated by the now robust status of the LMA-MSW solution of the solar neutrino problem indicated by recent data Maltoni:2002ni , solfits , and also by the prospects of its possible confirmation by the KamLAND experiment kamland , we take up this solution as the basis of our investigation. It has been shown in Ref. Miranda  that in this case the effect of a solar magnetic field can be safely neglected. Therefore, our present approach is justified and complementary to that in Ref. Miranda : there we were concerned with spin-flavour-conversions of solar neutrinos, sensitive to the product of $`\lambda `$ times a solar magnetic field, while here we are concerned with the effect of TMs in the neutrino detection process. Sizable TMs would contribute to the elastic neutrino–electron scattering in the Super-Kamiokande (SK) experiment Fukuda:2002pe , and hence, data from this experiment can be used to constrain electromagnetic neutrino properties Beacom:1999wx , joshipura .

Here we will tacitly employ the results of the latest analysis in Ref. Maltoni:2002ni  as well as the fit strategies used there, to perform a fit to the global solar neutrino data in terms of LMA-MSW oscillation parameters and neutrino TMs, in order to derive bounds on the latter. It is important to clarify for which combinations of the three complex moments the bounds are valid. We will show that all three moments can be bounded, however, not all at the same level. In addition, we will show that data from reactor neutrino experiments Derbin:wy , Derbin:ua , munu  which use elastic neutrino–electron scattering for neutrino detection are to some extent complementary to the solar neutrino data, so that by combining both data sets one can improve the bounds. We also explore the potential of the upcoming Borexino experiment borexino , Alimonti:2002xc  in probing the TMs, showing that it is substantially more sensitive to the TMs than present experiments.

The paper is organized as follows. In Section 2 we give our theoretical framework and we discuss in detail which elements of the TM matrix are constrained. In Section 3 we present the statistical method we apply to calculate bounds on the TMs, and in Section 4 we describe our analysis of solar and reactor neutrino data. The resulting bounds are discussed in Section 5 and the sensitivity of the upcoming Borexino experiment is considered in Section 6. A brief discussion concludes the paper in Section 7.

## 2 The theoretical framework

In this paper, we adopt the following basic working assumptions:

* Only three light neutrinos exist, doublets under $`SU(2)U(1)`$ (no light sterile neutrinos); this is well-motivated by recent global fits of all available neutrino oscillation data, including also those that follow from short baseline searches Maltoni:2002xd .
* The solution of the solar neutrino puzzle is given by the LMA-MSW solution, as indicated by recent global fits of solar data solfits , 3flavors .
* Neutrinos have Majorana nature as expected from theory rev , Schechter:1981hw .
* Neutrinos are endowed with TMs arising from some unspecified electroweak gauge theory.

Although the latter is probably the least motivated of our assumptions, there are some models of the electroweak interaction where the possibility of enhanced TMs has been discussed rev .

In experiments where the neutrino detection reaction is elastic neutrino–electron scattering, like in SK, Borexino and some reactor experiments Derbin:wy , Derbin:ua , munu , the electromagnetic cross section is Bardin:wr

$$\frac{d\sigma _{\mathrm{em}}}{dT}=\frac{\alpha ^2\pi }{m_e^2\mu _B^2}\left(\frac{1}{T}\frac{1}{E_\nu }\right)\mu _{\mathrm{eff}}^2,$$
(2)
where the effective MM square is given by GS

$$\mu _{\mathrm{eff}}^2=a_{}^{}\lambda ^{}\lambda a_{}+a_+^{}\lambda \lambda ^{}a_+.$$
(3)
The electromagnetic cross section adds to the weak cross section and allows to extract information on the TM matrix $`\lambda `$. In this cross section, $`T`$ denotes the kinetic energy of the recoil electron, $`E_\nu `$ is the energy of the incoming neutrino and the 3-vectors $`a_{}`$ and $`a_+`$ denote the neutrino amplitudes for negative and positive helicities, respectively, at the detector. The square of the effective MM given in Eq. (3) is independent of the basis chosen for the neutrino state GS . In what follows we consider both the flavour basis and the mass eigenstate basis. We will use the convention that $`a_{}`$ and $`\lambda =(\lambda _{\alpha \beta })`$ denote the quantities in the flavour basis, whereas in the mass basis we will use $`\stackrel{~}{a}_{}`$ and $`\stackrel{~}{\lambda }=(\lambda _{jk})`$. The two bases are connected via the neutrino mixing matrix $`U`$:

$$\stackrel{~}{a}_{}=U^{}a_{},\stackrel{~}{a}_+=U^Ta_+,\stackrel{~}{\lambda }=U^T\lambda U.$$
(4)
One usually decomposes the transition moment matrix as

$$\lambda =\mu id,$$
(5)
where $`\mu `$ and $`d`$ are hermitian matrices. On general grounds these matrices are, in addition, antisymmetric and imaginary for Majorana neutrinos Schechter:1981hw  (see also Ref. GS  and references cited therein). It is, furthermore, useful to define vectors $`𝚲=(\mathrm{\Lambda }_\alpha )`$ and $`\stackrel{~}{𝚲}=(\mathrm{\Lambda }_j)`$ in the flavour and mass basis, respectively, by

$$\lambda _{\alpha \beta }=\epsilon _{\alpha \beta \gamma }\mathrm{\Lambda }_\gamma \mathrm{and}\lambda _{jk}=\epsilon _{jkl}\mathrm{\Lambda }_l.$$
(6)
Thus, in the flavour basis we have $`\lambda _{e\mu }=\mathrm{\Lambda }_\tau `$, $`\lambda _{\mu \tau }=\mathrm{\Lambda }_e`$ and $`\lambda _{\tau e}=\mathrm{\Lambda }_\mu `$. Note also that

$$|𝚲|^2=\frac{1}{2}\mathrm{Tr}\left(\lambda ^{}\lambda \right)|𝚲|=|\stackrel{~}{𝚲}|.$$
(7)
This means if we are able to find a bound on $`|𝚲|`$, we have not only constrained the TMs in the flavour basis but also in the mass basis. In the numerical section of the paper, this will be exactly our strategy.

Let us now discuss the detailed form that the effective MM square $`\mu _{\mathrm{eff}}^2`$ in Eq. (3) takes when we assume the LMA-MSW solution of the solar neutrino problem, denoted by $`\mu _{\mathrm{LMA}}^2`$. Numerical integrations of the neutrino evolution equations performed in Ref. Miranda  showed that for an effective MM of $`10^{11}\mu _B`$ a solar magnetic field of the order of 80 KGauss has practically no effect on the LMA-MSW solution. The results of Ref. Miranda  imply that it is safe to neglect possible solar magnetic field effects on the evolution of the neutrino state in the sun in case of LMA-MSW. This follows essentially from the fact that the terms in the evolution Hamiltonian related to the LMA-MSW mass-squared difference and the matter potential are much bigger than the electromagnetic terms. Hence, helicity is conserved in solar neutrino propagation, so that we can set $`a_+=0`$, thereby eliminating the second term in Eq. (3). Another good approximation used in fitting solar data in a three-neutrino scenario 3flavors  is to set $`\mathrm{\Delta }m^2atm\mathrm{\Delta }m^2sol\mathrm{\Delta }m^2`$, i.e., oscillations with $`\mathrm{\Delta }m^2atm`$ are averaged out. Using the parameterization (6) we obtain from Eq. (3)

$$\mu _{\mathrm{LMA}}^2=|𝚲|^2P_{e3}^{3\nu }|\mathrm{\Lambda }_3|^2\underset{j,k=1}{\overset{2}{}}(\stackrel{~}{a}_{}^j)(\stackrel{~}{a}_{}^k)^{}\mathrm{\Lambda }_j^{}\mathrm{\Lambda }_k.$$
(8)
The brackets $`\mathrm{}`$ in the last term in Eq. (8) denote the average over the production point, earth-sun distance and zenith angle, and we have made use of the relation

$$P_{ej}^{3\nu }|\stackrel{~}{a}_{}^j|^2(j=1,2,3)$$
(9)
for the probability that the neutrino produced in the core of the sun as a $`\nu _e`$ arrives at the detector as a mass eigenstate $`\nu _j`$. From Eq. (8) one learns that besides the absolute values of the TM$`|\mathrm{\Lambda }_1|,|\mathrm{\Lambda }_2|,|\mathrm{\Lambda }_3|`$ only one complex phase $`\delta =\mathrm{arg}(\mathrm{\Lambda }_1^{}\mathrm{\Lambda }_2)`$ enters into $`\mu _{\mathrm{LMA}}^2`$.

Eq. (8) can be further simplified by making use of the large mass-squared difference $`\mathrm{\Delta }m^25\times 10^5`$ eV2 implied by the LMA-MSW solution. First, because of vacuum oscillations on the way from the sun to the earth the neutrino state arriving at the earth can be considered as an incoherent mixture of mass eigenstates GS , dighe . Note also that for the LMA-MSW solution, earth matter effects are very small earthmatter  so that we neglect them in our expression for $`\mu _{\mathrm{LMA}}^2`$.111Note, however, that the neglection of earth matter effects is not a very good approximation when performing a combined fit of TMs and oscillation parameters. Therefore, in this case we do take into account earth matter effects in the minimization of the $`\chi ^2`$ with respect to the solar oscillation parameters. Then we find that only the diagonal elements of the matrix $`\left(\stackrel{~}{\lambda }^{}\stackrel{~}{\lambda }\right)`$ contribute, and as a result only the probabilities (9) appear in determining the detected neutrino signal. Setting $`s_{13}\mathrm{sin}\theta _{13}`$, $`c_{13}\mathrm{cos}\theta _{13}`$ and $`U_{e3}=s_{13}`$ one has, to a very good approximation 3flavors ,

$$P_{e3}^{3\nu }=s_{13}^2\mathrm{and}P_{ej}^{3\nu }=c_{13}^2P_{ej}^{2\nu }(j=1,2)\mathrm{with}P_{e1}^{2\nu }+P_{e2}^{2\nu }=1.$$
(10)
The probabilities $`P_{ej}^{2\nu }`$ $`(j=1,2)`$ are the effective 2-neutrino probabilities for the solar neutrino problem, where all the averages mentioned above have been taken into account. Using now the probabilities in Eq. (10) and the approximations mentioned above, Eq. (8) becomes

$$\mu _{\mathrm{LMA}}^2=|\mathrm{\Lambda }_1|^2+|\mathrm{\Lambda }_2|^2+c_{13}^2|\mathrm{\Lambda }_3|^2c_{13}^2\underset{j=1,2}{}P_{ej}^{2\nu }|\mathrm{\Lambda }_j|^2.$$
(11)
It is well known from reactor neutrino data CHOOZ  that the mixing angle $`\theta _{13}`$ is rather small. Fits to the CHOOZ data and atmospheric neutrino data 3flavors  show that at 3$`\sigma `$ one has $`s_{13}^2\stackrel{<}{}\mathrm{\hspace{0.25em}0.05}`$. This allows us to replace $`c_{13}^2`$ in Eq. (11) by 1. Then we arrive at the final formula

$$\mu _{\mathrm{LMA}}^2=|𝚲|^2|\mathrm{\Lambda }_2|^2+P_{e1}^{2\nu }\left(|\mathrm{\Lambda }_2|^2|\mathrm{\Lambda }_1|^2\right),$$
(12)
which will be used in the numerical section evaluating solar neutrino data. The probability $`P_{e1}^{2\nu }`$ is a function of the ratio $`\mathrm{\Delta }m^2/E_\nu `$ and the solar mixing angle $`\theta \theta _{12}`$. In the following we will drop the super-script $`2\nu `$ and $`P_{e1}`$ or $`P_{ee}`$ refers always to a 2-neutrino probability. Eq. (12) naturally makes no distinction between MMs and electric dipole moments. Since we aim at constraining $`|𝚲|`$, all elements of the TM matrix will be bounded at the same time.

Now we move to the case of reactor neutrinos. There we have a pure $`\overline{\nu }_e`$ source. Therefore, we have $`a_{}=0`$ and $`a_+=(1,0,0)^T`$ since in these experiments the baseline is much too short for any oscillations to develop. The resulting $`\mu _{\mathrm{eff}}^2`$ relevant in reactor experiments is given as

$$\mu _\mathrm{R}^2=|\mathrm{\Lambda }_\mu |^2+|\mathrm{\Lambda }_\tau |^2.$$
(13)
From this relation it is clear that reactor data on its own cannot constrain all TMs contained in $`\lambda `$, since $`\mathrm{\Lambda }_e`$ does not enter in Eq. (13). In order to combine reactor and solar data it is useful to rewrite $`\mu _\mathrm{R}^2`$ in terms of the mass basis quantities. With Eq. (6) we readily derive

$$\left(\stackrel{~}{\lambda }\stackrel{~}{\lambda }^{}\right)_{jk}=|𝚲|^2\delta _{jk}\mathrm{\Lambda }_j\mathrm{\Lambda }_k^{}.$$
(14)
Then with Eq. (4) we obtain

$$\mu _\mathrm{R}^2=|𝚲|^2c^2|\mathrm{\Lambda }_1|^2s^2|\mathrm{\Lambda }_2|^22sc|\mathrm{\Lambda }_1||\mathrm{\Lambda }_2|\mathrm{cos}\delta ,$$
(15)
where $`c=\mathrm{cos}\theta `$ and $`s=\mathrm{sin}\theta `$, $`\theta `$ being the solar mixing angle. Further, we notice that the relative phase $`\delta =\mathrm{arg}(\mathrm{\Lambda }_1^{}\mathrm{\Lambda }_2)`$ between $`\mathrm{\Lambda }_1`$ and $`\mathrm{\Lambda }_2`$ appears in addition to $`|𝚲|`$, $`|\mathrm{\Lambda }_1|`$ and $`|\mathrm{\Lambda }_2|`$.

## 3 Statistical method and qualitative discussion

In the following we will use data from solar and reactor neutrino experiments to constrain neutrino TMs. The $`\chi ^2`$-function obtained from the data depends on the solar oscillation parameters $`\mathrm{tan}^2\theta `$ and $`\mathrm{\Delta }m^2`$ as well as on the elements of the TM matrix $`\lambda `$. Regarding the dependence on the oscillation parameters we will take two different attitudes. One is to assume that $`\mathrm{tan}^2\theta `$ and $`\mathrm{\Delta }m^2`$ will be determined with good accuracy at the KamLAND experiment, and hence, we will consider the $`\chi ^2`$ at fixed points in the $`\mathrm{tan}^2\theta \mathrm{\Delta }m^2`$ plane (method I). In the second approach we will derive a bound on the TMs by minimizing the $`\chi ^2`$ with respect to $`\mathrm{tan}^2\theta `$ and $`\mathrm{\Delta }m^2`$ (method II). This second procedure takes into account the present uncertainty of our knowledge of the oscillation parameters.

Let us describe in detail how we calculate a bound on the TMs. Our aim is to constrain $`|𝚲|`$, therefore it is convenient to consider the $`\chi ^2`$ as a function of the independent parameters $`|𝚲|,|\mathrm{\Lambda }_1|,|\mathrm{\Lambda }_2|`$ and $`\delta `$. As discussed in the previous section, $`\delta `$ appears only if reactor data are included. The $`\chi ^2`$-functions which we are using for the individual data sets (solar rates, SK recoil electron spectrum, reactor data) will be described in detail in the following sections. When performing the fit to the data, we find that in general the minimum of the $`\chi ^2`$ occurs close or outside the physical boundary of the parameters $`|𝚲|,|\mathrm{\Lambda }_1|,|\mathrm{\Lambda }_2|`$. To take this into account we apply Bayesian techniques to calculate an upper bound on $`|𝚲|`$. We minimize the $`\chi ^2`$ with respect to $`|\mathrm{\Lambda }_1|`$, $`|\mathrm{\Lambda }_2|`$ and $`\delta `$ for each value of $`|𝚲|`$, taking into account the allowed region $`0|\mathrm{\Lambda }_1|^2+|\mathrm{\Lambda }_2|^2|𝚲|^2`$:

$$\chi ^2(|𝚲|)=\mathrm{Min}\left[\chi ^2(|𝚲|,|\mathrm{\Lambda }_1|,|\mathrm{\Lambda }_2|,\delta )\right].$$
(16)
In method I we do this for fixed values of the oscillation parameters, scanning over the LMA-MSW region, whereas in method II we minimize also with respect to $`\mathrm{tan}^2\theta `$ and $`\mathrm{\Delta }m^2`$ in Eq. (16). Then the $`\chi ^2`$ is transformed into a likelihood function via

$$\mathrm{exp}\left(\frac{1}{2}\chi ^2\right).$$
(17)
Now we use Bayes’ theorem and a flat prior distribution $`p(|𝚲|)`$ in the physically allowed region, $`p(|𝚲|)=\mathrm{\Theta }(|𝚲|)`$, to obtain a probability distribution for $`|𝚲|`$:

$$f(|𝚲|)d|𝚲|=\frac{(|𝚲|)\mathrm{\Theta }(|𝚲|)d|𝚲|}{_0^{\mathrm{}}(|𝚲|^{})d|𝚲|^{}}.$$
(18)
An upper bound $`b_\alpha `$ on $`|𝚲|`$ at a C.L. $`\alpha `$ is given by the equation

$$_0^{b_\alpha }f(|𝚲|)d|𝚲|=\alpha .$$
(19)
Let us consider in more detail the minimization with respect to $`|\mathrm{\Lambda }_1|`$ and $`|\mathrm{\Lambda }_2|`$ in Eq. (16) in the case of solar data without reactor. As we will show later there is no evidence for a non-zero $`\mu _{\mathrm{LMA}}^2`$ in the data. Therefore the minimum of the $`\chi ^2`$ for a given $`|𝚲|`$ will occur if $`\mu _{\mathrm{LMA}}^2`$ is minimal. Departing from Eq. (12) it is easy to show that

$$\mathrm{Min}[\mu _{\mathrm{LMA}}^2]=\{\begin{array}{cc}|𝚲|^2P_{e1}\hfill & \mathrm{for}P_{e1}0.5,\hfill \\ |𝚲|^2(1P_{e1})\hfill & \mathrm{for}P_{e1}>0.5.\hfill \end{array}$$
(20)
The minimum occurs at

$$|𝚲|^2=\{\begin{array}{cc}|\mathrm{\Lambda }_2|^2\text{for}\hfill & P_{e1}<0.5,\hfill \\ |\mathrm{\Lambda }_1|^2+|\mathrm{\Lambda }_2|^2\text{for}\hfill & P_{e1}=0.5,\hfill \\ |\mathrm{\Lambda }_1|^2\text{for}\hfill & P_{e1}>0.5.\hfill \end{array}$$
(21)
From Eq. (20) follows that the bound on $`|𝚲|`$ is strongest if $`P_{e1}=0.5`$ because in this situation $`\mu _{\mathrm{LMA}}^2`$ is maximal. In Fig. 1 we show contours of constant $`P_{e1}`$ in the $`\mathrm{tan}^2\theta \mathrm{\Delta }m^2/E_\nu `$ plane. For definiteness, the probabilities in the figure are obtained by performing the averaging over the production distribution inside the sun for the 7Be flux most relevant for Borexino. However, the probabilities for the 8B flux relevant for SK are very similar.

The shaded region in Fig. 1 is the region relevant for Borexino, assuming a mass-squared difference in the range $`10^5`$ eV$`{}_{}{}^{2}<\mathrm{\Delta }m^2<10^4`$ eV2, whereas for SK the region below the dash-dotted line is most important, due to the higher energy of the 8B neutrinos. We can read off from the figure that in most part of the SK region $`P_{e1}`$ is very small. According to Eq. (21) this means that the sensitivity of SK for $`|𝚲|`$ is limited because of the small $`P_{e1}`$. In contrast, we expect a much better sensitivity of Borexino, because, in a large part of the relevant parameter space in this case, $`P_{e1}`$ is close to the optimal value of $`0.5`$.

The arguments presented here are valid as long as the $`\chi ^2`$ depends only on the parameter combination given in Eq. (12). The data from reactor experiments depend on a different combination of TM parameters (see Eq. (15)) implying that the combined analysis of SK and reactor data will potentially lead to more stringent limits than the SK data alone.

## 4 Analysis of solar and reactor neutrino data

In this section we describe briefly our analysis of solar neutrino data and the data which we are using from reactor neutrino experiments. Neglecting the electromagnetic contribution to the small elastic scattering signal in SNO, one notices that the neutrino TMs contribute directly to the observed signal only in SK. However, the uncertainty in the 8B flux leads to correlations between SK and all other experiments. Therefore, even for fixed oscillation parameters (method I), also other experiments give some information on TMs and it is important to include them in the analysis joshipura .

We divide the solar neutrino data into the total rates observed in all experiments and the SK recoil electron spectrum (with free overall normalization). Before combining rates with SK spectra, we derive bounds on $`|𝚲|`$ for each of these cases separately. Subsequently we consider the reactor data by themselves and in combination with the global sample of solar data.

### 4.1 Solar rates

In this analysis we use the event rates measured in the Chlorine experiment Homestake Cleveland:nv  (Cl), the Gallium experiments Sage Abdurashitov:2002nt , Gallex and GNO GaGNO  (GaGNO), the total event rate of SK Fukuda:2002pe  based on the 1496 days data sample and the recent result of the Sudbury neutrino observatory SNONC  for the charged current (SNO CC) and neutral current (SNO NC) solar neutrino rates. For the analysis we use the $`\chi ^2`$-function

$$\chi _{\mathrm{rates}}^2=\underset{j,k}{}(R_jD_j)(V_{jk}^{\mathrm{rate}})^1(R_kD_k).$$
(22)
Here the indices $`j,k`$ run over the 6 solar neutrino rates (Cl, Sage, GaGNO, SK, SNO CC, SNO NC), $`D_j`$ are the experimental rates, $`R_j`$ are the theoretical predictions and $`V_{jk}^{\mathrm{rate}}`$ is the covariance matrix.

The data we are using are given in Table 1. The rates of SK and SNO are the ratio of the measured flux and the 8B flux predicted by the standard solar model (SSM) BP00 BP00  and no disappearance of solar $`\nu _e`$ ($`\mathrm{\Phi }_\mathrm{B}^0=5.05\times 10^6`$ cm-2s-1). Note that the SNO CC and NC fluxes as given in Ref. SNONC  have been obtained from a fit to the SNO spectral data assuming no distortion of the neutrino energy spectrum, i.e., a survival probability constant in energy, which in general is not realized in the case of neutrino oscillations. However, in the case of the LMA-MSW solution it is a very good approximation to use the fluxes given in Ref. SNONC .

The theoretical predictions $`R_j`$ which appear in Eq. (22) are calculated as described in Ref. Maltoni:2002ni . However, the rate $`R_{\mathrm{SK}}`$ of the SK experiment includes an extra contribution from the electromagnetic scattering. This rate is given by

$$R_{\mathrm{SK}}=\frac{_0^{\mathrm{}}𝑑E_\nu \varphi _\mathrm{B}(E_\nu )\left\{P_{ee}(E_\nu )\sigma _e^W(E_\nu )+[1P_{ee}(E_\nu )]\sigma _\mu ^W(E_\nu )+\mu _{\mathrm{LMA}}^2(E_\nu )\sigma _{\mathrm{em}}^{}(E_\nu )\right\}}{_0^{\mathrm{}}𝑑E_\nu \varphi _\mathrm{B}(E_\nu )\sigma _e^W(E_\nu )},$$
(23)
where $`\varphi _\mathrm{B}(E_\nu )`$ is the shape of the 8B flux normalized to 1, which we take from Ref. bahcallpage , and

$$\sigma (E_\nu )=_0^{T_{\mathrm{max}}(E_\nu )}𝑑T^{}\frac{d\sigma (T^{},E_\nu )}{dT^{}}_{T_1}^{T_2}𝑑TR(T^{},T)\text{with}\sigma =\sigma _e^W,\sigma _\mu ^W,\sigma _{\mathrm{em}}^{}.$$
(24)
We use the SM weak cross sections $`\sigma _{e,\mu }^W`$ including radiative corrections Bahcall:1995mm , $`\mu _{\mathrm{LMA}}^2`$ is given in Eq. (12), and

$$\frac{d\sigma _{\mathrm{em}}^{}(T^{},E_\nu )}{dT^{}}=\frac{\alpha ^2\pi }{m_e^2\mu _B^2}\left(\frac{1}{T^{}}\frac{1}{E_\nu }\right).$$
(25)
Here $`E_\nu `$ is the neutrino energy, while $`T^{}`$ and $`T`$ denote the true and the measured recoil electron kinetic energies, respectively. The integration in $`T`$ is from the SK threshold 5 MeV (total electron energy) up to 20 MeV. The maximum kinetic energy of the recoil electrons is given by $`T_{\mathrm{max}}(E_\nu )=2E_\nu ^2/(2E_\nu +m_e)`$. The resolution function $`R(T^{},T)`$ is taken from Ref. Fukuda:1998ua .

The covariance matrix $`V_{jk}^{\mathrm{rate}}`$ in the $`\chi ^2`$ of Eq. (22) takes into account the experimental errors and theoretical uncertainties from detection cross sections $`V_{jk}^{\mathrm{CS}}`$ and SSM predictions $`V_{jk}^{\mathrm{flux}}`$ added in quadrature fogli ,

$$V_{jk}^{\mathrm{rate}}=V_{jk}^{\mathrm{exp}}+V_{jk}^{\mathrm{CS}}+V_{jk}^{\mathrm{flux}}.$$
(26)
The experimental errors for different experiments are uncorrelated. For $`j,k=`$ Cl, Sage, GaGNO, SK we have $`V_{jk}^{\mathrm{exp}}=\delta _{jk}\sigma _j^2`$, where the $`\sigma _j^2`$ are the experimental uncertainties, calculated by adding in quadrature the statistical and systematic uncertainties given in Table 1. However, the statistical and systematic errors for the SNO CC and NC rates are strongly correlated, with $`\rho ^{\mathrm{stat}}=0.518`$ (the same as the correlation for the CC and NC day-night asymmetry, given in Ref. SNODN ) and $`\rho ^{\mathrm{sys}}=0.508`$ (which we derive from Table II of Ref. SNONC ). The second and third terms in Eq. (26) take into account the uncertainties in the detection cross sections and in the SSM predictions of the neutrino fluxes, respectively. For details see Refs. Maltoni:2002ni , fogli  and references therein.

### 4.2 The Super-Kamiokande recoil electron spectrum

In this section we consider the shape of the SK recoil electron spectrum. The electromagnetic contribution to the neutrino–electron scattering cross section leads to a substantially different spectrum of the scattered electrons than expected from the SM weak interaction. Therefore, the SK spectral data provide a useful tool to constrain neutrino TMBeacom:1999wx .

We perform a fit to the latest 1496 days SK data presented in Ref. Fukuda:2002pe . There the event rates are given in 8 energy bins ranging from 5 to 20 MeV. The energy bins 2 to 7 are further divided into 7 bins of zenith angle, which makes up a total of 44 data points $`D_i`$. We define the $`\chi ^2`$-function

$$\chi _{\mathrm{spect}}^2=\underset{i,j=1}{\overset{44}{}}(\alpha R_iD_i)(V_{ij}^{\mathrm{spect}})^1(\alpha R_jD_j).$$
(27)
The theoretical prediction $`R_i(\mathrm{tan}^2\theta ,\mathrm{\Delta }m^2,\mu _{\mathrm{LMA}}^2)`$ for the $`i`$-th bin is calculated as in Eq. (23) with the integration interval $`[T_1,T_2]`$ in Eq. (24) chosen according to the energy interval of the given bin. The covariance matrix $`V_{ij}^{\mathrm{spect}}`$ in Eq. (27) contains statistical and systematic experimental errors. These are obtained from Ref. Fukuda:2002pe  taking into account the correct correlations of systematic errors. The $`\chi ^2`$ of Eq. (27) is minimized with respect to the normalization factor $`\alpha `$ in order to isolate only the shape of the spectrum.

### 4.3 Reactor data

Data from neutrino electron scattering at nuclear reactor experiments Derbin:wy , Derbin:ua , munu  can be used to constrain the combination of TMs given in Eq. (15). Here we use data from the Rovno nuclear power plant Derbin:wy  and from the Bugey reactor munu . To include this information in our analysis we make the following ansatz for the $`\chi ^2`$-function:

$$\chi _{\mathrm{reactor}}^2(\mu _\mathrm{R})=\underset{i}{}\left(\frac{N_{\mathrm{weak}}^i+N_{\mathrm{em}}^i(\mu _\mathrm{R})N_{\mathrm{obs}}^i}{\sigma ^i}\right)^2.$$
(28)
The sum is over the two experiments Rovno and Bugey, $`N_{\mathrm{obs}}^i`$ is the observed number of events with the one standard deviation error $`\sigma ^i`$, $`N_{\mathrm{weak}}^i`$ is the number of events expected in the case of no neutrino TMs (only the standard weak interaction) and $`N_{\mathrm{em}}^i(\mu _\mathrm{R})`$ is the number of events due to the electromagnetic scattering of neutrinos with an effective MM $`\mu _\mathrm{R}`$. We can write the latter as

$$N_{\mathrm{em}}^i(\mu _\mathrm{R})=C^i\left(\frac{\mu _\mathrm{R}}{10^{10}\mu _B}\right)^2.$$
(29)
For the Rovno reactor we obtain from Table I of Ref. Derbin:wy  the values $`N_{\mathrm{obs}}^{\mathrm{Rovno}}=41`$, $`\sigma ^{\mathrm{Rovno}}=20`$ and $`N_{\mathrm{weak}}^{\mathrm{Rovno}}=32`$. Furthermore, in this table the number of events due to the electromagnetic scattering of neutrinos with an effective MM $`\mu _\mathrm{R}=2\times 10^{10}\mu _B`$ is given: $`N_{\mathrm{em}}^{\mathrm{Rovno}}(2\times 10^{10}\mu _B)=54`$. Using this last information the constant in Eq. (29) can be determined as $`C^{\mathrm{Rovno}}=54/4`$. Results of the MUNU experiment performed at the Bugey reactor have been presented recently at the Neutrino 2002 conference munu . These lead to the bound $`\mu _\mathrm{R}<1.3\times 10^{10}\mu _B`$ at 90 % C.L., if an energy cut of $`E_e>1`$ MeV is used. From Ref. munu  we get (in events per day) $`N_{\mathrm{obs}}^{\mathrm{Bugey}}=0.37`$, $`\sigma ^{\mathrm{Bugey}}=0.25`$ and $`N_{\mathrm{weak}}^{\mathrm{Bugey}}=0.45`$. The constant in Eq. (29) can be calculated by making use of the 90% C.L. bound cited above: we solve the equation $`\chi _{\mathrm{Bugey}}^2(\mu _\mathrm{R}=1.3\times 10^{10}\mu _B)=2.7`$. This leads to $`C^{\mathrm{Bugey}}=0.196`$. In Fig. 2 we show the $`\mathrm{\Delta }\chi ^2`$ of the two experiments and also the sum. Note that with our ansatz Eqs. (28) and (29) we can reproduce to good accuracy the 1 $`\sigma `$ upper bound $`\mu _\mathrm{R}<1.5\times 10^{10}\mu _B`$ from the Rovno data Derbin:wy .

## 5 Bounds on $`|𝚲|`$ from solar and reactor data

First we discuss the bounds on $`|𝚲|`$ from solar data alone. Fixing the oscillation parameters at the current best fit point $`\mathrm{tan}^2\theta =0.44`$, $`\mathrm{\Delta }m^2=6.6\times 10^5`$ eV2 Maltoni:2002ni  we obtain with the Bayesian methods described in Section 3 the 90% C.L. bound

$$|𝚲|<4.0\times 10^{10}\mu _B\text{(best fit point, solar data)}.$$
(30)
However, such a bound substantially depends on the values of the neutrino oscillation parameters. In Fig. 3 we show contours of the 90% C.L. bound on $`|𝚲|`$ in the $`\mathrm{tan}^2\theta \mathrm{\Delta }m^2`$ plane. We find that the bound gets weaker for small values of $`\mathrm{\Delta }m^2`$, whereas in the upper left part of the LMA-MSW region a bound of the order $`2\times 10^{10}\mu _B`$ is obtained.

By combining solar and reactor data we obtain considerably stronger bounds. At the best fit point we get at 90% C.L.

$$|𝚲|<2.0\times 10^{10}\mu _B\text{(best fit point, solar + reactor data)}.$$
(31)
In Fig. 4 we show the contours of the bound in the $`\mathrm{tan}^2\theta \mathrm{\Delta }m^2`$ plane for the combination of solar and reactor data. By comparing Fig. 4 with Fig. 3 we find that reactor data drastically improve the bound for low $`\mathrm{\Delta }m^2`$ values.

In order to better understand these results we give in Tab. 2 bounds on $`|𝚲|`$ for five characteristic points in the $`\mathrm{tan}^2\theta \mathrm{\Delta }m^2`$ plane, including the current best fit point, which is labeled with a star. The other points are marked in Figs. 3 and 4 by crosses and are labeled A, B, C, D, and the corresponding values of the oscillation parameters are given in Tab. 2. Besides the bound on $`|𝚲|`$ we give also bounds on the individual $`|\mathrm{\Lambda }_i|(i=1,2,3)`$. These bounds are obtained by setting the other two $`\mathrm{\Lambda }_i`$ to zero. We show the bounds for the different data samples: solar rates, the shape of the SK spectrum with free overall normalization (SK spect), solar rates combined with the SK spectrum (solar global), reactor data only, solar data combined with reactor data (sol + react).

Notice that, in contrast with the bound on $`|𝚲|`$ which, thanks to Eq. (7), is valid in any basis, the bounds on individual $`|\mathrm{\Lambda }_i|`$ do not have a basis-independent meaning. The bounds on $`|\mathrm{\Lambda }_i|`$ given in Tab. 2 refer to the mass-eigenstate basis, also used in Ref. Schechter:1981hw .

By comparing the first and the second column of Tab. 2 we notice that both solar rates and SK spectrum are important: although the bound from the SK spectrum is slightly stronger, the one which follows from solar rates is in general of the same order. The bound we obtain by combining rates and SK spectrum is displayed in the third column. Concerning the reactor data, clearly these bounds do not depend on $`\mathrm{\Delta }m^2`$, and the bound on $`|\mathrm{\Lambda }_3|`$ is even completely independent of the neutrino oscillation parameters (see Eq. (15)). As mentioned above, reactor data on their own cannot be used to constrain $`|𝚲|`$.

In the upper parts of the LMA-MSW region, the solar data alone give already a strong bound on $`|𝚲|`$. From the bounds given in Tab. 2 at point A we see that for such high values of $`\mathrm{\Delta }m^2`$ all $`|\mathrm{\Lambda }_i|`$ are strongly constrained by solar data. This is due to the fact that there the probability $`P_{e1}`$ relevant in SK is close to the optimal value of 0.5 (see Fig. 1). In this parameter region the combination with reactor data does not improve the bound significantly. This is also true for the region of small $`\mathrm{tan}^2\theta `$ around point C.

The behavior of the bounds for low $`\mathrm{\Delta }m^2`$ values can be understood by looking at point B. In this $`\mathrm{\Delta }m^2`$ region the probability $`P_{e1}`$ is very small for the neutrino energies relevant for SK, as seen in Fig. 1. Therefore, solar data give only a very weak constraint on $`|\mathrm{\Lambda }_2|`$ as can be seen from Tab. 2, and the bound on $`|𝚲|`$ is dominated by this weak bound on $`|\mathrm{\Lambda }_2|`$ (see Eq. (21)) leading to the rather poor bound of $`11.6\times 10^{10}\mu _B`$ from solar data. In this case the combination with reactor data improves the bound drastically to $`2.3\times 10^{10}\mu _B`$.

Up to now we have calculated bounds on neutrino TMs for fixed values of the oscillation parameters $`\mathrm{tan}^2\theta `$ and $`\mathrm{\Delta }m^2`$ (method I described in Section 3). These results will be especially useful after KamLAND will have determined the oscillation parameters with good accuracy. In the following we change our strategy and minimize the $`\chi ^2`$ for each value of $`|𝚲|`$ with respect to $`\mathrm{tan}^2\theta `$ and $`\mathrm{\Delta }m^2`$ (method II). This will lead to a bound on $`|𝚲|`$ taking into account the current knowledge concerning the oscillation parameters. To this end we make use of the global solar neutrino data (including earth matter effects) as described in Ref. Maltoni:2002ni  in order to obtain the correct $`\chi ^2`$ behavior as a function of $`\mathrm{tan}^2\theta `$ and $`\mathrm{\Delta }m^2`$. Only in the expression for $`\mu _{\mathrm{LMA}}^2`$ we neglect earth matter effects as described in Section 2, since they are known to be small in the LMA-MSW region earthmatter . Performing this analysis we obtain the following bounds at 90% C.L.:

$$|𝚲|<\{\begin{array}{cc}6.3\times 10^{10}\mu _B\hfill & \text{(unconstrained, solar data)}\hfill \\ 2.0\times 10^{10}\mu _B\hfill & \text{(unconstrained, solar + reactor data).}\hfill \end{array}$$
(32)
Together with Figs. 3 and 4 the bounds given in Eq. (32) constitute an important result of this paper. They show that, assuming the solar LMA-MSW solution, current solar neutrino data can be used to constrain all elements of the Majorana neutrino TM matrix. A combination with data from reactor experiments significantly strengthens the bound on $`|𝚲|`$.

## 6 Simulation of the Borexino experiment

Here we investigate the sensitivity of the Borexino experiment borexino , Alimonti:2002xc  to neutrino TMs. This experiment is mainly sensitive to the solar 7Be neutrino flux, which will be measured by elastic neutrino–electron scattering. Therefore, Borexino is similar to SK, the main difference is the mono-energetic line of the 7Be neutrinos, with an energy of 0.862 MeV, which is roughly one order of magnitude smaller than the energies of the 8B neutrino flux relevant in SK.

### 6.1 The Borexino $`\chi ^2`$-function

To estimate the sensitivity of Borexino we consider the following $`\chi ^2`$-function:

$$\chi _{\mathrm{borexino}}^2=\underset{i,j=1}{\overset{N_{\mathrm{bins}}}{}}(N_i^{\mathrm{th}}D_i)(V_{ij}^{\mathrm{borex}})^1(N_j^{\mathrm{th}}D_j).$$
(33)
Here $`N_i^{\mathrm{th}}`$ is the theoretical prediction for the number of events in the electron recoil energy bin $`i`$, and $`D_i`$ is the (hypothetical) observed number of events. Borexino will observe recoil electrons with kinetic energy in the range 0.25 to 0.8 MeV. As a realistic example we adopt $`N_{\mathrm{bins}}=8`$ bins in electron recoil energy and an energy resolution of 0.058 MeV borexino , Alimonti:2002xc , Berezhiani:2001rt . We have checked that our results are essentially independent of the exact value of the energy resolution.

The theoretical prediction for the number of events in the bin $`i`$ after $`N_{\mathrm{years}}`$ years of Borexino data taking is given by

$$N_i^{\mathrm{th}}=N_i^{\mathrm{weak}}+N_i^{\mathrm{em}}(\mu _{\mathrm{LMA}})+N_i^{\mathrm{bg}},N_i^x=365N_{\mathrm{years}}n_i^x$$
(34)
with $`x=`$ th, weak, em, bg, and $`n_i^{\mathrm{weak},\mathrm{em},\mathrm{bg}}`$ is the number of events per day from weak scattering, electromagnetic scattering and the background, respectively. If not stated otherwise we consider a running time of $`N_{\mathrm{years}}=3`$ years. In order to estimate the sensitivity of Borexino for neutrino TMs we assume that the data are generated by neutrinos without TMs:

$$D_i=N_i^{\mathrm{weak}}+N_i^{\mathrm{bg}}.$$
(35)
Hence, we obtain

$$\chi _{\mathrm{borexino}}^2=\underset{i,j=1}{\overset{N_{\mathrm{bins}}}{}}N_i^{\mathrm{em}}(V_{ij}^{\mathrm{borex}})^1N_j^{\mathrm{em}}.$$
(36)
The minimum of this $`\chi ^2`$ occurs for $`\mu _{\mathrm{LMA}}=0`$ and is always zero. With the $`\chi ^2`$ of Eq. (36) we calculate a bound on $`|𝚲|`$ at a given C.L. This bound corresponds to the maximum allowed value of $`|𝚲|`$ which cannot be distinguished from $`|𝚲|=0`$, and is therefore a measure for the obtainable sensitivity to $`|𝚲|`$ at Borexino.

The numbers of events due to the weak and electromagnetic scattering are calculated by

$`n_i^{\mathrm{weak}}`$ $`=`$ $`𝒩{\displaystyle \underset{f}{}}\mathrm{\Phi }_f^0{\displaystyle _0^{\mathrm{}}}𝑑E_\nu \varphi _f\left\{P_{ee}^f[\sigma _e^W]_i+(1P_{ee}^f)[\sigma _\mu ^W]_i\right\},`$
$`n_i^{\mathrm{em}}`$ $`=`$ $`𝒩{\displaystyle \underset{f}{}}\mathrm{\Phi }_f^0{\displaystyle _0^{\mathrm{}}}𝑑E_\nu \varphi _f\mu _{\mathrm{LMA}}^2[\sigma _{\mathrm{em}}^{}]_i,`$ (37)
respectively. The effective cross sections in the bin $`i`$ are calculated like in Eqs. (24) and (25), by choosing the lower and upper bounds of the recoil energy interval $`T_1`$ and $`T_2`$ according to the given bin $`i`$. We sum over the four most important neutrino fluxes ($`f`$ = Be, pep, N, O) in Borexino (see Ref. Alimonti:2002xc ). For the mono-energetic lines we have $`\varphi _{\mathrm{Be}}(E_\nu )=\delta (E_\nu 0.862\mathrm{MeV})`$ and $`\varphi _{\mathrm{pep}}(E_\nu )=\delta (E_\nu 1.442\mathrm{MeV})`$, and we use the normalized spectra $`\varphi _\mathrm{N}(E_\nu )`$ and $`\varphi _\mathrm{O}(E_\nu )`$ given in Ref. bahcallpage . The absolute values of the neutrino fluxes $`\mathrm{\Phi }_f^0`$ from Ref. BP00  are given by $`\mathrm{\Phi }_{\mathrm{Be}}^0=0.429,\mathrm{\Phi }_{\mathrm{pep}}^0=0.014,\mathrm{\Phi }_\mathrm{N}^0=0.055,\mathrm{\Phi }_\mathrm{O}^0=0.048`$ in units of $`10^{10}`$ cm-2s-1. Note that from the two 7Be lines the line at 0.38 MeV, which constitutes 10% of the total 7Be flux, does not contribute to the signal in Borexino. Therefore the value for $`\mathrm{\Phi }_{\mathrm{Be}}^0`$ given above has been obtained by multiplying the value given in Ref. BP00  with 0.9.

In Table 3 we show the expected number of events per day in Borexino for the kinetic energy of the recoil electrons in the range 0.25 to 0.8 MeV, resulting from different solar neutrino fluxes, for the standard solar model without neutrino oscillations. We do not consider other neutrino fluxes, because they contribute to the signal with less than 0.2 events per day. We fix the normalization constant $`𝒩`$ in Eq. (37) in such a way that the number of events per day from the 7Be line for $`P_{ee}=1`$ and $`\mu _{\mathrm{LMA}}=0`$ is 43.3. Then we can reproduce the number of events arising from the other fluxes as given in Table 3 to good accuracy.

The number of background events $`n_i^{\mathrm{bg}}`$ is obtained in the following way. We read off the shape of the background as a function of the electron energy from Fig. 15 of Ref. Alimonti:2002xc  and normalize this function to 19 background events per day in the energy range 0.25 to 0.8 MeV borexino . Then we can calculate $`n_i^{\mathrm{bg}}`$ by integrating over the corresponding energy interval for each bin.

We use the following covariance matrix in Eq. (36):

$$V_{ij}^{\mathrm{borex}}=\delta _{ij}N_i^{\mathrm{th}}+\underset{f_1,f_2}{}N_{if_1}N_{jf_2}\underset{\beta =1}{\overset{12}{}}\alpha _{f_1\beta }\alpha _{f_2\beta }\left(\mathrm{\Delta }\mathrm{ln}X_\beta \right)^2+N_i^{\mathrm{bg}}N_j^{\mathrm{bg}}(\mathrm{\Delta }\mathrm{ln}N^{\mathrm{bg}})^2.$$
(38)
The first term is the statistical uncertainty, taken as the square-root of the predicted number of events. The second term describes the uncertainty in the solar neutrino fluxes. Here $`N_{if}`$ are the contributions of the individual solar neutrino fluxes to the total event numbers in each bin: $`N_i^{\mathrm{th}}=_fN_{if}+N_i^{\mathrm{bg}}`$, and the sum is over the four fluxes relevant in Borexino (for definition and values of the quantities $`\alpha _{f\beta }`$ and $`X_\beta `$ see Ref. fogli ). The last term takes into account a systematic uncertainty in the number of background events. We assume full correlation between the bins and that the relative uncertainty is the same for all bins: $`\mathrm{\Delta }N_i^{\mathrm{bg}}/N_i^{\mathrm{bg}}=\mathrm{\Delta }\mathrm{ln}N^{\mathrm{bg}}`$. We adopt a standard value for our calculations of $`\mathrm{\Delta }\mathrm{ln}N^{\mathrm{bg}}=10\%`$.

### 6.2 Sensitivity on $`|𝚲|`$ of Borexino

Using the statistical method described in Sections 3 and 6.1 we obtain for the current best fit point $`\mathrm{\Delta }m^2=6.6\times 10^5`$ eV2, $`\mathrm{tan}^2\theta =0.44`$ the upper bound (sensitivity)

$$|𝚲|0.28\times 10^{10}\mu _B\text{at}90\%\text{C.L.}$$
(39)
after three years of Borexino data taking. This bound is about one order of magnitude stronger than the bound from existing data. We have checked that a combined analysis of Borexino with existing data (solar as well as reactor data) does not improve the bound of Eq. (39). The bound depends on the actual value of the oscillation parameters; in Fig. 5 we show contours of the 90% C.L. bound in the $`\mathrm{tan}^2\theta \mathrm{\Delta }m^2`$ plane. The strongest attainable limit is roughly $`0.24\times 10^{10}\mu _B`$. In agreement with the discussion we gave in connection with Eq. (21) we find that the strongest bound is obtained when $`P_{e1}=0.5`$, as shown by the dotted line in Fig. 5.

We have performed several tests concerning our assumptions for Borexino. Fig. 6 (left panel) shows the dependence of the bound on the number of bins in the recoil electron energy. We can see from this figure that a measurement of the total rate alone ($`N_{\mathrm{bins}}=1`$) gives a rather weak bound. However a few bins are already enough for a considerable improvement of the bound on $`|𝚲|`$. It is interesting to see that the error on the background is only important for few bins, and becomes completely unimportant for more than 3 or 4 bins, because the main information comes from the spectral shape.

In the right panel in Fig. 6 we show the bound as a function of the Borexino running time. Here we adopt again our nominal value of 8 energy bins. After one year a sensitivity of $`0.4\times 10^{10}\mu _B`$ can be obtained with Borexino, which is already about one order of magnitude stronger than the current bound. After three years of data-taking only a minor improvement of the bound is possible, mainly due to the uncertainty in the solar neutrino fluxes.222Note that we do not include any experimental systematic uncertainties in our Borexino analysis.

## 7 Conclusions

In this paper we have presented stringent bounds on electromagnetic Majorana transition moments (TMs). Such TMs—if present—would contribute to the elastic neutrino–electron scattering signal in solar neutrino experiments like Super-Kamiokande or the upcoming Borexino experiment, as well as in reactor neutrino experiments. Motivated by the robust status of the LMA-MSW solution of the solar neutrino problem indicated by recent data, and also by the prospects of its possible confirmation by the KamLAND experiment, we have taken this solution as basis for our investigation. In this context we have clarified which effective magnetic moments are probed in solar and reactor experiments—see Eqs. (12) and (13), respectively. Using most recent global solar neutrino data we have performed a fit in terms of the oscillation parameters and the elements of the complex TM matrix $`\lambda `$ of three active Majorana neutrinos. Taking into account the antisymmetry of the TM matrix by the parameterization $`\lambda _{jk}=\epsilon _{jkl}\mathrm{\Lambda }_l`$, we have shown that solar neutrino data allow to constrain $`|𝚲|`$, i.e., all elements of the TM matrix are bounded at the same time. We want to stress that the bounds on $`|𝚲|`$ hold in any basis, because $`|𝚲|`$ is an intrinsic neutrino property; in terms of TMs it is given by

$$|𝚲|^2=|\lambda _{e\mu }|^2+|\lambda _{\mu \tau }|^2+|\lambda _{\tau e}|^2=|\lambda _{12}|^2+|\lambda _{23}|^2+|\lambda _{31}|^2$$
(40)
in the flavour and mass-eigenstate bases, respectively. On the other hand, the bounds on the individual $`|\mathrm{\Lambda }_i|`$ ($`i=1,2,3`$) in Tab. 2 refer only to the mass-eigenstate basis.

A fit to the global solar neutrino data leads to the bound $`|𝚲|<6.3\times 10^{10}\mu _B`$ at 90% C.L. We have also considered the role of reactor neutrino data on neutrino TMs, shown to be complementary to solar neutrino data. A combined fit of reactor and solar data leads to significantly improved bounds: at 90% C.L. we get $`|𝚲|<2.0\times 10^{10}\mu _B`$. In the very near future the KamLAND experiment will crucially test the LMA-MSW solution. If KamLAND confirms the LMA-MSW solution then a precise determination of the parameters $`\mathrm{tan}^2\theta `$ and $`\mathrm{\Delta }m^2`$ might be possible. This motivated us to scan the $`\mathrm{tan}^2\theta \mathrm{\Delta }m^2`$ plane and to calculate the corresponding bound on $`|𝚲|`$ in each point. These results are shown in Figs. 3 and 4 for solar data and solar + reactor data, respectively.

We have also investigated the potential of the upcoming neutrino–electron scattering solar neutrino experiment Borexino to constrain neutrino TMs. Performing a detailed simulation of the experiment we find that it will improve the bound on $`|𝚲|`$ by about one order of magnitude with respect to present bounds.

Last, but not least, let us mention that alternative solutions to the solar neutrino problem based on non-standard neutrino matter interactions Guzzo:2001mi , which may arise in models of neutrino masses NSImodels , can also be constrained in an analogous way Berezhiani:2001rt . If, as seems likely, the LMA-MSW solution is finally borne out, the improved determination of the oscillation parameters expected, say at KamLAND, will provide an extremely useful stepping stone for probing other features of neutrino physics beyond the Standard Model.

## Acknowledgments

This work was supported by Spanish grant PB98-0693, by the European Commission RTN network HPRN-CT-2000-00148 and by the European Science Foundation network grant N. 86. T.S. has been supported by the DOC fellowship of the Austrian Academy of Science and, in the early stages of this work, by a fellowship of the European Commission Research Training Site contract HPMT-2000-00124 of the host group. M.M. was supported by contract HPMF-CT-2000-01008 and M.A.T. was supported by the M.E.C.D. fellowship AP2000-1953.


format: text

text:
# Homogeneous isotropic turbulence in dilute polymers: scale by scale budget E. De Angelis Dip. Mecc. Aeron., Università di Roma La Sapienza, via Eudossiana 18, 00184, Roma, Italy., C.M. Casciola∗, R. Benzi Dipartimento di Fisica e INFM, Università di Roma Tor Vergata, Via della Ricerca scientifica 1, 00133 Roma, Italy, & R. Piva∗

## Abstract

The turbulent energy cascade in dilute polymers solution is addressed here by considering a direct numerical simulation of homogeneous isotropic turbulence of a FENE-P fluid in a triply periodic box. On the basis of the DNS data, a scale by scale analysis is provided by using the proper extension to visco-elastic fluids of the Karman-Howarth equation for the velocity. For the microstructure, an equation, analogous to the Yaglom equation for scalars, is proposed for the free-energy density associated to the elastic behavior of the material. Two mechanisms of energy removal from the scale of the forcing are identified, namely the classical non-linear transfer term of the standard Navier-Stokes equations and the coupling between macroscopic velocity and microstructure. The latter, on average, drains kinetic energy to feed the dynamics of the microstructure. The cross-over scale between the two corresponding energy fluxes is identified, with the flux associated with the microstructure dominating at small separations to become sub-leading above the cross-over scale, which is the equivalent of the elastic limit scale defined by De Gennes-Tabor on the basis of phenomenological assumptions.

## I Introduction

Turbulence in dilute polymers solutions is still an open issue, despite the growing interest on the subject. Most of the efforts have been addressed to the comprehension of practical aspects, such as drag reduction, related to the modification of turbulence by long chain polymers in wall bounded flows, . However, in most cases, the attempts have proven inconclusive, suggesting that, even from the point of view of applications, a more fundamental approach is required. Under this respect, basic studies should focus on simpler flow conditions where turbulence-polymers interaction occurs.

A possible choice is to exploit homogeneity and isotropy to simplify the statistical treatment of the data and to make use of the available, relatively well established, rheological models for dilute polymers, see  for a comprehensive review.

In this context, a recent experimental investigation on decaying grid turbulence of dilute polymers and surfactants  shows, for the polymers, a substantial alteration of the small scales, though still in presence of a substantial degree of anisotropy, ascribed to the grid shear layers.

In this and other related experiments, see e.g.  and , the measurements were mostly aimed at the velocity, but also at the pressure field see e.g.  for a recent example. Actually, as one can imagine, the experimental analysis of turbulent flows of polymeric liquids is particularly difficult, and hardly one can proceed beyond the global characterization of the flow and the analysis of certain aspects of the macroscopic field, see ,  as additional examples in the more complex configuration of a channel flow. This makes the actual mechanism of interaction between polymers and turbulence particularly obscure, and leads to consider numerical simulations as a viable tool to address the problem. For triply-periodic boundary conditions, in particular, fast and accurate spectral methods can easily be developed to achieve the direct numerical simulation of a strictly homogeneous isotropic flow, once a reasonably accurate and sufficiently simple rheological model has been selected. Among these, the so-called FENE-P model  is possibly the best compromise between accurate representation of polymers dynamics and minimal computational complexity, see . In fact, recent results show that this model is able to reproduce the drag reducing behavior of dilute polymers in wall bounded flows, see , , .

As often in rheology, one can easily get stuck in a discussion about the proper model for a specific application. We take here the model for granted and try to grasp the mechanisms by which it affects the turbulence. To do this we take freedom to chose parameter values which may not correspond too well to experimental conditions, but which may help revealing the underlying physics.

In summary, our model for the solution considers an ensemble of elastic dumbbells attached to each material point of the continuum. The dumbbells are stretched by the flow due to friction and react through a non-linear elastic spring. The probability density function of the end-to-end vector of the dumbbells is described through a second order tensor, the conformation tensor, assumed to represent the covariance matrix of the end-to-end vectors. It obeys a transport equation which accounts for material stretching and for the elastic reaction. The force exerted on the continuum by the dumbells contributes an extra-stress term in the macroscopic momentum balance. On average, the work done by the extra-stress against the macroscopic velocity field amounts to a draining of energy from the macroscopic field. The corresponding power feeds the dynamics of the microstructure, is partially accumulated as free-energy of the ensemble and is eventually dissipated, see .

In terms of turbulence dynamics, several basic questions arise. First of all, one should understand if the nature of this kind of viscoelastic turbulence is substantially similar, or, on the contrary, essentially different from standard Newtonian turbulence. According to the common view, a turbulent flow is described as a superposition of different scales of motion. In homogeneous isotropic turbulence of ordinary fluids the interaction of the different scales originates a flux of energy from the large towards the small scale, with energy injected by an external forcing at large scales and dissipated by viscosity at small scales, the so-called direct, or forward, energy cascade (see e.g. the book by Frisch ). The presence of the microstructure is expected to alter this process, since it corresponds to an alternative path the energy flux may take to achieve the eventual dissipation. If this is the case, one may wonder if the alteration is localized at certain spatial scales or if the entire spectrum is affected by the energy draining of the polymers. This poses the additional issue of determining a possible scale which delimits the range of scales where the presence of the polymers is effective from those where they are substantially passive, see . This scale may be close to the classical Kolmogorov dissipation scale, but it may even be substantially larger, and it is not even clear if the effect of polymers may be reduced to an additional dissipation, , or if the role of elastic energy is crucial, , see also .

All these issues can be addressed systematically by analyzing the data of a DNS of homogeneous isotropic turbulence once proper diagnostic tools are available. To this purpose we present and discuss a scale by scale budget based on the extension to viscoelastic fluids of the classical Karman-Howarth  equation of ordinary turbulence, see e.g. ,  and  for similar extensions in the context of shear flows. By this equation we can address the effect of polymers on the macroscopic kinetic field. To consider the fluctuations in the polymers we would need to extend the classical Yaglom equation for a scalar, see , to a second order active tensor. In a more straightforward way instead we consider, in the proper thermodynamic framework, the Yaglom equation for the free-energy of the polymers.

## II The evolution equations for dilute polymer solutions

As discussed in the introduction, the momentum balance for a dilute solution of long chain polymers in otherwise Newtonian, incompressible solvent is described by a slightly modified form of the Navier-Stokes equation, namely

$`{\displaystyle \frac{Du_i}{Dt}}`$ $`=`$ $`{\displaystyle \frac{p}{x_i}}+\nu {\displaystyle \frac{^2u_i}{x_jx_j}}+{\displaystyle \frac{T_{ij}}{x_j}}+f_i,`$ (1)
where $`D/Dt=/t+u_k/x_k`$ is the usual substantial derivative, $`f_i`$ the external forcing, $`p`$ is the pressure normalized by the constant density of the solution, $`\nu `$ is the kinematic viscosity of the solvent and $`T_{ij}`$ is the additional contribution to the stress tensor - extra-stress - due to the polymers chains. The macroscopic velocity field $`u_i`$ is solenoidal and the constitutive relation for the extra-stress is

$`T_{ij}`$ $`=`$ $`\nu _p/\tau (f(R_{kk};\rho _m,\rho _0)R_{ij}/\rho _0^2\delta _{ij}),`$ (2)
where $`\nu _p`$ is a constant typically of the order of a small but finite fraction of $`\nu `$. The conformation tensor, $`R_{ij}`$, is a second order tensor taken to characterize the statistical behavior of the population of polymers chains attached to a given point, and $`f(R_{kk};\rho _m,\rho _0)`$ is a dimensionless nonlinear spring coefficient specified as $`f=(\rho _m^2\rho _0^2)/(\rho _m^2R_{kk})`$, with $`\rho _0`$ the equilibrium length of the chains and $`\rho _m`$ their maximum allowed length, as standard in the context of the finite extensibility nonlinear elastic model (FENE, see ). The dynamics of the population of chains is described by an evolution equation for the conformation tensor,

$`{\displaystyle \frac{DR_{ij}}{Dt}}`$ $`=`$ $`{\displaystyle \frac{u_i}{x_k}}R_{kj}+R_{ik}{\displaystyle \frac{u_k}{x_j}}{\displaystyle \frac{1}{\tau }}(fR_{ij}\rho _0^2\delta _{ij}),`$ (3)
where $`\tau `$ is the principal relaxation time of the chains. The system (1, 2, 3) completed with the continuity equation,

$`{\displaystyle \frac{u_i}{x_i}}`$ $`=`$ $`0,`$ (4)
for the macroscopic velocity forms the FENE-P model for dilute polymers solutions, see e.g. .

Phenomenologically, stretching of the chains by the flow, described by the terms involving the velocity gradient in (3), is counteracted by an elastic reaction, the term in brackets, while the polymers are transported and re-oriented by the flow. The reaction of the polymers alters the macroscopic force balance, as seen in eq. (1), where the divergence of the extra-stress enters as an additional source or sink of momentum.

The extra-stress is able to make work against the velocity as follows from the balance of kinetic energy

$`{\displaystyle \frac{D(u^2/2)}{Dt}}={\displaystyle \frac{(p\delta _{ij}+\mathrm{\Sigma }_{ij})u_i}{x_j}}+{\displaystyle \frac{(T_{ij}u_i)}{x_j}}ϵ_N+S,`$ (5)
where $`\mathrm{\Sigma }_{ij}=2\nu e_{ij}`$ is the standard viscous component of the stress tensor and $`ϵ_N=2\nu e_{ij}e_{ij}`$ is the Newtonian component of the dissipation, with $`e_{ij}`$ the macroscopic velocity deformation tensor. With respect to the kinetic energy balance of the standard Navier-Stokes equation, two additional terms arise here, to describe the interaction of the macroscopic field with the microstructure. The one in divergence form contributes to the spatial redistribution of kinetic energy. The other,

$`S=T_{ij}{\displaystyle \frac{u_i}{x_j}},`$ (6)
called hereafter the stress-power, represent energy per unit time that the microstructure can drain, or, in principle, even release to the macroscopic kinetic field. The presence of this term is physically due to the fact that the polymers chains can both dissipate, via friction with the solvent, or store energy, due to their elasticity. The correct thermodynamical setting refers to the free-energy $`a`$ of the polymers ensemble. Given the expression, see e.g. ,

$`a={\displaystyle \frac{1}{2}}{\displaystyle \frac{\nu _p}{\tau }}[(\rho _{max}^2/\rho _0^21)\mathrm{log}\left({\displaystyle \frac{\rho _{max}^2R_{kk}}{\rho _{max}^2\rho _0^2}}\right)`$ (7)
(8)
$`+1/3\mathrm{log}(detR/\rho _0^2)],`$ (9)
it follows directly from the evolution equation (3) that the energy balance for the polymers is

$`{\displaystyle \frac{Da}{Dt}}`$ $`=`$ $`Sϵ_P,`$ (10)
where

$`ϵ_P={\displaystyle \frac{1}{2}}{\displaystyle \frac{\nu _p}{\tau ^2}}f\left[(fR)_{kk}+(fR)_{kk}^16\right]`$ (11)
a positive definite quantity, is the rate of energy dissipation associated with the polymers. Adding eq. (5) and (10) gives the balance of the total free-energy of the system, $`E=u^2/2+a`$, as

$`{\displaystyle \frac{DE}{Dt}}`$ $`=`$ $`f_iu_iϵ_T+{\displaystyle \frac{J_k}{x_k}}`$ (12)
where $`J_i`$ is the spatial flux of total free-energy, whose expression follows immediately by inspection of the kinetic energy balance. Note that $`S`$ drops from the global balance since it corresponds to a conservative exchange of energy between the two constituents, namely the macroscopic kinetic energy and the polymers free-energy.

## III The Karman-Howarth equation

In a statistically stationary turbulent system, the external forcing $`f_i`$ provides the energy to maintain the turbulent fluctuations which otherwise would be damped by the dissipation. For the present case both the macroscopic velocity and the microstructure exhibits turbulent fluctuations. Characteristic feature of the system is that, while the external forcing acts only on the macroscopic field, both the macroscopic flow and the microstructure contribute to the dissipation of the injected mechanical power. This is possible since the extra-stress captures part of the available mechanical energy to feed the fluctuations of the microstructure and the associated dissipation. In order to try to disentangle these complex processes, the simplest setting is provided by homogeneous and isotropic turbulence once suitable balance equation are available for the fluctuations of both kinds.

Purpose of the present section is to present the equation appropriate to address the scale by scale budget of the macroscopic velocity fluctuations. This is the extention to dilute polymer solutions of the well known Karman-Howarth equation of homogeneous isotropic turbulence for Newtonian fluids. As discussed in the previous section, where the FENE-P model has been briefly recalled, the solenoidal velocity field, $`u_i`$, obeys a modified form of the Navier-Stokes equation, where the divergence of the extra-stress $`T`$, denoted here by

$`g_i`$ $`=`$ $`{\displaystyle \frac{T_{ij}}{x_j}},`$ (13)
appears as an additional source of momentum.

Starting with the modified Navier-Stokes equation, the equation for the correlation tensor $`C_{i,j}=<u_iu_j^{}>`$, where unprimed and primed variables are evaluated at $`x_i`$ and at $`x_i^{}`$, respectively, with $`x_i^{}`$ the displacement of $`x_i`$ by the separation vector $`r_i`$, easily follows from the standard procedure, see e.g. , . The trace of this equation reads

$`{\displaystyle \frac{C_{i,i}}{t}}+{\displaystyle \frac{}{r_k}}<u_iu_i^{}u_k^{}u_iu_i^{}u_k>=`$ (14)
$`2\nu {\displaystyle \frac{^2C_{i,i}}{r_kr_k}}+<g_i^{}u_i+g_iu_i^{}>+<f_i^{}u_i+f_iu_i^{}>,`$ (15)
where $`f_i`$ is here specified as a random stirring force acting on the large scales. To obtain the above equation we have used standard properties of homogeneous incompressible turbulence, which allow to drop the pressure term and to recast spatial derivatives with derivatives with respect to the separation, see e.g. .

Equation (14) can be re-expressed in terms of velocity increments, $`\delta V_i=u_i(x_r^{})u_i(x_r)`$, $`\delta V^2=\delta V_i\delta V_i`$, once recalled that $`C_{i,i}=<u_iu_i>1/2<\delta V^2>`$, and that, by standard manipulations, we have

$`<\delta V^2\delta V_k>=2<u_iu_i^{}u_k^{}u_iu_i^{}u_k>+`$ (17)
$`<u_iu_iu_k^{}><u_i^{}u_i^{}u_k>.`$ (18)
Hence, for a solenoidal velocity field, the equation in terms of velocity increments follows as

$`{\displaystyle \frac{}{r_k}}<\delta V^2\delta V_k+2T_{ki}^{}\delta V_i>=4<f_iu_i>+`$ (19)
$`2<\delta V_i\delta f_i>+\mathrm{\hspace{0.17em}2}\nu {\displaystyle \frac{^2}{r_kr_k}}<\delta V^2>.`$ (20)
In deriving eq. (19) a statistically stationary field has been assumed and the term $`<f_i^{}u_i+f_iu_i^{}>`$ has been rearranged in the form $`2<f_iu_i><\delta V_i\delta f_i>`$ to extract explicitly the power injected by the external forcing. Finally, the integration by parts

$`<g_i^{}u_i+g_iu_i^{}>={\displaystyle \frac{}{r_k}}(T_{ki}^{}u_iT_{ki}u_i^{})={\displaystyle \frac{}{r_k}}(T_{ki}^{}\delta u_i),`$ (22)
where $`T_{ki}^{}=T_{ki}^{}+T_{ki}`$, leads to eq. (19).

Equation (19) can be integrated in $`r`$-space over a ball $`B_r`$ of radius $`r`$ (see e.g.  for a similar approach) to yield

$`{\displaystyle \frac{1}{4\pi r^2}}{\displaystyle _{B_r}}<\delta V^2\delta V_k+2T_{ki}^{}\delta V_i>n_k𝑑S_r=4/3\overline{ϵ}_Tr`$ (24)
$`+{\displaystyle \frac{\nu }{2\pi r^2}}{\displaystyle \frac{d}{dr}}{\displaystyle _{B_r}}<\delta V^2>dS_r`$ (25)
$`+{\displaystyle \frac{1}{2\pi r^2}}{\displaystyle _{B_r}}<\delta V_i\delta f_i>dV_r,`$ (26)
where the global energy balance (12), which for a stationary homogeneous field reads

$`<f_iu_i>`$ $`=`$ $`\overline{ϵ}_T,`$ (28)
has been accounted for.

The Karman-Howarth equation (24) highlights the main difference between Newtonian and viscoelastic turbulence. For either case, in the range of scales where the direct effect of viscosity and the correlation between velocity and forcing increments are both negligible, the energy flux occurring through scale $`r`$ equals the total dissipation. For Newtonian fluids, all the dissipation is provided by the viscosity and the flux through scale $`r`$ only occurs due to the classical nonlinearity associated to the advection. For viscoelastic fluids, an additional dissipative process takes place in the polymers, and, consistently, the flux through scale $`r`$ also presents a new component in addition to that of the standard Navier-Stokes dynamics. In fact, the additional flux term describes the amount of energy intercepted by the microstructure. Though the classic component is expected to be associated to a forward energy cascade, i.e. the corresponding flux enters the ball of radius $`r`$, concerning the viscoelastic contribution we can, at the moment, only conjecture that energy intercepted at scale $`r`$ is passed forward to the microstructure at smaller scale, consistently with the idea of polymer dissipation as a small scale process. However, since, to our knowledge, the process has never been investigated before, we have to examine and discuss the numerical simulation before drawing any conclusion.

## IV The Yaglom equation for the free-energy

The Karman-Howarth equation for dilute polymer solution may be used to analyze the scale by scale budget for the velocity field. Its Fourier transform actually corresponds to the budget of turbulent kinetic energy in the different modes of the spectrum. To address the complete dynamics of dilute polymers solution we need a similar equation concerning the microstructure. Under this respect, in section 1 we have shown that a free-energy $`a`$ can be associated to the ensemble of polymers attached to a given point in space and that, locally, the energy balance for the microstructure is provided by equation (10).

The free-energy $`a`$ is a measure of the average elongation of the polymers as given by the trace of the conformation tensor $`R_{ij}`$. In fact one can derive an equation for the scale by scale budget of the trace though we prefer to discuss the Yaglom equation for $`a`$, whose meaning in terms of energetics of the system is more definite. To this purpose, let us denote by

$`\psi =Sϵ_p`$ (29)
the instantaneous excess power, i.e. the amount of power transfered to the polymers in excess to the dissipation they originate. Strictly speaking this is an improper naming, since on the average the excess power is zero. Nonetheless, it measures the amount of energy stored per unit time in the polymers. From (10) we can derive the equation for the correlation $`<aa^{}>`$,

$`{\displaystyle \frac{<aa^{}>}{t}}+{\displaystyle \frac{}{r_k}}<u_k^{}a^{}au_kaa^{}>`$ $`=`$ (30)
$`<a^{}\psi +a\psi ^{}>,`$ (31)
which, considering that

$`<a\psi ^{}+a^{}\psi >=2<a\psi ><\delta a\delta \psi >,`$

and noting that, for a stationary state ($`<aDa/Dt>=0`$),

$`<a\psi >=0,`$ (32)
can be recast in the form

$`{\displaystyle \frac{}{r_k}}<\delta u_k\delta a^2>=2<\delta a\delta \psi >.`$ (33)
Integration over the sphere of radius $`r`$ then yields

$`{\displaystyle \frac{1}{4\pi r^2}}{\displaystyle _{B_r}}<\delta a^2\delta V_k>n_k𝑑S_r`$ $`=`$ (34)
$`{\displaystyle \frac{2}{4\pi r^2}}{\displaystyle _{B_r}}<\delta a\delta \psi >𝑑V_r,`$ (35)
which, by using isotropy reduces to

$`<\delta a^2\delta V_{}>`$ $`=`$ $`2{\displaystyle _0^r}<\delta a\delta \psi >r^2𝑑r,`$ (36)
with $`\delta V_{}=\delta V_ir_i/r`$ the longitudinal velocity increment.

Equation (36) is the appropriate form of the Yaglom equation for the polymers. It states the equivalence between the flux of free energy through the different spatial scales, as described by $`<\delta a^2\delta V_{}>`$ and the integral of the correlation between free energy fluctuations, $`\delta a`$, and fluctuation of the excess of power, $`\delta \psi `$. Classically for scalars the Yaglom equation provides a term linear in the separation associated to the average scalar dissipation. In the present case, such term is zero as consequence of eq. (32).

## V Global features of the flow and DNS

The system of equations (1,3,4), corresponding to $`9`$ scalar equation, $`3`$ velocity components plus $`6`$ components of the conformation tensor, has been integrated using a Fourier $`\times `$ Fourier $`\times `$ Fourier spectral method for spatial discretization and a four-stages/third order Runge-Kutta method for time advancement with all the non-linear terms fully de-aliased by the three-halves rule. An exactly solenoidal velocity is achieved by the standard projection method. The side of the cubic computational box is $`2\pi `$ and $`64`$ modes are used in each direction. A fully developed field of Newtonian turbulence has been used as initial condition for the viscoelastic calculation and the random forcing is applied to the first shell of wave-vectors, with constant amplitude and uniformly distributed phases. The random phase extracted at each step is used subject to the condition of positive energy injection, otherwise it is changed by $`\pi `$.

The nominal Reynolds number of the simulation can be constructed from the intensity of the forcing, $`\widehat{f}_0`$, the length of the box, $`l_x`$, and the zero shear kinematic viscosity, $`\nu _T=\nu +\nu _P`$, as $`Re=(l_x^3\widehat{f}_0)^{(1/2)}/\nu _T=676`$. The parameters concerning the polymers are $`\eta _p=\nu _p/\nu =.1`$, see e.g. , $`\rho _m^2/\rho _0^2=1000`$ and $`De=\tau (\widehat{f}_0/l_x)^{1/2}=.38`$. When comparing viscoelastic and Newtonian simulations, we keep constant the nominal Reynolds number, i.e. for the Newtonian case we use the same forcing intensity and a value for the viscosity equal to the value of $`\nu _T`$ for the viscoelastic fluid.

After the transient is elapsed, statistically stationary conditions are achieved. The general impression on the field may be gained by considering figure 1, where instantaneous vortical structures are shown for the Newtonian and the viscoelastic simulation. The eduction criterion for vortical structures is based on the discriminant of the velocity gradient tensor, one of the now classical approaches to extract vortical structures from a complex field, see e.g.  for details. Algebraically, the idea is to identify with vortical structures those regions where a couple of eigenvalues of the velocity gradient are complex conjugate, implying a locally helicoidal relative motion of nearby particles. In principle, other criteria could be employed and different field quantities could be visualized. In any case, the inspection of the instantaneous fields always leads to the conclusion that, with respect to Newtonian, much less but wider structures are present in the viscoelastic case.

Table I presents the global parameters of a typical viscoelastic simulation, $`De=.38`$, together with a corresponding Newtonian case. In the table, the Taylor microscale is defined as $`\lambda =(15/2u_{rms}^2/\mathrm{\Omega })^{1/2}`$, with $`\mathrm{\Omega }=1/2<\zeta ^2>`$, $`\zeta `$ being the vorticity. The average value of the free-energy is denoted by $`\overline{a}`$, while the other symbols have been already defined in the previous sections. Concerning the Kolmogorov scale, $`\eta `$, for the viscoelastic case this quantity is not uniquely defined. Actually, for Newtonian turbulence one assumes that the small scale dynamics is controlled by dissipation and viscosity, originating the classical definition of the Kolmogorov length as $`\eta =(\nu ^3/ϵ)^{1/4}`$. In the present case, several other parameters related to the more complex rheology may in principle affect the small scale dynamics. In the table we have conventionally assumed $`\eta _T=(\nu _T^3/ϵ_T)^{1/4}`$.

By discussing the data in table I, the macroscopic effects of viscoelasticity can be summarized as follows. For fixed forcing amplitude, i.e. for fixed nominal Reynolds number in identical boxes and with same total viscosity, $`\nu _T`$, we observe a substantial increase of the total dissipation. Clearly, in the steady state, this corresponds to an equivalent increase in the power extracted from the external forcing. For the viscoelastic case, most of the dissipation, about $`71\%`$, occurs in the polymers. The $`u_{rms}`$ increases of about $`20\%`$, with respect to the Newtonian simulation. The increase of the energy content in the flow is even more apparent if we consider that, in addition to kinetic energy, a non-negligible contribution arises from the free energy $`\overline{a}`$. In terms of total energy, $`\overline{E}=3/2u_{rms}^2+\overline{a}`$, the increase with respect to Newtonian turbulence is of the order of $`48\%`$. Clearly, the increase in the total dissipation for fixed total viscosity amounts to a reduction of the conventional Kolmogorov scale.

A further effect is observed when comparing the Taylor microscale for the two simulations. The increase in turbulent kinetic energy of about $`44\%`$ is accompanied by

a reduction of the order of $`50\%`$ in the enstrophy, giving an increase of $`73\%`$ for the Taylor microscale. This suggests a substantial depletion in the energy content of the high wave-number modes. Consistently a substantial increase in the Taylor-Reynolds number, $`Re_\lambda =u_{rms}\lambda /\nu _T`$, is achieved.

From the numerical point of view, a comparable value for $`Re_\lambda `$ in Newtonian fluids cannot be reached with only $`64`$ spatial modes. For the viscoelastic case, however, the requirements on the grid are less strict, since most of the dissipation, the polymeric contribution, is not associated to spatial gradients. We observe in passing, that spectral calculation with the FENE-P model may be prone to high-wavenumber instability, which is maintained under control by adding a small amount of artificial viscosity in the algorithm for the conformation tensor , see also  for an alternative approach.

The evolution of the system, followed for about $`50`$ large-eddy turnover times, $`T=(l_x/\widehat{f}_0)^{1/2}=7.9`$, is globally described in figure 2, where the time history of the instantaneous integral scale $`L=\pi /2_0^{K_{max}}1/k\widehat{e}(k)𝑑k/u_{rms}^2`$, where $`\widehat{e}`$ is the instantaneous energy spectrum, is shown on the top. From the comparison between viscoelastic and Newtonian case, the average increase in $`L`$ is apparent. Clearly this is associated to the increase in the amount of energy in the energy containing range, and corresponds to the observation of larger structures in the field. The presence of fewer but more energetic coherent structures is also consistent with the existence of more pronounced intermittency cycles, more evident in the history of the kinetic energy shown on the bottom part of the same figure.

Figure 3 describes the evolution of enstrophy, top, and Taylor microscale, bottom. On average, the enstrophy reduces with respect to Newtonian turbulence, and presents large spikes following with slight delay the corresponding peaks found in the kinetic energy. In figure 4 we report on top the time behavior of the spatial average of the free energy. Apparently, the amount of energy stored by the polymers is substantially lower than the kinetic energy, already given in figure 2 and reproduced

here for convenience. From the comparison one may note the more irregular behavior of the free-energy with reduced peak to peak oscillations. Finally, on the bottom of the figure we plot the evolution of the two forms of dissipation, Newtonian and viscoelastic respectively. Here the trend is opposite, with a larger portion of dissipation occurring in the polymers.

## VI Scale by scale budget

The Karman-Howarth equation (24) provides the appropriate tool to analyze the scale by scale budget for the velocity field. To ease its interpretation one should keep in mind the term

$`{\displaystyle \frac{1}{\pi r^2}}{\displaystyle \frac{d}{dt}}{\displaystyle _{B_r}}C_{ii}𝑑V_r,`$

originally present at the left-hand side of the equation, that has been dropped due to the assumed stationarity of the field.

A negative flux through the boundary of $`B_r`$,

$`\mathrm{\Phi }(r)={\displaystyle \frac{1}{4\pi r^2}}{\displaystyle _{B_r}}<\delta V^2\delta V_k+2T_{ki}^{}\delta V_i>n_k𝑑S_r,`$ (37)
corresponding to a flux towards the inside of $`B_r`$, implies a negative contribution to the integral up to $`r`$ of the correlation, i.e. an average depletion of correlation in the scales smaller than $`r`$. For stationary flows, the depletion balances the feeding operated by the forcing,

$`F(r)=4/3<f_iu_i>r=4/3\overline{ϵ}_Tr,`$

consistently with a steady correlation at all scales.

For polymers, the flux is split into two parts, the classical part which, for isotropic conditions, can be expressed in terms of the third order longitudinal structure function,

$`\mathrm{\Phi }_c(r)={\displaystyle \frac{1}{4\pi r^2}}{\displaystyle _{B_r}}<\delta V^2\delta V_k>n_k𝑑S_r={\displaystyle \frac{3}{5}}<\delta V_{}^3>`$ (38)
and the additional term associated with the viscoelastic reaction,

$`\mathrm{\Phi }_p(r)={\displaystyle \frac{1}{4\pi r^2}}{\displaystyle _{B_r}}<T_{ki}^{}\delta V_i>n_k𝑑S_r.`$ (40)
Figure 5 shows all the contributions appearing in the steady state Karman-Howarth equation, and, in particular, on the top, it reports the budget of the Newtonian simulation for comparison. Considering the Newtonian case first, the straight-line with negative slope corresponds to $`4/3\overline{ϵ}_Nr`$, filled squares. The positive correlation between forcing increments and velocity increments is represented by the filled diamonds. The nonlinear transfer term, $`\mathrm{\Phi }_c`$, is given by the filled circles, while the viscous correction corresponds to the filled triangles. All terms sum up to zero within the accuracy of the available statistics. Clearly, the viscous correction is present at all scales, confirming that the Reynolds number of the simulation is too small to achieve a proper inertial range. A further point to be highlighted is the non-negligible contribution of the correlation between forcing increments and velocity increments. Nonetheless, for a fair range of scales, up to $`r/\eta _T30`$, the balance involves only the nonlinear transfer term , $`\mathrm{\Phi }_c`$, and the forcing $`F`$. $`\mathrm{\Phi }_c`$ displays a behavior in terms of separation which, a part from minor corrections due to the other terms, follows the linear trend of $`F`$.

The same kind of budget is discussed for the viscoelastic case in the bottom part of the figure. The filled symbols have the same meaning defined for the Newtonian case, except that, now, the transfer term, filled circles, corresponds to the entire flux $`\mathrm{\Phi }`$, i.e. it accounts also for the polymeric contribution. As we may see, the picture grossly reproduces that already described for Newtonian turbulence. One should note that the forcing term, $`F`$, now involves the entire dissipation $`\overline{ϵ}_T`$, i.e. the sum of the

Newtonian and the polymeric contributions. The slope is then substantially larger than in the Newtonian case displayed on the top. Not unexpectedly, the global flux term is qualitatively similar to the corresponding Newtonian one, with a substantial range of scales where the linear trend defined by $`F`$ is reasonably well reproduced by $`\mathrm{\Phi }`$. This emerges quite clearly from the plots in figure 6 where the budgets for Newtonian and viscoelastic turbulence are superimposed the one to the other, after a rescaling proportional to the respective values of the global dissipation. Also for the polymers, we are in presence of a direct cascade, with a negative global flux.

Returning to the plots on the bottom of figure 5, the most interesting feature is represented by the open symbols, circles and squares, which give the splitting of the flux into the classical contribution associated to the third order velocity structure function, $`\mathrm{\Phi }_c`$ (circles), and to the polymer contribution, $`\mathrm{\Phi }_p`$ (triangles), respectively. The larger scales are less affected by the polymers, with $`\mathrm{\Phi }_p`$ becoming sub-leading with respect to $`\mathrm{\Phi }_c`$. One may conjecture that, at sufficiently large scale and for Reynolds numbers large enough, the turbulence scalings are substantially unaffected by the polymers. Under this assumption, we see that the total dissipation $`\overline{ϵ}_T`$ may be taken to measure a purely inertial energy transfer at large scales. This attaches to our conventional definition of Kolmogorov scale the obvious meaning of the length-scale where dissipation would occur in the absence of polymers.

Within the scenario according to which the large scales are unaffected by the polymers, the time criterion of Lumley  fixes the scale $`r^{}`$ below which the polymers begin to feel the turbulent fluctuations. This happens when the eddy-turnover time of classical turbulent eddies, i.e. without the effect of the viscoelastic reaction, equals the principal relaxation time of the polymer chains. Assuming Kolmogorov scaling for the velocity increments, this leads to

$`r^{}`$ $`=`$ $`(\overline{ϵ}_T\tau ^3)^{1/2}`$ (41)
or, scaled with the conventional Kolmogorov length, to

$`r^{}/\eta _T`$ $`=`$ $`\tau ^{3/2}(\overline{ϵ}_T/\nu _T)^{3/4}.`$ (42)
From the data given in table I, for our simulation we have $`r^{}/\eta _T70`$, within a computational box of $`l_x/\eta _T174`$, suggesting that the forcing, applied to the first shell of wavenumbers, acts just above the scale where the polymers already feel the turbulent fluctuations. In fact, a better way to apply the time criterion is to identify the eddy-turnover time as function of separation from the actual data provided by the DNS. This is done in figure 7, where as characteristic time for the scale $`r`$ we assume $`r/\mathrm{\Phi }_c^{1/3}`$. From the figure we infer that our viscoelastic turbulence achieve an eddy-turnover time comparable with the viscoelastic relaxation time, horizontal line in the figure, at a scale considerably smaller then predicted on purely dimensional grounds, leaving a reasonable range of between the forcing and Lumley scale.

Coming back to our comparison between the two components of the flux, at small scales we find that the component originated by the polymers takes over, entailing the reduction of the third-order structure function. The two curves, open circles and squares, clearly identify a cross-over scale $`\mathrm{}_p`$ according to the equation

$`\mathrm{\Phi }_p(\mathrm{}_p)=\mathrm{\Phi }_c(\mathrm{}_p).`$ (43)
For our case the balance is achieved at $`r/\eta _T22`$, while $`\lambda /\eta _T=33`$. Observe for comparison that, in the corresponding Newtonian case, $`\lambda /\eta `$ is about half the value. By inspection of the figure, we also find that all along the available range of scales the polymers reaction is dynamically relevant for the turbulent fluctuations, with $`\mathrm{\Phi }_p`$ always comparable in magnitude with $`\mathrm{\Phi }_c`$.

Dimensionally, the right-hand side of eq. (43) can be estimated as $`\overline{ϵ}_T\mathrm{}_p`$, as follows by roughly assuming a purely classical inertial scaling for $`r\mathrm{}_p`$. On the other hand, $`\mathrm{\Phi }_p(\mathrm{}_p)`$ should be somehow proportional to $`\sigma ^{}\nu _p/\tau \overline{ϵ}_T^{1/3}\mathrm{}_p^{1/3}`$, where, again, Kolmogorov scaling is used to estimate $`\delta V`$ in eq. (40). Accounting for the constitutive relation (2), $`T^{}`$ is expressed in terms of the, as yet, undetermined factor $`\sigma ^{}`$ which should account for the level of stretching achieved in the polymers at the cross-over scale. This yields the expression

$`\mathrm{}_p={\displaystyle \frac{(\sigma ^{}\nu _p/\tau )^{3/2}}{\overline{ϵ}_T}},`$

which can be used to build the dimensionless quantity

$`\sigma ^{}`$ $`=`$ $`\tau /\nu _p(\mathrm{}_p\overline{ϵ}_T)^{2/3}.`$ (44)
For the present case we have $`(\nu _p/\tau )^{3/2}/\overline{ϵ}_T\mathrm{1.4\hspace{0.17em}10}^5`$, which, after considering that $`\mathrm{}_p/\eta _p22`$,

yields $`\sigma ^{}1400`$. The parameter $`\sigma ^{}`$ allows to quantify the effective stretching of the chains at cross-over by assuming $`\sigma ^{}f^{}R_{kk}^{}/\rho _0^2`$, see eq. (2).

The quantity

$`\sigma =f(R_{kk};\rho _{max},\rho _0)R_{kk}/\rho _0^2`$ (45)
is plotted versus the elongation $`R_{kk}/\rho _0^2`$ on the right of figure 8. We see clearly the finite extensibility effect, with the unbound growth of $`\sigma `$ as the critical extension $`\rho _{max}^2/\rho _0^2`$ is approached from the left. Note that, at small elongation, $`\sigma `$ goes almost linearly with $`R_{kk}/\rho _0^2`$. The horizontal line plotted in the figure corresponds to the value $`\sigma ^{}`$ we have estimated from the numerical simulation. The intercept of this line with the curve representing $`\sigma `$ lays relatively far from the critical extension, though well within the region of non-linear elastic response.

It may be instructive to compare the elongation associated to the cross-over scale to the pdf of $`R_{kk}/\rho _0^2`$.

To this purpose, we superimpose to the plot of $`\sigma `$ given in figure 8 the histogram of the elongation as obtained by the numerical simulation. The effective elongation at cross-over, denoted by the vertical line, lies to the right of the maximum of the histogram, indicating a level of stretching substantially larger then either the most probable or the average elongation. Clearly the estimate here proposed should be understood as an order of magnitude analysis, to allow the comparison of different flow conditions.

To complete the analysis, in figure 9 we address the Yaglom-like equation for the free-energy. As seen in the top part of the figure, the left and the right hand side of eq. (34) balance within the accuracy of the available statistics.

The fluctuations of excess-power are converted into a free-energy flux with a negative sign, as wee see from the figure. Hence, also for the microstructure, we are in presence of a direct cascade. We like to stress that the direction of the cascade is not a priori obvious in this case, since we miss the standard term, linear in the

separation and related to the dissipation of the scalar, which usually determines the sign of the energy flux.

From the figure, the free-energy flux,

$`\mathrm{\Phi }_a=<\delta a^2\delta V_{}>,`$

manifests a pronounced minimum in correspondence with the cross-over scale $`\mathrm{}_p`$, suggesting that, near cross-over the activity of the polymers is particularly intense and the polymers tend to become less active as the scale is increased. The bottom part of the figure gives the decomposition of the right-hand side of equation (34), with the positive term

associated to the stress power denoted by filled squares and that associated with the polymeric dissipation indicated by the filled triangles.

Order of magnitude considerations confirm that the data shown in the bottom part of figures 5 and in figure 9 are consistent. At cross-over, $`r=\mathrm{}_p`$, we find from figure 5 that $`\mathrm{\Phi }_c^{}<\delta V_{}^{}{}_{}{}^{3}>.12`$, suggesting $`\delta V^{}.5`$. From figure 9 we also have $`<\delta a_{}^{}{}_{}{}^{2}\delta V^{}>.00015`$. It follows $`\delta a^{}.017`$ a value considerably close to $`\delta V_{}^{}{}_{}{}^{2}.025`$. In other words, at the cross-over scale fluctuations of free-energy, or more loosely of elastic energy, and fluctuation of turbulent kinetic energy become comparable, as it should be expected.

## VII Concluding remarks

The interaction between long chain polymers in dilute solution and turbulence is known to have a profound influence on the structure of turbulence. Several possible explanations have been proposed, with the eventual aim to understand the drag reducing properties manifested by such fluids in wall bounded turbulent flows. Even the most powerful approach, however, has to face the elusive mechanism by which polymers dynamics and turbulence are mutually coupled. Among others, specially noteworthy is the theory of De Gennes and Tabor, who proposed a cascade model for the viscoelastic turbulence of dilute polymers solutions. Their elegant theory issues from considering homogeneous isotropic conditions and rests upon suitable assumptions on the self-similar stretching of polymers chains under a fluctuating velocity field. The entire model is essentially phenomenological in nature, and should thus be contrasted with experimental observations.

In this context, we believe to have provided here a reasonable framework based on two essential ingredients. The first ingredient is the direct numerical simulation of the natural flow where the cascade issue should be addressed, namely homogeneous isotropic turbulence numerically modeled as a triply-periodic box with random forcing at large scales. Beside the obvious implementation in terms of spectral algorithms, this step required the selection of a good rheological model for the polymers. Based on previous experience, our natural choice has been the FENE-P model, known to capture reasonably well the drag-reducing effects still maintaining an affordable level of computational complexity.

The second ingredient is the extension to viscoelastic fluids of two of the, let’s say, few exact results known in classical turbulence theory, namely the Karman-Howarth equation for the velocity and the Yaglom equation for scalars. As a characteristic feature of viscoelastic turbulence, the coupling between velocity and extra-stress in the Karman-Howarth equation prevents recasting the result in the familiar form involving only longitudinal structure functions. Hence we have adopted a natural generalization in terms of suitable surface and volume integrals in the space of separations. Concerning the microstructure, the descriptor field is a second order tensor, making the derivation of the appropriate equation exceedingly cumbersome. We decided then to fucus our attention on the single most important scalar quantity, namely the free-energy, whose Yaglom-like equation follows straightforwardly.

The scale by scale budget based on the Karman-Howarth and Yaglom equation has been used to understand the respective role of the traditional energy transfer term versus polymeric transfer. We find a direct cascade occurring both in the kinetic field and in the microstructure, for which, as already discussed, the sign of the energy flux is not a priori obvious.

The budget has allowed us to identify a cross-over scale below which the polymers contribution becomes dominant, being sub-leading with respect to inertial transfer at larger scales. This is the analogous of the scale defining the elastic limit in the theory of De Gennes and Tabor. It is here derived in a rigorous context, and evaluated on the basis of a direct numerical simulation of a drag-reducing visco-elastic material as described by the FENE-P model.

We find that the polymers substantially deplete the energy content of the smallest scales, in agreement with recent experimental results, but they affect the dynamics of the fluctuations at all the scales of our simulation. The latter effect could be interpreted as an artifact of the limited extension of our computational domain which forced us to apply the external excitation relatively close to the scale defined by Lumley’s time criterion. Nonetheless the present results make questionable the existence of a purely passive range, where polymers are deformed without back-reaction on the velocity field.

Based on the balance of terms appearing in the exact form of the Karman-Howarth equation we have proposed a dimensionless parameter which may give a quantitative measure of the stretching the polymers experience at cross-over, when elastic effects become comparable with the corresponding kinetic energy.

Finally, a substantial level of intermittency is observed in the system, with fluctuations in turbulent kinetic energy considerably larger than usually observed in homogeneous isotropic turbulence of Newtonian fluids. This effect is presumably associated with the existence of the additional mechanism of energy removal from the scales of the forcing. In particular, a substantial part of the energy income does not follow the classical route cascading towards viscous dissipation. Instead, it is moved to the microstructure, to feed an additional cascade process seemingly requiring a larger level of fluctuations.