format: text

text:
# Formal Solution of a Class of Reaction-Diffusion Models: Reduction to a Single-Particle Problem

## Abstract

We consider the trapping reaction $`A+BB`$ in space dimension $`d2`$. By formally eliminating the $`B`$ particles from the problem we derive an effective dynamics for the $`A`$ particles from which the survival probability of a given $`A`$ particle and the statistics of its spatial fluctuations can be calculated in a rather general way. The method can be extended to the study of annihilation/coalescence reactions, $`B+B0`$ or $`B`$, in $`d=2`$.

First-passage problems involving more than a few degrees of freedom are notoriously difficult to solve Sid ; SatyaRev . In this Letter we introduce a technique that enables one to solve a class of first-passage problems involving an infinite number of degrees of freedom. For definiteness we develop the method in the context of the “trapping reaction”, $`A+BB`$, but the applications are more general, as emphasized in the latter part of the Letter. The main result of our approach is to reduce the problem to one described by a single degree of freedom whose late-time behavior can be extracted analytically.

The asymptotic dynamics of the trapping reaction has been a long-standing puzzle. The main question is how the density of $`A`$ particles decreases with time. A related problem, much studied in the context of chemical kinetics TW ; BL  is the two-species annihilation reaction, $`A+B0`$, with initial densities $`\rho _A(0)<\rho _B(0)`$. This is equivalent to the trapping reaction at late times, when $`\rho _A(t)\rho _B(t)`$ and $`\rho _B(t)`$ is essentially constant. Again, the standard problem is to compute the asymptotic form of the $`A`$-particle density $`\rho _A(t)`$ or, equivalently, the probability $`Q(t)`$ that a given $`A`$ particle survives until time $`t`$. Since the particles do not interact with other particles of the same species, to compute $`Q(t)`$ it suffices to consider a single $`A`$ particle moving in an infinite sea of $`B`$ particles with density $`\rho `$ ($`=\rho _B`$).

Traditional approaches to this type of problem treat the $`A`$ particle as an absorbing boundary for the $`B`$ particles. Unfortunately, for an arbitrary $`A`$-particle trajectory this absorbing-boundary problem cannot be solved. In this Letter we introduce a different approach in which we treat the $`A`$ and $`B`$ particles as if they were non-interacting. We exploit the initial condition that each $`B`$ particle is randomly located anywhere in the system to show that certain ‘events’, where a $`B`$ particle meets the $`A`$ particle for the first time (remember that we are treating them as non-interacting, so they can meet more than once) have a Poisson distribution, i.e. the probability $`p_n`$ that $`n`$ such events have occurred up to time $`t`$ is given by $`p_n=(\mu ^n/n!)\mathrm{exp}(\mu )`$, where the mean, $`\mu `$, of the distribution is a functional, $`\mu [\stackrel{}{z}]`$ of the trajectory $`\stackrel{}{z}(\tau )`$, $`0\tau t`$, of the $`A`$ particle. The probability that the trajectory $`\stackrel{}{z}(\tau )`$ has survived, in the original interacting problem, is simply $`p_0[\stackrel{}{z}]=\mathrm{exp}(\mu [\stackrel{}{z}])`$. Finally, $`Q(t)`$ is obtained by averaging $`\mathrm{exp}(\mu [\stackrel{}{z}])`$ over all possible $`A`$-particle trajectories $`\stackrel{}{z}(\tau )`$ with the appropriate (Wiener) measure, $`\mathrm{exp}[(1/4D^{})_0^t𝑑\tau (d\stackrel{}{z}/d\tau )^2]`$, where $`D^{}`$ is the $`A`$ particle’s diffusion constant. In this way, the $`B`$ particles have been eliminated from the problem, and one has an effective $`A`$-particle dynamics described by the Wiener measure and the functional $`\mu [\stackrel{}{z}]`$. The final step, that makes further analytical progress possible, is the observation that the path integral over $`\stackrel{}{z}(\tau )`$ is dominated at late times by a single $`A`$-particle trajectory.

The main results of this approach are: (i) the trajectory where the $`A`$ particle is stationary is proved to be the dominant trajectory and determines the asymptotic form of the $`A`$ particle’s survival probability BL , $`Q(t)\mathrm{exp}(\lambda _dt^{d/2})`$ for $`d<2`$ (with a logarithmic correction in $`d=2`$), where $`\lambda _d`$ is a calculable constant and $`d`$ is the dimensionality of space; (ii) typical fluctuations of surviving $`A`$-particle trajectories around this dominant path have variance $`z^2(t)t^{(2d)/2}`$ for $`d<2`$ (iii) exact results are obtained for $`Q(t)`$ and the form of the dominant path in a system with a non-uniform initial density of $`B`$ particles; (iv) this approach provides a powerful method for calculating first-passage properties for a deterministically moving boundary $`\stackrel{}{z}(t)`$; (v) the method also provides a formalism for calculating $`Q(t)`$ in the highly nontrivial situation where the $`B`$ particles themselves interact, e.g. $`B+B0`$, at least in $`d=2`$ where the density correlations induced by these reactions are negligible.

We begin by deriving the Poisson property that plays a central role in the analysis. We consider a finite volume, $`V`$, containing $`N=\rho V`$ $`B`$-particles (diffusion constant $`D`$), randomly distributed within it, and a single $`A`$ particle (diffusion constant $`D^{}`$), initially located at the origin. Let $`\stackrel{}{z}(t)`$ be the $`A`$ particle’s trajectory, and let $`P(\stackrel{}{x},t)`$ be the probability that a given $`B`$ particle, starting at $`\stackrel{}{x}`$, has met the $`A`$ particle before time $`t`$. The average of this quantity over the initial position, $`\stackrel{}{x}`$, is $`(1/V)_V𝑑VP(\stackrel{}{x},t)=R(t)/V`$, where $`R(t)`$ is an implicit functional of $`\stackrel{}{z}(t)`$. The probability that $`n`$ distinct $`B`$ particles have met the $`A`$ particle, averaged over their initial positions, is $`p_n(t)=\left(\genfrac{}{}{0pt}{}{N}{n}\right)(R/V)^n(1R/V)^{Nn}`$. Taking the limit $`N\mathrm{}`$, $`V\mathrm{}`$, with $`\rho =N/V`$ and $`n`$ held fixed, one recovers the Poisson distribution, $`p_n=(\mu ^n/n!)\mathrm{exp}(\mu )`$, with $`\mu =\rho R`$.

One can derive an equation for the functional $`\mu [z]`$ by calculating, in two ways, the probability density to find a $`B`$ particle at the point $`\stackrel{}{z}(t)`$ at time $`t`$. First, since the particles are treated as non-interacting, and the $`B`$ particles start in a steady-state configuration of uniform density $`\rho `$, this probability density is just $`\rho `$. Secondly, from the Poisson property, the probability that a $`B`$ particle (i.e. any $`B`$ particle) meets the $`A`$ particle for the first time in the time interval $`(t^{},t^{}+dt^{})`$ is $`\dot{\mu }(t^{})dt^{}`$. The probability density for such a particle to subsequently arrive at $`\stackrel{}{z}(t)`$ at time $`t`$ is given by the diffusion propagator $`G(\stackrel{}{z}(t),t|\stackrel{}{z}(t^{}),t^{})=[4\pi D(tt^{})]^{d/2}\mathrm{exp}\{[\stackrel{}{z}(t)\stackrel{}{z}(t^{})]^2/4D(tt^{})\}`$. Equating the results from these two methods gives our fundamental equation,

$$\rho =_0^t𝑑t^{}\dot{\mu }(t^{})G(\stackrel{}{z}(t),t|\stackrel{}{z}(t^{}),t^{}),$$
(1)
which is an implicit equation for the functional $`\mu [\stackrel{}{z}]`$ (noting that $`\mu (t=0)=0`$, since no $`B`$ particle can meet the $`A`$ particle in zero time). Finally, $`Q(t)=\mathrm{exp}(\mu [\stackrel{}{z}])_z`$, where the average is over all paths $`\stackrel{}{z}(t)`$ weighted with the Wiener measure.

As a first application of this equation we prove that the trajectory $`\stackrel{}{z}=0`$ is the dominant path, i.e. that it gives the smallest possible value of $`\mu [\stackrel{}{z}]`$ for all $`t`$. This function, $`\mu _0(t)`$, satisfies Eq. (1) with $`\stackrel{}{z}=0`$:

$$\rho =_0^t𝑑t^{}\dot{\mu }_0(t^{})[4\pi D(tt^{})]^{d/2}.$$
(2)
By inspection, $`\mu _0(t)`$ must have the form $`\mu _0(t)=\lambda _dt^{d/2}`$ (for $`d<2`$), in order that the right-hand side be independent of $`t`$. Substituting this form in (2), and evaluating the integral, gives

$$\lambda _d=\rho \left(\frac{2}{\pi d}\right)\mathrm{sin}\left(\frac{\pi d}{2}\right)\left(4\pi D\right)^{d/2},d<2,$$
(3)
while for $`d=2`$ one finds, for $`t\mathrm{}`$, $`\mu _0(t)4\pi \rho Dt/\mathrm{ln}t`$ d=2 . The corresponding A-particle survival probability is $`Q_0(t)=\mathrm{exp}[\mu _0(t)]`$. This simple case of a static $`A`$ particle is sometimes called the ‘target annihilation problem’, and our method reproduces the known results for that problem target  in a very simple way. To prove that $`\stackrel{}{z}(t)=0`$ gives the global minimum of $`\mu [\stackrel{}{z}]`$ we write $`\mu =\mu _0+\mu _1`$ in (1). This equation can then be rearranged, with the help of Laplace transform techniques, to give an implicit equation for $`\mu _1[\stackrel{}{z}]`$:

$`\mu _1[\stackrel{}{z}]`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}(\pi d/2)}{\pi }}{\displaystyle _0^t}{\displaystyle \frac{dt_1}{(tt_1)^{(2d)/2}}}`$ (4)
$`\times {\displaystyle _0^{t_1}}{\displaystyle \frac{dt_2}{(t_1t_2)^{d/2}}}\dot{\mu }(t_2)K(t_1,t_2),`$
where $`K(t_1,t_2)=1\mathrm{exp}\{[\stackrel{}{z}(t_1)\stackrel{}{z}(t_2)]^2/4D(t_1t_2)\}`$.

Eq. (4) is ‘implicit’ because the full $`\mu `$ appears on the right-hand side. Now note that $`K(t_1,t_2)0`$ and $`\dot{\mu }0`$ (because $`\mu (t)`$ is the mean number of different $`B`$ particles that have met the $`A`$ particle up to time $`t`$ – clearly a non-decreasing function). Therefore $`\mu _1[\stackrel{}{z}]0`$ for all paths $`\stackrel{}{z}(t)`$, with equality when $`\stackrel{}{z}(t)=0`$ for all $`t`$. It follows that $`Q(t)\mathrm{exp}(\mu _0\mu _1)_\stackrel{}{z}\mathrm{exp}[\mu _0(t)]`$. This rigorous upper bound for $`Q(t)`$, combined with the identical rigorous lower bound derived in BB , proves that the asymptotic form of $`Q(t)`$ is the same as for the ‘target’ problem, where the $`A`$ particle is stationary, for all $`d2`$. The interpretation of this result is that, since $`\mu `$ is large for $`t\mathrm{}`$ ($`\mu t^{d/2}`$), the path-integral for $`Q(t)`$ is dominated by the path that minimizes $`\mu `$, i.e. we are essentially evaluating the path integral by the method of steepest descents. Small fluctuations around the dominant path will determine the corrections to the asymptotic form.

We next compute the probability distribution, $`P(z,t)`$, of the position $`z`$ of the $`A`$ particle at time $`t`$, given that it survives. Numerical studies MG  suggest that, in $`d=1`$, $`z^2(t)^{1/2}t^\varphi `$, with $`\varphi =0.250.3`$, while similar studies in $`d=2`$ are inconclusive. Our methods give $`\varphi =1/4`$ in $`d=1`$ and $`\varphi =(2d)/4`$ for all $`d<2`$. The technique is to expand $`\mu _1[\stackrel{}{z}]`$, given by Eq. (4), to order $`z^2`$ to compute the variance of the Gaussian fluctuations around the dominant trajectory $`z(t)=0`$. To this order one can replace $`\mu (t_2)`$ on the right-hand side by $`\mu _0(t_2)=\lambda _dt_2^{d/2}`$, and expand the function $`K(t_1,t_2)`$ to order $`z^2`$. Specializing to $`d=1`$, the result is, at time $`t`$,

$$\mu _1[z]=\frac{\lambda _1}{8\pi D}_0^t\frac{dt_1}{\sqrt{tt_1}}_0^{t_1}𝑑t_2\frac{[z(t_1)z(t_2)]^2}{\sqrt{t_2}(t_1t_2)^{3/2}}.$$
(5)
The probability distribution for $`z`$ at time $`t`$ is given, up to an overall normalization, by the path integral

$$P(z,t)=𝒟z(t)\mathrm{exp}\left(\frac{1}{4D^{}}_0^t𝑑\tau \dot{z}^2(\tau )\mu _1[z]\right),$$
(6)
where the integral is over all paths satisfying the boundary conditions $`z(0)=0`$, $`z(t)=z`$.

The path integral has the form $`𝒟z(t)\mathrm{exp}(S[z])`$, where the ‘action’ $`S[z]=S_T[x]+S_V[z]`$ is a quadratic functional of $`z(\tau )`$, where $`S_T[z]=(1/4D^{})_0^t𝑑\tau \dot{z}^2(\tau )`$, and $`S_V[z]=\mu _1[z]`$. Since the integrand is Gaussian, the $`z`$-dependence of the path integral is exactly captured by the extremal path connecting $`z(0)=0`$ and $`z(t)=z`$. Although we have been unable to find this path analytically, power counting on the two terms suggests $`S_V\lambda _1z^2/D\sqrt{t}`$ and $`S_Tz^2/D^{}t`$, so that $`S_V`$ dominates at large $`t`$. A more careful analysis details  confirms the dominance of $`S_V`$ and the form $`S_V=c_V\lambda _1z^2/D\sqrt{t}`$, where $`c_V`$ is a pure number. The probability weight for the fluctuations $`z(t)`$ of surviving trajectories is therefore Gaussian, with variance $`z^2(t)=(D/2c_V\lambda _1)\sqrt{t}`$ for large $`t`$, with the value of $`c_V`$ determined by the extremal path. Note that this result is independent of the $`A`$-particle diffusion constant $`D^{}`$. If $`D^{}=0`$, of course, there are no fluctuations, but for any $`D^{}>0`$ the variance becomes independent of $`D^{}`$ at large enough $`t`$. Also, since typically $`zt^{1/4}`$, the next term in the expansion of the function $`K(t_1,t_2)`$ in (4) is of relative order $`z^2/Dtt^{1/2}`$ so it is negligible at large $`t`$. The distribution of $`z`$ is therefore exactly Gaussian, at least in the ‘scaling limit’ $`z\mathrm{}`$, $`t\mathrm{}`$, with $`z/t^{1/4}`$ fixed. Similar arguments details  give the generalization $`z^2(t)(D/\lambda _d)t^{(2d)/2}`$ for $`d<2`$.

We turn now to a related problem with a non-trivial dominant path that can be exactly determined. Consider, in $`d=1`$, a system where the density of $`B`$ particles at $`t=0`$ has different values, $`\rho _L`$ and $`\rho _R`$, to the left and right of the $`A`$ particle. The derivation of an equation for $`\mu [z]`$ proceeds exactly as before, except that the probability density to find a $`B`$ particle at the point $`z`$ at time $`t`$ in the noninteracting system, which appears on the left of Eq. (1), has to be recalculated. In terms of the diffusion propagator $`G`$ introduced earlier, this probability is (quite generally) $`P_B(z,t)=_{\mathrm{}}^{\mathrm{}}𝑑x\rho (x)G(z,t|x,0)`$ where $`\rho (x)`$ is the initial $`B`$-particle density at position $`x`$. When $`\rho (x)=\rho `$, a constant, one finds $`P_B(z,t)=\rho `$, as before. When $`\rho (x)=\rho _L`$ for $`x<0`$ and $`\rho _R`$ for $`x>0`$, the generalized version of (1) becomes

$$\rho \left[1\mathrm{\Delta }\mathrm{erf}\left(\frac{z(t)}{\sqrt{4Dt}}\right)\right]=_0^t𝑑t^{}\dot{\mu }(t^{})G(z(t),t|z(t^{}),t^{}),$$
(7)
where $`\rho =(\rho _L+\rho _R)/2`$ is now the mean density, $`\mathrm{\Delta }=(\rho _L\rho _R)/(\rho _L+\rho _R)`$ is a measure of the left-right asymmetry, and $`\mathrm{erf}(x)`$ is the error function.

Physical intuition suggests that, because of the asymmetry, surviving $`A`$-particle trajectories will tend to be those that drift into the region (the right, say) where the $`B`$-particle density is initially smaller. Upper and lower bounds have been derived earlier BB2  for the asymptotics of the $`A`$-particle survival probability, $`Q(t)`$, which show that it has the asymptotic form $`Q(t)\mathrm{exp}[g(\mathrm{\Delta })\lambda _1\sqrt{t}]`$, where $`g(0)=1`$ for consistency with the symmetric case, $`\rho _L=\rho _R`$. This form for $`Q(t)`$ shows that $`\mu [z]`$ for the optimal path has the time-dependence $`\mu \sqrt{t}`$. Both sides of (7) can then be rendered time-independent by the choice $`z(\tau )=\alpha \sqrt{4D\tau }`$ for all $`\tau t`$, where $`\alpha `$ is a constant to be determined. A more detailed analysis details  shows that the dominant path is indeed of this form.

Putting $`\mu (t)=g(\mathrm{\Delta })\lambda _1\sqrt{t}`$ and $`z(\tau )=\alpha \sqrt{4D\tau }`$, in (7), and evaluating the integral on the right-hand side gives

$$g(\mathrm{\Delta })=\frac{\mathrm{exp}(\alpha ^2)}{1\mathrm{erf}^2(\alpha )}[1\mathrm{\Delta }\mathrm{erf}(\alpha )].$$
(8)
The final step is to minimize the right-hand side with respect to $`\alpha `$ to obtain the optimal path. This can be done numerically. The resulting $`g(\mathrm{\Delta })`$ is shown for $`\mathrm{\Delta }0`$ as the inset in Figure 1 (note that, by symmetry, $`g(\mathrm{\Delta })`$ is symmetric around $`\mathrm{\Delta }=0`$). It clearly satisfies the bounds $`1|\mathrm{\Delta }|g(\mathrm{\Delta })1`$ derived in BB2 .

In the main part of Fig. 1, numerical results for $`g(t)=\mathrm{ln}Q(t)/\lambda _1\sqrt{t}`$, obtained using the algorithm of ref.MG , are displayed, with $`D=1/2`$ and $`\rho _L=0.5`$ in all cases, while $`\rho _R=0.25,0.125,0.0625`$ (top to bottom on the right), corresponding to $`\mathrm{\Delta }=1/3,3/5,7/9`$ respectively. The horizontal lines on the right show the asymptotic values obtained from the inset for the corresponding values of $`\mathrm{\Delta }`$. The slow approach to asymptopia is similar to that observed MG ; BB2  in the symmetric case ($`\mathrm{\Delta }=0`$).

As a bonus, the same calculation solves the first-passage problem of $`B`$ particles with density $`\rho _L`$ for $`x<0`$ and $`\rho _R`$ for $`x>0`$ moving in the presence of a deterministically moving absorbing boundary located at $`x(t)=\alpha \sqrt{4Dt}`$. The probability that no particle has reached the boundary up to time $`t`$ has the simple form $`Q(t)=\mathrm{exp}[g(\mathrm{\Delta })\lambda _1\sqrt{t}]`$, with $`g(\mathrm{\Delta })`$ given by (8). We are not aware of any other way of obtaining this result. Extensions to deterministically moving absorbing boundaries in dimension $`d>1`$ are also possible details .

In the final part of this Letter, we apply this approach to a nontrivial problem with $`d=2`$. Consider the annihilation/coalescence reaction $`B+B0`$ with probability $`1/(q1)`$ and $`B+BB`$ with probability $`(q2)/(q1)`$. The density of $`B`$ particles is known to decay as $`\rho (t)=a_d[(q1)/q](Dt)^{d/2}`$ for $`d<2`$, where $`a_d`$ is a universal constant Lee ; Cardy  equal to $`1/2\pi ϵ`$ for $`d2`$, with $`ϵ=2d`$. For $`d=2`$, a logarithmic correction is obtained, $`\rho (t)\mathrm{ln}t/t`$. Now suppose that one of the $`B`$ particles is tagged and relabeled an $`A`$ particle, with diffusion constant $`D^{}`$. We consider the probability (the ‘walker persistence’ probability WP ), $`Q(t)`$, that the $`A`$ particle has not met any $`B`$ particle up to time $`t`$. The limit $`d2`$ provides a simplification because it is the borderline dimension above which the rate equation approach, $`d\rho /dt\rho ^2`$, which gives $`\rho (t)1/t`$, is qualitatively correct because density fluctuations can be ignored Lee ; Cardy . Eq.(1) is readily adapted to this case. As before, we treat the $`A`$ particle as if it does not interact with the $`B`$ particles, while the interactions of the $`B`$ particles with each other give rise to their decreasing density, $`\rho (t)`$. The Poisson distribution for the number of first crossings of the $`A`$ particle by $`B`$ particles still holds for this system. The left-hand side of Eq. (1), i.e. the probability density to find a $`B`$ particle at the point $`z(t)`$ at time $`t`$, becomes $`\rho (t)`$, while on the right-hand side the propagator $`G(z(t),t|z(t^{}),t^{})`$ has to be multiplied by a factor $`\rho (t)/\rho (t^{})`$, being the probability of a given $`B`$ particle surviving till time $`t`$ given that it survives till time $`t^{}`$. The required generalization of (1) then reads

$$1=_0^t\frac{dt^{}}{\rho (t^{})}\dot{\mu }(t^{})G(\stackrel{}{z}(t),t|\stackrel{}{z}(t^{}),t^{}).$$
(9)
It is convenient to approach the limit $`d2`$ from below. Consider first the case $`D^{}=0`$, for which $`Q(t)`$ becomes the ‘site-persistence’ probability, i.e. the probability that a given point in space (the location of the $`A`$ particle) has not been visited by any $`B`$ particle. The static $`A`$ particle corresponds to $`z(\tau )=0`$, for all $`\tau `$, and with $`\rho (t^{})=a_d[(q1)/q](Dt^{})^{d/2}`$ for large $`t`$, (9) becomes

$$(4\pi )^{d/2}a_d(q1)/q=_0^t𝑑t^{}\dot{\mu }(t^{})t_{}^{}{}_{}{}^{d/2}(tt^{})^{d/2}$$
(10)
for large $`t`$. In order that the right-hand side be time-independent for large $`t`$, $`\mu (t)`$ must have the asymptotic form $`\mu (t)\theta \mathrm{ln}t`$. Inserting this form into (10), and evaluating the integral gives $`\theta =2^d\pi ^{d/21}\mathrm{sin}(\pi ϵ/2)a_d(q1)/q`$. Taking the limit $`ϵ0`$, using $`a_d=1/2\pi ϵ`$ in this limit, gives $`\theta =(q1)/q`$ for $`d=2`$. Finally, $`Q(t)=\mathrm{exp}[\mu (t)]t^\theta `$ for $`t\mathrm{}`$. The result $`\theta =(q1)/q`$ in $`d=2`$ agrees with that obtained by Cardy Cardy  using field-theoretic methods.

If one now considers the case $`D^{}>0`$, i.e. a diffusing $`A`$ particle, one sees immediately that for $`D^{}=D`$ and $`q=2`$ (so that $`B+B0`$ always), the $`A`$ particle is equivalent to another $`B`$ particle, and $`Q(t)=\rho (t)`$, the density, since every surviving particle has not met any other particle. Hence $`Q(t)=\rho (t)\mathrm{ln}t/t`$ for $`D^{}=D`$ and $`q=2`$. This suggests that, for general $`D^{}`$, $`Q(t)`$ will decay as $`t^\theta `$ where $`\theta `$ is a nontrivial function of $`D^{}/D`$.

The calculation of $`\theta `$ can readily be extended to $`D^{}>0`$ within the present formalism, using a power-series expansion in $`D^{}/D`$. The method is to use a cumulant expansion to write $`Q(t)=\mathrm{exp}(\mu )_z=\mathrm{exp}(\mu _z+\{\mu ^2_z\mu _z^2\}/2!+\mathrm{})`$. To first order in $`D^{}`$, one needs only the first cumulant. The result is $`\theta =(1+D^{}/D)(q1)/q`$. This class of problems, where $`Q(t)`$ decays as a power-law instead of a stretched exponential, marks the borderline where the path-integral for $`Q(t)`$ is no longer dominated by a single path (and small fluctuations about it), giving a leading large-$`t`$ form independent of $`D^{}`$, but has to be evaluated exactly, with results that depend on $`D^{}`$ even for $`t\mathrm{}`$. Full details of this calculation, together with results to higher order in $`D^{}/D`$, will be presented elsewhere.

In conclusion, we have introduced an analytic approach to a class of reaction-diffusion models that reduces them to one-particle systems. We hope to use this method in future to address, inter alia, the vexing problem of the very slow approach to asymptopia in the trapping reaction.

RAB thanks the EPSRC for financial support under grant GR/R53197.


text:
# Abstract

#### Abstract

: Recently the Sudbury Neutrino Observatory (SNO) has measured both the rates as well as the day-night asymmetries in the flux of neutrinos from the Sun in the charged-current (CC), neutral-current (NC) and elastic scattering (ES) channels. Motivated by the small but negative day-night asymmetry in the NC channel at SNO (1.2 standard deviations away from zero), we consider whether such a non-zero asymmetry can be generated in a framework where active neutrinos can both oscillate and decay, since the NC day-night asymmetry is zero in a frame-work that includes pure oscillations of active flavours. We find that the neutral current day-night asymmetry is strictly zero when we include both oscillations and decay. The result holds for arbitrary number of active neutrino species, with some reasonable assumptions on the decay parameters. Hence, a non-zero asymmetry in the neutral current sector, if it survives, can only come from oscillations/decay into sterile flavours. We also examine the allowed region of parameter space coming from a combined analysis of the total event rates in the charged-current, neutral current and elastic scattering sectors in a two-flavour framework, and find that the neutral current rates are most sensitive to (and hence restrict considerably) the allowed values of the decay parameter involving the life-time of the heavier mass eigenstate to $`\tau _{2_0}>10^4`$ s for eV mass neutrinos.

PACS numbers: 14.60.Pq, 26.65.+t

## 1 Introduction

There is now very good evidence for neutrino flavour oscillation from solar neutrino data . Combined analyses of these data seem to indicate  a preference for the so-called Large Mixing Angle (LMA) solution in the mass-squared difference–mixing angle ($`\delta `$-$`\theta `$) parameter space (in the case of two-flavour mixing); on inclusion of a third flavour, these results still apply to the 1–2 mixing angle and corresponding mass-squared difference ($`\delta _{12}`$-$`\omega `$), $`\omega \theta _{12}`$, since the 1–3 mixing angle ($`\varphi \theta _{13}`$) is known to be small .) This LMA solution is sensitive to solar matter effects on the propagating neutrino in a mildly energy-dependent way. The Small Mixing Angle (SMA) solution which is also a viable solution  to the solar neutrino problem, on the other hand, has significant energy dependence due to solar matter effects. A detailed study of the various available data on solar neutrinos in order to isolate the allowed parameter regions in oscillation parameter space has been the subject of many rigourous analyses . Analyses of data including decay have also been done for example in Refs. . Detailed analysis of models including oscillation and decay, applicable to atmospheric and supernova neutrinos as well, can be found in Ref. . In this paper, we concentrate on some constraints specifically arising from solar neutrino data (both from the rates and day-night asymmetries), particularly from the neutral current (NC) sector, at the Sudbury Neutrino Observatory (SNO) .

Information on Earth matter effects on neutrino mixing and hence flavour survival and conversion probabilities is available through a study of the day-night asymmetry, $`A^{ND}`$, defined as

$$A^{ND}2\frac{𝒩^N𝒩^D}{𝒩^N+𝒩^D},$$

where $`𝒩^i`$ is the event rate measured during $`i=`$ Night, Day. Here data from elastic scattering (ES) of neutrinos on electrons at Super-K  has been augmented by data from SNO . The SNO data includes the day-night asymmetry in ES, in charged current (CC) $`\nu _ed`$ interaction as well as in the NC $`\nu _xd`$ interaction.

The SNO and Super-K ES data are consistent with each other, within error-bars. They indicate a small day-night effect. The SNO CC asymmetry is positive (about 2 $`\sigma `$ away from zero) while the NC asymmetry is barely negative (1.2 $`\sigma `$ away from zero).

It is well-known that oscillations among purely active neutrino flavours cannot generate a NC day-night asymmetry. In this paper, we examine whether it is possible to include interactions that generate even a small (whether positive or negative) day-night asymmetry in the NC sector. We find that this is impossible even with the inclusion of decay into the oscillation framework. In other words, if, for example, neutrino decay occurs through Majoron emission,

$$\nu _2\overline{\nu }_1\varphi ,$$

even the inclusion of both decay effects in the neutrino sector as well as of the antineutrinos that are produced in the decay does not yield rates that depend on the Earth-matter effects (that is, whether it is day or night at the time of detection). Hence the NC day-night asymmetry remains zero, as in the pure oscillation case. This is true for arbitrary number of (active) flavours, with some assumptions on the decay lifetimes.

We then analyse the rates in the NC, CC, and ES channels in a two-flavour framework, which is justified since the 1–3 mixing angle is small . The parameters are the two mixing parameters $`(\delta _{12}`$ and $`\omega )`$, and a decay parameter $`\alpha _2`$ that involves the life-time of the heavier mass eigenstate. It turns out that the contribution of the antineutrinos to the rates in any channel is small (depending on the life-time of $`\nu _2`$). Also, there is not much variation in the allowed oscillation parameter space for a large range of the decay parameter $`\alpha _2`$ in the CC and ES channels. However, the NC channel is very sensitive to $`\alpha _2`$ and thereby to the oscillation parameters. In fact, it turns out that the NC channel is the most sensitive measurement that constrains the decay parameter $`\alpha _2`$.

In summary, the day-night NC asymmetry is strictly zero in a framework involving both neutrino oscillations and decay (of active neutrino species). Of course, pure oscillations cannot also generate a day-night NC asymmetry. The contribution of anti-neutrinos to CC and ES day-night asymmetry is also small. Hence a clean signature for antineutrinos and hence for decay cannot come from the observation of such asymmetries: the best method remains to look for such signals in the rates rather than in the asymmetries: this will be signalled by two neutrons in coincidence with positron in the CC sector at SNO. However, the decay parameter itself is rather tightly constrained from the NC rates data.

In Section II we calculate the day-night asymmetries in a scenario that includes both neutrino oscillation and decay, for the NC, CC and ES channels. We show that the NC asymmetry is zero (leaving details to the Appendix), with some reasonable assumption on the decay life-times. Section III contains some numerical results on the allowed mixing parameter space from analysis of the CC, NC and ES rates at SNO  and the corresponding day-night asymmetries. Section IV contains a short discussion and conclusions. The Appendix includes some details of the day-night asymmetry calculations.

## 2 Neutrino and antineutrino survival probabilities with decay

We consider a framework where neutrino flavours $`\nu _\alpha `$ are mixtures of the mass eigenstates, $`\nu _i`$, and the heavier mass eigenstates $`\nu _i`$, $`i1`$, decay with life time $`\tau _i`$ ($`\nu _1`$ is assumed to be stable). We have,

$$|\nu _\alpha =\underset{i}{}U_{\alpha i}|\nu _i,$$
(1)
where the mixing matrix $`U_{\alpha i}`$ can be expressed in terms of the mixing angles. We ignore any phase dependence in the mixing matrix. The time evolution of the heavier mass eigenstates includes a decay term apart from the phase factor, $`\mathrm{exp}(iEt)`$:

$$|\nu _j(t)=\mathrm{exp}(iE_jt)\mathrm{exp}[t/(2\gamma _j\tau _{j_0})]|\nu _j(0),$$

where $`\tau _{j_0}`$ is the life-time in the neutrino rest frame and $`\gamma _j=E_j/m_j`$ is the usual Lorentz factor.

The (electron-type) neutrinos are produced in the core of the Sun. On their way out, they can interact with solar matter, and then reach the Earth. Depending on whether it is day or night at the detector, the neutrinos may further pass through Earth matter. The neutrino survival (or conversion) probabilities thus depend on the mixing matrix in the Sun, in the Earth, and the life-time of the heavier eigenstates. We will begin with the well-known expressions for the day-time probabilities. Details are included in the Appendix.

### 2.1 Day rates

The neutrino produced in the core at time $`t_0`$ propagates as a mass eignestate $`\nu _i`$ in the Sun until the time $`t_R`$ when it can undergo non-adiabatic transitions to another eigenstate, $`\nu _j`$. It then continues to propagate outwards, leaving the Sun at time $`t_1`$ and reaching the Earth at time $`t_2`$, where it is detected. The probability of a flavour $`\alpha `$ produced in the core being detected as flavour $`\beta `$ on Earth in the day-time can be written in a factorised form  as,

$$P_{\alpha \beta }^D=\underset{j}{}P_{\alpha j}^SP_{j\beta }f_j^D,$$
(2)
where the superscript $`S`$ indicates the dependence on solar parameters. The braces over $`P`$ indicate an averaging over the resonance region (or equivalently, an average over the production region ). We will always assume such an averaging and drop the braces in what follows. Note that the assumption of such a phase averaging restricts the neutrino mass-squared difference to be at least $`10^8`$ eV2.

The first term in the expression above represents the probability that the produced flavour $`\nu _\alpha `$ will traverse the solar matter and exit it as the mass eigenstate $`\nu _j`$. The second term is the (vacuum) probability of detection of this state as the flavour $`\nu _\beta `$ in the detector. The factor $`f_j^D`$ contains the decay information which also factors out, provided we make a single assumption about the decay life-times, as follows.

From Eq. 1, we see that, apart from the phase factor that arises on propagation, the probability amplitude includes decay terms that involve factors of the type,

$``$ $`\mathrm{exp}\left[{\displaystyle _I^F}{\displaystyle \frac{\mathrm{d}t}{2\tau _j}}\right],`$
$`=`$ $`\mathrm{exp}\left[{\displaystyle _I^F}{\displaystyle \frac{m_j}{2\tau _{0_j}}}{\displaystyle \frac{\mathrm{d}t}{E_j}}\right],`$
where $`(I,F)`$ refer to any of the time intervals $`(t_0,t_R)`$, $`(t_R,t_1)`$, or $`(t_1,t_2)`$. Now, if the decay rate is very large (life-time is small compared to the time taken to traverse the Sun-Earth distance), most of the neutrinos would have decayed by the time they reached the Earth. This is certainly ruled out by observations. Hence, the relevant life-times are at least of the order of the Sun-Earth distance. Due to this, only the time interval $`(t_1,t_2)`$ (when the neutrino propagates from the Sun to the Earth) is relevant for decay, and the other decay factors involving shorter time intervals can be neglected. With this single assumption, the decay factor $`f_j^D`$ factorises out of the probability expression in Eq. 2. Note that $`f^D=1`$ for the case of pure oscillations as well as the lightest (stable) neutrino. Detailed expressions for these terms are given in the Appendix. We note here that $`f_j^D`$ can be written in terms of the parameter $`\alpha _j`$ as $`f_j^D=\mathrm{exp}(\alpha _j/E_j)`$, where $`\alpha _j`$ involves the average Sun-Earth distance $`R`$ and a matter-independent term: $`\alpha _j=Rm_j/\tau _{j_0}`$  where $`m_j`$ are the (unknown) masses of the heavier eigenstates.

Similarly, the probabilities for an electron neutrino to be detected as an antineutrino, $`P_{\alpha \overline{\beta }}`$, at the detector can be computed and are also dependent on the decay parameters $`\alpha _j`$. An equation similar to Eq. 2 holds, with the matter-dependent mixing angle, $`\omega _m`$, being redefined with the opposite sign for the matter dependent term. For details, see the Appendix.

Finally, the event rate in the day time can be obtained by multiplying the incident flux, $`P_{e\beta }\mathrm{\Phi }_e`$ with the corresponding cross-section for detecting $`\nu _\beta `$ over the detector volume.

### 2.2 Night rates

It turns out that the conversion probabilities at night can also be written in a factorised form as the product of probabilities. While this is straightforward for the case of pure oscillations, in the case of decay, the factorisation again needs the assumption on the decay life-times used earlier. We have (see the Appendix for details),

$$P_{\alpha \beta }^N=\underset{j}{}P_{\alpha j}^SP_{j\beta }^Ef_j^N,$$
(3)
where $`P_{\alpha j}^S`$ is the same as for the day probability, but the second term includes an Earth matter effect that depends on the path length traversed through the Earth. Furthermore, within the approximations used, it turns out that $`f_j^N=f_j^D`$. Again, a similar equation (with $`AA`$ in the matter-dependent terms) holds for the case when $`\nu _\beta `$ is an antineutrino. It is important that the factorised form of the probabilities holds in both the day and night cases, for particle and antiparticle production; this is a crucial ingredient in showing that the NC asymmetry is zero.

Number conservation, or unitarity, implies that, for any $`j`$,

$`{\displaystyle \underset{\alpha }{}}P_{j\alpha }`$ $`=`$ $`1,`$
$`{\displaystyle \underset{\alpha }{}}P_{j\alpha }^E`$ $`=`$ $`1.`$
The first expression is valid in vacuum, the second in Earth. Similar expressions hold for the antineutrino case as well. In the neutral current channel, all flavours contribute equally. The expression that occurs in the NC night rates calculation below is,

$`{\displaystyle \underset{\alpha }{}}P_{e\alpha }^N`$ $`=`$ $`{\displaystyle \underset{j}{}}P_{ej}^Sf_j^N{\displaystyle \underset{\alpha }{}}P_{j\alpha }^E,`$ (4)
$`=`$ $`{\displaystyle \underset{j}{}}P_{ej}^Sf_j^N,`$
where we have used unitarity in the second line. This is independent of the Earth-matter effects and hence is the same as the day rate expression for the same sum that involves the vacuum probabilities, $`P_{j\alpha }`$:

$`{\displaystyle \underset{\alpha }{}}P_{e\alpha }^D`$ $`=`$ $`{\displaystyle \underset{j}{}}P_{ej}^Sf_j^D{\displaystyle \underset{\alpha }{}}P_{j\alpha },`$ (5)
$`=`$ $`{\displaystyle \underset{j}{}}P_{ej}^Sf_j^N,`$
where we have used $`f_j^D=f_j^N`$. This equality also holds for the sum $`_\alpha P_{e\overline{\alpha }}`$. This is a key ingredient in our result.

Trivially, the following sum rule is satisfied:

$$\underset{\alpha }{}\left[P_{e\alpha }+P_{e\overline{\alpha }}\right]=1.$$

### 2.3 Day-Night Asymmetries

We are now ready to write down expressions for the asymmetries observed at SNO. There are three asymmetries, corresponding to the CC, NC and ES events. We begin with the neutral current asymmetry. The cross-sections for the three cases are known, both for incident neutrinos and antineutrinos . It turns out that the cross-section for neutrinos (or antineutrinos) on deuterium is flavour independent in the NC channel. This is another crucial ingredient to our null result on the pure NC asymmetry including decay.

The total number of events in the day-time for the NC case is given by,

$`𝒩_{NC}^D`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}\left[𝒩_\alpha ^D+𝒩_{\overline{\alpha }}^D\right],`$
$`=`$ $`K_{NC}\mathrm{\Phi }_e{\displaystyle \underset{\alpha }{}}\left[P_{e\alpha }^D\sigma ^\alpha (E_\nu )+P_{e\overline{\alpha }}^D\sigma ^{\overline{\alpha }}(E_{\overline{\nu }})\right],`$
where $`P_{e\overline{\alpha }}`$ involves both the energy of the produced $`\nu _e`$ and that of the decay product, $`\overline{\nu }_\alpha `$. A similar expression holds for the night-time rate, with the superscript $`D`$ being replaced by $`N`$ everywhere. The overall constant factor $`K_{NC}`$ includes the detector volume (number of target atoms) as well as the time over which the measurement is made111Due to differences in lengths of day over the year, the day and night data are usually accumulated over different times; hence, the rates are to be normalised appropriately when calculating the asymmetry.. We now make use of the equality of the NC cross-sections, $`\sigma ^\alpha =\sigma `$ for all flavours, $`\alpha `$, and write,

$`𝒩_{NC}^D`$ $`=`$ $`K_{NC}\mathrm{\Phi }_e\left[\sigma {\displaystyle \underset{\alpha }{}}P_{e\alpha }^D+\overline{\sigma }{\displaystyle \underset{\alpha }{}}P_{e\overline{\alpha }}^D\right],`$
$`𝒩_{NC}^N`$ $`=`$ $`K_{NC}\mathrm{\Phi }_e\left[\sigma {\displaystyle \underset{\alpha }{}}P_{e\alpha }^N+\overline{\sigma }{\displaystyle \underset{\alpha }{}}P_{e\overline{\alpha }}^N\right],`$ (6)
$`=`$ $`K_{NC}\mathrm{\Phi }_e\left[\sigma {\displaystyle \underset{\alpha }{}}P_{e\alpha }^D+\overline{\sigma }{\displaystyle \underset{\alpha }{}}P_{e\overline{\alpha }}^D\right].`$
We see from Eqs. 4 and 5 that the day and night rates (when accumulated over the same period) are the same. Hence the NC day-night asymmetry,

$$A_{NC}^{ND}=0,$$

in a frame-work including both oscillation and decay of (active) neutrino flavours. This result does not hold if the decay occurs to a sterile flavour, since neutrinos are then “lost” and number conservation (the unitarity argument used above) does not hold. Finally, note that in deriving this result, no assumptions have been made on the details of the matter-dependence either in the Sun or in the Earth.

For completeness, we list the other asymmetries in the Appendix. The CC asymmetry involves only the electron neutrino interactions while the ES asymmetry includes the muon neutrino as well (and of course the antiparticles). We now move on to a numerical analysis of the data.

## 3 Numerical Calculations for rates and asymmetries

In this section, we compute the event rates in the CC, NC, and ES channels in an oscillation frame-work including decay. For simplicity, we work in a 2-flavour scenario so the parameters involved are the mixing parameters $`\delta _{12}`$ and $`\omega `$ and the decay parameter $`\alpha _2`$. In a 3-flavour frame-work, apart from additional mixing angles, $`\varphi =\theta _{13}`$ and $`\psi =\theta _{23}`$, and the corresponding mass-squared differences, there is another decay parameter, $`\alpha _3`$. However, there is very little sensitivity to the 1–3 mixing angle $`\varphi `$  which is known to be small, and no dependence at all of solar neutrino rates on the 2–3 angle $`\psi `$. Furthermore, in the limit of small $`U_{e3}`$ (and hence small $`\varphi `$), which is driven by the CHOOZ result , it was shown in Ref.  that the $`\alpha _3`$ dependence drops out of the expressions for the day-time rates so the decay is still governed only by the second-heaviest eigenstate through $`\alpha _2`$. In short, the 2-flavour analysis is sufficiently sensitive to most of the features of the model; we shall therefore present numerical results within such a frame-work.

The detailed expressions for the rates in a two-flavour frame-work are listed in the Appendix. We shall use these formulae to explicitly compute the rates in the CC, NC, and ES channels. We use a chi-squared minimisation program to best-fit the free parameters and obtain the allowed regions in the parameter-space. We also compute, for completeness, the day-night asymmetry in these channels corresponding to the best-fit values of the parameters.

### 3.1 Allowed parameter region from rates analysis

As stated before, we use the cross-sections given by Ref.  for the neutrino-deuteron CC and NC cross-sections and that given by Ref.  for the ES cross-sections. We assume a nearly massless majoron so that the antineutrino produced by decay has half the energy of the initial neutrino.

We use the fiducial volume of 0.77 ktons of heavy water and the zenith angle dependence for 128.5 “days” of day data and 177.9 “days” of night data as given by SNO . We use the normalisation of NC events (with kinetic energy of detected electron $`T_e>5`$ MeV) as indicated by SNO . We also use the resolution function for the CC, ES and NC case as given by SNO. For the theoretical predictions, we use the standard solar model of Bahcall, Pinsonneault and Basu . We use the Boron neutrino flux and ignore the hep contribution. Note that there are about 20% errors on the Boron flux of neutrinos.

We use a simple one-slab model of the Earth with an average density of 5 gm/cc. Due to the latitude location of SNO, it detects neutrinos that have passed through the (outer) core of the Earth roughly 8.7% of the time. (SNO does not see any solar neutrinos that pass through the inner core of the Earth). By assuming a one-slab model that only incorporates the mantle, we have effectively neglected the contribution from this fraction of events. We use the 2-flavour expressions for the probabilities as given in the Appendix; the parameters to be determined by comparison with the SNO data are the mass squared difference $`\delta _{12}`$, the mixing angle $`\omega `$ and the decay parameter $`\alpha _2`$.

Before we begin the actual numerical calculation, we note that the NC event rate is amenable to some approximations. Since the antiparticle cross-section is small (the cross-section depends non-linearly on the energy and the average energy of the antiparticles is half that of the particles), the contribution from the antineutrino sector can be neglected. (This is also true in the CC case). We have,

$`𝒩`$ $``$ $`\mathrm{\Phi }\left[(P_{e1}^S+P_{e2}^S\mathrm{exp}(\alpha _2/E_2))\sigma +P_{e2}^S(1\mathrm{exp}(\alpha _2/E_2))\overline{\sigma }\right],`$ (7)
$`=`$ $`\mathrm{\Phi }\left[\sigma P_{e2}^S(1\mathrm{exp}(\alpha _2/E_2))\left\{\sigma \overline{\sigma }\right\}\right],`$
$``$ $`\mathrm{\Phi }\sigma \left[1P_{e2}^S(1\mathrm{exp}(\alpha _2/E_2))\right].`$ (8)
Note that the NC events are independent of mixing in the absence of decay. Hence, when $`\alpha _2=0`$, the second term is zero and the NC rate is independent of the mixing parameters. In fact, it is proportional to $`\mathrm{\Phi }\sigma `$, which is the total rate in the absence of neutrino oscillation/decay. As decay is turned on (so that $`\alpha _20`$), the NC rate begins to depend on the mixing parameters through $`P_{e2}^S`$ ($`=1P_{e1}^S`$ for two-flavours). Since the NC rates are well within that predicted by the standard solar model, to counteract the effect of a non-zero $`\alpha _2`$, the parameter space in mixing gets squeezed in order to minimise $`P_{e2}^S`$ (or to maximize $`P_{e1}^S`$) so that the second term remains small. Since it is not possible to squeeze $`P_{e1}^S`$ beyond a point (as it is bounded by unity), larger values of $`\alpha _2`$ get ruled out by the SNO result that the NC event rates are within 1 $`\sigma `$ of the standard solar model . Hence it turns out that the NC channel is the best suited to constrain the decay parameter. Finally, such a squeezed solution will be a bad fit for the CC and ES rates, since here the rates are depleted compared to expectations. Hence a combined analysis of the CC, NC, and ES data will severely constrain $`\alpha _2`$.

With this heuristic discussion, we proceed to calculate the rates from different channels at SNO. We have performed a combined chi-squared analysis of the total rates at SNO. The CC and NC rates are strongly anticorrelated while the CC/NC and ES rates are mildly anticorrelated; however, we have not included this correlation in our analysis.

While SNO CC data favours a LOW solution (low mass squared difference, large mixing angle) we know that when combined with data from other experiments such as GALLEX , this solution is ruled out . Hence we concentrate on solutions, including decay, around the “pure oscillation” LMA and SMA solutions. When $`\alpha _2=0`$ the solutions are determined by the CC and ES rates alone. In fact this turns out to have the lowest $`\chi ^2`$ value when $`\alpha _2`$ is turned on. The best-fit values to the parameters for different values of $`\alpha _2`$ are shown in Table 1. It is seen that $`\alpha _2=5`$ MeV is disallowed at 2 $`\sigma `$.

The 1 $`\sigma `$ and 2 $`\sigma `$ allowed parameter space from fitting the combined CC, NC and ES rates in $`(\delta _{12},\mathrm{tan}^2\omega )`$ for $`\alpha _2=0,1,5`$ MeV are shown in Fig. 1. There is very little difference between the allowed regions for $`\alpha _2=0,1`$ MeV. However, there are no solutions at 1 $`\sigma `$ for $`\alpha _2=5`$ MeV. This is a much stronger result than the limit $`\alpha _2<18`$ MeV obtained by considering only the CC and ES data . The 2 $`\sigma `$ allowed regions are shown on the right hand side of Fig. 1. We see that the “usual” LMA region for $`\alpha _2=5`$ MeV is disallowed by the pinching effect on the NC rates described earlier.

This result is easy to understand if we examine the allowed regions from the individual channels. We find that the 1 $`\sigma `$ allowed parameter space in mixing for a given $`\alpha _2`$ from CC data and ES data overlap almost completely, while the allowed region from NC data is very different. The allowed regions coming from individual analyses of NC and CC rates are shown in Fig. 2.

We know that the NC rate is independent of mixing parameters for $`\alpha _2=0`$. Even when $`\alpha _2=1`$ MeV, the entire parameter space is allowed at 1 $`\sigma `$ level as can be seen from the left hand side of Fig. 2. For $`\alpha _2=5`$ MeV, there is no overlap region between the allowed regions from CC/ES and NC data; hence there are no 1 $`\sigma `$ allowed regions for such large values of $`\alpha _2`$. In the LMA region, for instance, the allowed region from the NC channel shifts towards smaller $`\omega `$ in an attempt to maximise $`P_{e1}^S`$, as discussed earlier. This results in eliminating the region in parameter space that is preferred by the CC rate.

Note that we have fitted the total rates and not the energy-dependent spectrum. Inclusion of larger amounts of decay will distort the spectrum significantly and so it seems as if analysis of the spectrum may limit the decay parameters better than the rates. However, as $`\alpha _2`$ becomes large, there is no allowed parameter space and hence no observable spectral distortion; a rates analysis is therefore unlikely to further constrain this result. This is borne out by the consistency of our results with the detailed spectral analysis, including decay, found in Ref. , where a limit $`\alpha _2<5.7`$ MeV is obtained.

### 3.2 The day-night asymmetries

We calculate the day-night asymmetries for CC and ES channels using the best-fit values of the parameters for different values of $`\alpha _2`$ as given in Table 1. We find that the CC asymmetry is roughly 0.2 and positive for allowed $`\alpha _2`$ in the LMA region of mixing parameters, while the ES asymmetry is slightly smaller and also positive. The asymmetries are much smaller in the SMA case. When $`\alpha _2`$ exceeds 4 MeV or so, the asymmetries are nearly zero.

The corresponding day-night asymmetries are shown as a function of the observed kinetic energy of the electron in Fig. 3. There is a negligible contribution from the antiparticles, which contribute only at low energies, $`T<5`$ MeV. The energy-averaged values of the SNO CC and ES asymmetries are shown alongside for comparison. Most of the events correspond to lower energies; hence the average asymmetry is closer to the value at the lower energy end. There is no significant constraint on the decay parameter coming from the day-night asymmetries; this of course may change with the availability of more data from SNO.

## 4 Discussion and Conclusions

We have shown that the neutral current day-night asymmetry is zero, in a frame-work of neutrino mixing that includes oscillations plus decay, with a reasonable assumption on the size of the neutrino life-time. We discussed the effect of decay on the neutral-current, charged-current and elastic scattering rates in a two-flavour frame-work. In all cases, it turns out that the contribution from antiparticles is small in the region of allowed decay parameter space. This is because of small cross-sections owing to small (roughly half) energies of antineutrinos compared to those of neutrinos. Hence, even if there are large decay rates, very few events will be seen. We find that the neutral-current rates constrain the decay parameter quite stringently. On comparing with SNO data for the total rates in the CC, NC, and ES channels, we get $`\alpha _2<5`$ MeV, which translates to a bound on the life-time of the heavier mass eigenstate of $`\tau _{2_0}>10^4`$ s for eV mass neutrinos (equivalently, this is the bound on $`\tau _{2_0}/m`$ for an unknown neutrino mass $`m`$ in eV). This is very close to the bound that already been obtained in a detailed analysis of solar neutrino data in Ref. . This is much stronger than the bound $`\alpha _2<18`$ MeV obtained in Ref.  using only the charged-current SNO data; in fact, it is the limiting bound that can be obtained, given the Sun-Earth distance . It is interesting that this limit is obtained only on including the neutral-current data.

Here we have mainly highlighted the dual role of the NC SNO data: it is most constraining of the decay parameter as far as total event rates are concerned; however, the NC day-night asymmetry is completely oblivious to the parameters of either mixing or decay. An in-depth analysis of available solar neutrino data, in a 2- and 3-flavour frame-work including decay, has been done in Ref. , where detailed energy correlations as well as rates from other experiments such as the Chlorine and Gallium experiments , have been analysed along with the energy-dependent spectra and day-night asymmetries from Super-K and SNO . Our results are consistent with this detailed analysis.

Given the allowed range of $`\alpha _2`$ coming from the rates analyses, we see that the corresponding day-night asymmetries are not very sensitive to the mixing parameters or indeed to the decay parameter. Hence, the best way to put limits on, or observe neutrino decay, will be by a direct measurement of the rates and not the asymmetries. A distinctive signature of an “electron” detected in coincidence with two neutrons will signal the capture of an antineutrino. This occurs through the direct observation via the CC interaction in heavy water at SNO:

$$\overline{\nu }_ede^+nn.$$
(9)
Two neutrons must be seen in coincidence with the positron. If the neutrons are not observed, then this process is indistinguishable from the dominant one of

$$\nu _ede^{}pp,$$

since water Cerenkov detectors cannot distinguish electrons from positrons. Furthermore, both the $`\nu _e`$ and the $`\overline{\nu }_e`$ interactions have precisely the same angular dependence (which is $`2:1`$ backward-forward asymmetric) so that angular distribution cannot be used to separate the two signals. Finally, even if the antineutrino energy is half that of the decaying neutrino, SNO will still be able to detect them through the neutron coincidence measurement. It is expected that there will soon be a substantial improvement in the neutron detection efficiency (with the inclusion of salt in the SNO detector); it may then be possible to detect these channels with good efficiency.

Note here that while the interaction of (electron) antineutrinos on protons in water is very large (about 100 times the cross-section of neutrinos in water at solar neutrino energies), Super-Kamiokande will not see these antineutrinos since their angular distribution is roughly isotropic and Super-K cuts out non-forward events.

Finally, a remark about the NC day-night asymmetry measured by SNO. After all, this was the motivating factor for this calculation. Currently, this asymmetry is barely negative (about 1.2$`\sigma `$ away from zero). If, with improved statistics, this asymmetry persists, it will be difficult to find a source for the effect. It is possible then that there are sterile neutrinos involved in the oscillation frame-work; however, the extent of conversion to steriles is limited by the total NC event rate (which is already consistent with the BP2000 standard solar model prediction) which admittedly has large theoretical errors (of about 20%). Another possible source for a non-zero NC asymmetry could be that the NC cross-sections for $`\nu _e`$ and $`\nu _\mu `$, $`\nu _\tau `$ scattering on deuteron are not the same.

So, with the improved detection of neutrons, and the accumulation of more data, SNO will be in a good position to set strong bounds on neutrino decay and oscillations into sterile flavours.

Acknowledgements: I thank M.V.N. Murthy and G. Rajasekaran for many discussions and a critical reading of the manuscript.

## Appendix

We derive the neutrino conversion probabilities in a frame-work that includes both oscillation and decay. A flavour $`|\nu _\alpha `$ that is produced in the core of the sun (actually, an electron type neutrino), can be expressed in terms of the mass eigenstates $`|\nu _i`$, as

$$|\nu _\alpha =\underset{i}{}U_{\alpha i}^S|\nu _i^S,$$

where $`U_{\alpha i}`$ is the usual mixing matrix that can be expressed in terms of a mixing angle $`\omega `$ in the case of two-flavours, or three angles ($`\omega `$, $`\varphi `$, $`\psi `$) for the three-flavour case (we ignore any CP violating phases). The superscript $`S`$ indicates that these are the mass eigenstates in the Sun. The effect of decay is to modify the phase factor that determines the time-evolution of the mass eigenstates so that the time evolution is given by,

$$|\nu _i_t=\mathrm{exp}[iE_it]\mathrm{exp}[t/(2\tau _i)]|\nu _i_0.$$

The lifetime can be expressed in terms of that in the rest frame: $`\tau _i=\gamma _i\tau _{0_i}`$, $`\gamma _i=E_i/m_i`$.

The neutrino, produced at a time $`t_0`$ in the core of the Sun, propagates adiabatically upto the resonance point $`t_R`$, where it can undergo non-adiabatic conversion to a different mass eigenstate, and then continues to propagate adiabatically to the edge of the Sun which it reaches at $`t_1`$. Then, depending on whether it is night or day, the neutrino propagates upto the Earth and gets detected at time $`t_2`$, or traverses a further distance $`L`$ through the Earth and gets detected at time $`t_3`$.

Hence, the probability of a flavour $`\alpha `$ being detected as a flavour $`\beta `$ on Earth in the day-time is given by

$$P_{\alpha \beta }^D=\underset{ij}{}|U_{\beta j}|^2|U_{\alpha i}^S|^2|M_{ji}^S|^2f_j^D,$$

where $`U`$ is the vacuum mixing matrix, and we have averaged over the resonance region . The decay factor $`f_j^D`$ will be defined below. Here $`M_{ji}^S`$ represents the non-adiabatic transition probability amplitude from $`ij`$ so that $`|M_{ji}^S|^2`$ is the usual Landau Zener transition probability  in the Sun.

As has been pointed out earlier , it is possible to express the probability as a factor that is solar-matter dependent, times a solar-matter independent factor:

$$P_{\alpha \beta }^D=\underset{j}{}P_{\alpha j}^SP_{j\beta }f_j^D,$$
(1)
where

$$P_{\alpha j}^S=\underset{i}{}|U_{\alpha i}^S|^2|M_{ji}^S|^2,$$

and

$$P_{j\beta }=|U_{\beta j}|^2.$$

As explained in the main part of the text, the decay factor $`f_j^D`$ in the expression for the probability amplitude involves a product of terms:

$``$ $`\mathrm{exp}\left[{\displaystyle _I^F}{\displaystyle \frac{\mathrm{d}t}{2\tau _i}}\right],`$
$`=`$ $`\mathrm{exp}\left[{\displaystyle _I^F}{\displaystyle \frac{m_i}{2\tau _{0_i}}}{\displaystyle \frac{\mathrm{d}t}{E_i}}\right],`$
with one such term coming from every propagation from $`t_I`$ to $`t_F`$ of a neutrino $`\nu _i`$. In general, we can assume the modification of the energy eigenvalue due to the matter-dependent effects is small, so that the integral above reduces to $`(t_Ft_I)[m_i/(2\tau _{0_i}E_i)]`$ . If the time interval of propagation is short, then the numerator is small compared to the life-time unless the latter is itself very small. In the latter case, we would expect most of the produced neutrinos to have decayed before reaching the Earth, but this is clearly not true. We therefore assume that the life-time of the heavier mass eigenstate is at least the Sun-Earth distance (divided by the speed of light). Then the decay exponentials are all small and can be dropped, except the terms containing the time interval $`(t_Ft_I)=(t_2t_1)`$ which is the Sun-Earth separation. The relevant integral contains the vacuum energy eigenvalue and hence is anyway matter-independent. The decay term that appears in the expressions for the probability is then the square of the above, with no cross terms, for the energy eigenvalue $`E_j`$ corresponding to $`\nu _j`$ that propagates in vacuum from the Sun to the Earth,

$$f_j^D=\mathrm{exp}(\alpha _j/E_j),$$

where we have used $`\alpha _j=Rm_j/\tau _{0_j}`$  where $`R`$ is the Sun-Earth distance.

A similar formula holds at night-time, except that the vacuum probability $`P_{j\beta }`$ is replaced by an earth-matter dependent one:

$$P_{\alpha \beta }^N=\underset{j}{}P_{\alpha j}^SP_{j\beta }^Ef_j^N,$$
(2)
with the solar probability the same as before, and with

$$P_{j\beta }^E=\underset{kk^{}}{}U_{\beta k}^EU_{\beta k^{}}^EM_{kj}^EM_{k^{}j}^E\mathrm{exp}\left[i_{t_2}^{t_3}(E_k^EE_k^{}^E)dt\right].$$

Here $`M_{kj}^E=\nu _k^E|\nu _j`$ is the transition amplitude from $`jk`$, where $`\nu _k^E`$ is the mass eigenstate just inside the Earth. We have assumed only one non-adiabatic transition occurs, when the neutrino encounters the discontinuity in density at the Earth’s surface. However, it is straightforward to generalise this to an arbitrary number of discontinuous jumps in Earth’s density, using the so-called slab model of the Earth. Details may be found in Ref. . The decay factor is the same as in the day-time case since the matter effects (whether in Earth or in the Sun) due to decay are small due to the small distance of propagation in these. So the vacuum contribution is the only one that remains here as well. So we have

$$f_j^N=f_j^D.$$
(3)
Finally, the event rate can be obtained by multiplying the incident flux, $`P_{e\beta }\mathrm{\Phi }_e`$ with the corresponding cross-section for detecting $`\nu _\beta `$ over the detector volume.

Similar formulae hold for antineutrinos produced during decay, with all matter-dependent terms being redefined with $`AA`$, where $`A`$ is the matter dependent term in the Sun or Earth. It is important that the factorised form of the probabilities holds in both the day and night cases, for particle and antiparticle production; this is a crucial ingredient in showing that the NC asymmetry is zero.

Number conservation, or unitarity, implies that, $`j`$,

$`{\displaystyle \underset{\alpha }{}}P_{j\alpha }`$ $`=`$ $`1,`$
$`{\displaystyle \underset{\alpha }{}}P_{j\alpha }^E`$ $`=`$ $`1.`$
The first expression is valid in vacuum, the second in Earth. Note that the second equation is valid even when Eq. 2 for the night-time probability is replaced by a more complicated one using a many-slab model of the Earth.

Similar expressions hold for the antineutrino case as well. In the neutral current channel, all flavours contribute equally. The expression that occurs in the NC night rates calculation is,

$`{\displaystyle \underset{\alpha }{}}P_{e\alpha }^N`$ $`=`$ $`{\displaystyle \underset{j}{}}P_{ej}^Sf_j^N{\displaystyle \underset{\alpha }{}}P_{j\alpha }^E,`$ (4)
$`=`$ $`{\displaystyle \underset{j}{}}P_{ej}^Sf_j^N,`$
where we have used unitarity in the second line. This is independent of the Earth-matter effects and hence is the same as the day rate expression for the same sum that involves the vacuum probabilities, $`P_{j\alpha }`$:

$`{\displaystyle \underset{\alpha }{}}P_{e\alpha }^D`$ $`=`$ $`{\displaystyle \underset{j}{}}P_{ej}^Sf_j^D{\displaystyle \underset{\alpha }{}}P_{j\alpha },`$ (5)
$`=`$ $`{\displaystyle \underset{j}{}}P_{ej}^Sf_j^N,`$
where we have used $`f_j^D=f_j^N`$. This equality also holds for the sum $`_\alpha P_{e\overline{\alpha }}`$. From this it is clear that the NC rates are the same during the day or night for arbitrary numbers of active neutrino species; hence the NC day-night asymmetry is zero as explained in the main part of the text. If the decay is into a sterile flavour, of course, this is no longer true.

We now specialise to the two-flavour case where we list the detailed expressions for the probabilities needed for the computation of the rates and hence the asymmetries in the different channels at SNO. Then with the choice $`m_2>m_1`$, $`\nu _1`$ is stable so that $`\tau _2`$ is the only relevant decay parameter. We have,

$$f_2^D=f_2^N=\mathrm{exp}[\alpha _2/E_2].$$

We now list the various probabilities. Since the Earth’s density is much smaller than that in the solar core, there is only a small modification of the energy due to Earth-matter effects. Then we can define

$$X=(t_3t_2)\delta _{12}^E/(2E),$$
(6)
where $`\delta _{12}^E`$ is the mass-squared difference in (Earth) matter:

$$\delta _{12}^E=\delta _{12}\mathrm{cos}(2(\omega \omega _e))A_E\mathrm{cos}2\omega _e,$$

where the Earth matter dependent term is, $`A_E=2\sqrt{2}G_FN_e^EE`$ with $`N_e^E`$ the electron number density in the Earth and $`(t_3t_2)=2R_E\mathrm{cos}\theta _z`$ and the matter mixing angle is $`\omega _e`$. Here $`2R_E`$ is the diameter of the Earth ($`=5.1R_E(m)\text{MeV/eV}^2`$) and $`\theta _z`$ is the zenith angle of propagation of the neutrino with respect to the instantaneous direction of the Sun.

The Day-time probabilities: In the two-flavour case, $`P_{e2}=1P_{e1}`$; furthermore, we have,

$`P_{ee}^D`$ $`=`$ $`c^2P_{e1}^S+s^2P_{e2}^S\mathrm{exp}(\alpha _2/E_2),`$
$`P_{e\mu }^D`$ $`=`$ $`s^2P_{e1}^S+c^2P_{e2}^S\mathrm{exp}(\alpha _2/E_2),`$
$`P_{e\overline{e}}^D`$ $`=`$ $`c^2P_{e2}^S[1\mathrm{exp}(\alpha _2/E_2)],`$
$`P_{e\overline{\mu }}^D`$ $`=`$ $`s^2P_{e2}^S[1\mathrm{exp}(\alpha _2/E_2)].`$ (7)
Here $`c`$ and $`s`$ refer to $`\mathrm{cos}\omega `$ and $`\mathrm{sin}\omega `$ respectively. Defining the Landau-Zener probability $`P_{LZ}`$  for transition between the $`\nu _1`$ and $`\nu _2`$ eigenstates (note that the third mass eigenstate, even if included, does not participate in non-adiabatic transitions in the Sun, provided the relevant mass squared difference is that which solves the so-called atmospheric neutrino problem), we have

$$P_{e1}^S=s_m^2P_{LZ}+c_m^2(1P_{LZ}),$$
(8)
where $`s_m`$ and $`c_m`$ refer to $`\mathrm{sin}\omega _m`$ and $`\mathrm{cos}\omega _m`$ in the solar matter, given by

$$\mathrm{cos}2\omega _m=\frac{\delta _{12}\mathrm{cos}2\omega A}{\sqrt{(\delta _{12}^2\mathrm{sin}^22\omega +(\delta _{12}\mathrm{cos}2\omega A)^2)}}.$$
(9)
The above formula holds for antineutrinos with $`AA`$ and for neutrinos and antineutrinos in Earth matter as well, with $`AA_E`$ so that $`\omega _m(\pm A)\omega _e(\pm A_E)`$.

Note that

$`P_{ee}^D+P_{e\mu }^D`$ $`=`$ $`P_{e1}^S+P_{e2}^S\mathrm{exp}(\alpha _2/E_2),`$
$`P_{e\overline{e}}^D+P_{e\overline{\mu }}^D`$ $`=`$ $`P_{e2}^S(1\mathrm{exp}(\alpha _2/E_2)),`$ (10)
so that the total conversion probability into neutrinos and antineutrinos is conserved. When $`\alpha _2=0`$ we recover the usual formulae for pure oscillation case.

Night-time probabilities: Here we need to compute, in addition, the Earth-matter effects. We have

$`P_{1e}^E`$ $`=`$ $`12c_e^2s_e^2s^2\mathrm{cos}^22\omega _e+2csc_es_e\mathrm{cos}2\omega _e`$
$`+2c_es_e(c_es_e\mathrm{cos}2\omega cs\mathrm{cos}2\omega _e)\mathrm{cos}X,`$
$`P_{2e}^E`$ $`=`$ $`12c_e^2s_e^2c^2\mathrm{cos}^22\omega _e2csc_es_e\mathrm{cos}2\omega _e`$
$`2c_es_e(c_es_e\mathrm{cos}2\omega cs\mathrm{cos}2\omega _e)\mathrm{cos}X,`$
$`P_{1\mu }^E`$ $`=`$ $`P_{2e}^E,`$
$`P_{2\mu }^E`$ $`=`$ $`P_{1e}^E,`$ (11)
where $`X`$ has been defined above and $`\omega _e`$ is the Earth-matter dependent mixing angle. Notice that, along with the usual number conservation, $`P_{1e}^E+P_{1\mu }^E=1`$, we have also $`P_{1e}^E+P_{2e}^E=1`$ . Substituting in the transition probabilities, $`P_{\alpha \beta }^N`$, this gives us,

$`P_{ee}^N+P_{e\mu }^N`$ $`=`$ $`P_{e1}^S(P_{1e}^E+P_{1\mu }^E)+P_{e2}^S(P_{2e}^E+P_{2\mu }^E)\mathrm{exp}(\alpha _2/E_2)`$ (12)
$`=`$ $`P_{e1}^S+(1P_{e1}^S)\mathrm{exp}(\alpha _2/E_2).`$
Notice that the sum is independent of Earth matter effects and is exactly equal to the “Day” formula; this was already shown in the general case of arbitrary number of flavours.

For the detection of antiparticles at night, we have similar expressions to those for the neutrinos, with $`\omega _e(A_E)\omega _e(A_E)`$:

$`P_{e\overline{e}}^N+P_{e\overline{\mu }}^N`$ $`=`$ $`P_{e\overline{1}}^S(P_{1\overline{e}}^E+P_{1\overline{\mu }}^E),`$
$`=`$ $`(1P_{e1}^S)(1\mathrm{exp}(\alpha _2/E_2)).`$
Again, the total probability is conserved and is independent of the Earth-matter effects.

The event rates: With these probabilities, we now define the rates for CC, NC and ES events at SNO. We begin with the neutral current channel. The event rate at SNO in the NC channel is given by

$$𝒩^{NC}=𝒩^e+𝒩^\mu +𝒩^{\overline{e}}+𝒩^{\overline{\mu }}.$$

Hence,

$`𝒩_{Day}^{NC}`$ $`=`$ $`\mathrm{\Phi }\left[P_{ee}^D\sigma ^e+P_{e\mu }^D\sigma ^\mu +P_{e\overline{e}}^D\sigma ^{\overline{e}}+P_{e\overline{\mu }}^D\sigma ^{\overline{\mu }}\right],`$ (13)
$`=`$ $`\mathrm{\Phi }\left[\sigma P_{e2}^S(1\mathrm{exp}(\alpha _2/E_2))\left(\sigma \overline{\sigma }\right)\right];`$
$`𝒩_{Night}^{NC}`$ $`=`$ $`\mathrm{\Phi }\left[P_{ee}^N\sigma ^e+P_{e\mu }^N\sigma ^\mu +P_{e\overline{e}}^N\sigma ^{\overline{e}}+P_{e\overline{\mu }}^N\sigma ^{\overline{\mu }}\right],`$ (14)
$`=`$ $`\mathrm{\Phi }\left[\sigma P_{e2}^S(1\mathrm{exp}(\alpha _2/E_2))\left(\sigma \overline{\sigma }\right)\right],`$
$`=`$ $`𝒩_{Day}^{NC}.`$
These expressions leads us to the general result that we had shown earlier: that the NC day-night asymmetry is zero. Here $`\sigma `$ ($`\overline{\sigma }`$) is the flavour-independent neutral current cross-section of neutrinos (antineutrinos) on deuteron. For calculating the cross-section and hence the rates, we make the simple assumption that in the decay, the antiparticle carries away about half the energy of the parent neutrino. Hence, although the NC cross-sections for neutrinos and antineutrinos are roughly the same at the same energy, the contribution from antineutrinos is very small since the cross-section rises exponentially with energy. Indeed, for most values of $`\alpha _2`$ the bulk of the contribution to the rates is still from neutrinos.

Similar expressions can be written for the CC and ES rates. In the former, only the $`\nu _e`$ (or $`\overline{\nu }_e`$) contribute. We have,

$`𝒩_{Day}^{CC}`$ $`=`$ $`\mathrm{\Phi }\left[P_{ee}^D\sigma _{CC}^e+P_{e\overline{e}}^D\sigma _{CC}^{\overline{e}}\right],`$
$`𝒩_{Night}^{CC}`$ $`=`$ $`\mathrm{\Phi }\left[P_{ee}^N\sigma _{CC}^e+P_{e\overline{e}}^N\sigma _{CC}^{\overline{e}}\right],`$ (15)
where the probabilities involved have been listed earlier and the cross-sections are taken from Ref. . There is no simplification of the expressions although the anti-neutrino contribution is small here as well.

The ES rates involves all particles and antiparticles, but with flavour-dependent cross-sections as given in Ref. . We have,

$`𝒩_{Day}^{ES}`$ $`=`$ $`\mathrm{\Phi }\left[P_{ee}^D\sigma _{ES}^e+P_{e\mu }^D\sigma _{ES}^\mu +P_{e\overline{e}}^D\sigma _{ES}^{\overline{e}}+P_{e\overline{\mu }}^D\sigma _{ES}^{\overline{\mu }}\right],`$
$`𝒩_{Night}^{ES}`$ $`=`$ $`\mathrm{\Phi }\left[P_{ee}^N\sigma _{ES}^e+P_{e\mu }^N\sigma _{ES}^\mu +P_{e\overline{e}}^N\sigma _{ES}^{\overline{e}}+P_{e\overline{\mu }}^N\sigma _{ES}^{\overline{\mu }}\right].`$ (16)
The cross-section for elastic scattering of the $`\nu _e`$ is about 6 times larger than than of $`\nu _\mu `$ (and even larger than the antineutrino cross-sections), although the values are energy dependent.

The day-night asymmetries: The CC and ES asymmetries are not zero; the former depends only on $`P_{ee}`$ and $`P_{e\overline{e}}`$. These are different for the day and night-time. Hence, the CC asymmetry is non-zero by an amount proportional to

$$𝒩_{Night}^{CC}𝒩_{Day}^{CC}=\mathrm{\Phi }\underset{i}{}\left[(P_{ei}^NP_{ei}^D)\sigma _{CC}^i\right],$$
(17)
where $`i=e,\overline{e}`$. It is straightforward to similarly write an expression for the ES asymmetry, which has both charged-current and neutral-current contributions:

$$𝒩_{Night}^{ES}𝒩_{Day}^{ES}=\mathrm{\Phi }\underset{i}{}\left[(P_{ei}^NP_{ei}^D)\sigma _{ES}^i\right],$$
(18)
where the summation includes contributions from electrons and muons (and their antiparticles as well).