text:
# A numerical study of the Schrödinger-Newton equation 2: the time-dependent problem

## 1 Introduction

Ths is the second in a series of papers presenting a numerical study of the Schrödinger-Newton (or SN) equations. In the first () we reviewed known results on the SN equations and then analysed the linear stability of the spherically-symmetric stationary states. We found that the ground state is linearly stable, in that all the eigenvalues for linear perturbations were purely imaginary, while all the higher states are unstable. The $`(n+1)`$-th state, or equivalently the $`n`$-th excited state, has $`n`$ quadruples $`(\lambda ,\lambda ,\overline{\lambda },\overline{\lambda })`$ of complex eigenvalues with nonzero real part. We now want to go further in two directions. On the one hand we shall consider the full nonlinear evolution of the SN equations and on the other we shall move away from spherical symmetry, allowing for two space variables.

In the evolution of spherically-symmetric data we shall find that the picture from linear theory is confirmed: the ground state is stable but slight perturbations of the higher states decay, with probability density escaping off the grid to leave a ‘nugget’ consisting of the ground state rescaled to have total probability less than one and situated at the origin. For more general spherically-symmetric data the picture is confirmed: solutions typically disperse leaving a nugget of rescaled ground state. (By a rescaled ground state we mean the following: if $`\psi _0(r)`$ is the spatial part of the wave function for the ground state, normalised to have total probability equal to one and with conserved energy $`_0`$ in the terminology of , then $`p^2\psi _0(pr)`$ with $`p`$ a positive constant, is a stationary solution with conserved energy $`p^3_0`$ and total probability $`p`$.)

We can introduce a new space variable in two ways. We may consider axially-symmetric solutions in 3-dimensions. In this case there are new 3-dimensional stationary states, somewhat analogous to the axially-symmetric solutions of the hydrogen atom, and with nonzero expectation for the angular momentum. These turn out to be unstable under evolution, as one would expect. For the solution most analogous to a pure dipole, the evolution clearly shows the two ‘lumps’ of probability falling from rest into each other and evolving towards the ground state, with some scatter.

Alternatively we may consider translation-invariant solutions in 3-dimensions or equivalently the SN equations in $`2+1`$-dimensions. (Of course the solutions are not normalisable as $`3+1`$-dimensional solutions but this is still an interesting problem.) In this case we can find rigidly-rotating stationary states like two lumps of probability rotating around each other. When these are evolved with the time-dependent SN equations they prove to be unstable, so that it is possible to arrange for two lumps of probability to be in orbit around each other at least for a while before they merge.

In summary, the picture that we find is of a system with dispersive tendencies because of the Schrödinger equation and attractive or concentrating tendencies from the gravitational attraction. There appear to be infinitely many stationary states, all of them unstable except for the ground state. General data evolve to leave some residual probability in a rescaled ground state with the rest of the probability dispersing. Lumps of probability in the initial data can attract each other and even orbit each other but eventually the dispersive tendencies win.

The plan of the paper is as follows. We shall end this Introduction with an analytic result bounding the residual probability left on the grid at late times. In Section 2 we introduce the numerical method which we shall use to evolve the spherically-symmetric SN equations. The results of this evolution are presented in Section 3. In Section 4 we present the results of numerically evolving 3-dimensional axisymmetric solutions and in Section 5 we present the results of numerically evolving the 2-dimensional SN equations.

For the analytic calculation then, suppose that the picture presented above holds in general so that arbitrary initial data evolve to give a scattering solution which disperses to infinity and leaves a rescaled ground state. Then we can obtain a bound on the probability remaining in the ground state provided the initial energy is negative. For suppose the initial value of the conserved energy is $`_I`$ then since this is conserved, we can calculate at late times to find

$$_I=_S+p^3_0$$

where $`_S`$ is the energy in the scattering solution, which we suppose to be positive, $`_0`$ is the (negative) energy of the ground state and $`p`$ is the probability left in the ground state. (The ground state is strongly peaked at the origin so that the cross-term in $`_I`$ is zero.) Taking account of all the signs, if $`_I<0`$ we can rearrange this to read

$$p^3>\frac{|_I|}{|_0|}$$
(1)
which is the desired bound.

## 2 Numerical methods for the spherically-symmetric evolution

Recall from  that the SN equations are the system

$`i{\displaystyle \frac{\psi }{t}}`$ $`=`$ $`^2\psi +\varphi \psi ,`$ (2)
$`^2\varphi `$ $`=`$ $`|\psi |^2.`$
In this section and the next, we are concerned only with the spherically-symmetric case and in this case (2) can be simplified as

$`i{\displaystyle \frac{u}{t}}`$ $`=`$ $`{\displaystyle \frac{{}_{}{}^{2}u}{r^2}}+\varphi u,`$ (3)
$`{\displaystyle \frac{^2(r\varphi )}{r^2}}`$ $`=`$ $`{\displaystyle \frac{|u|^2}{r}},`$ (4)
where $`u=r\psi `$.

To solve this system, we shall use a Crank-Nicholson method for the time-evolution of the Schrödinger equation and the spectral method of  for the Poisson equation. Schematically, this takes the form

$`2i{\displaystyle \frac{u^{n+1}u^n}{\delta t}}`$ $`=`$ $`D^2(u^{n+1})D^2(u^n)`$ (5)
$`+\varphi ^{n+1}u^{n+1}+\varphi ^nu^n,`$
where, with the spectral method of , $`D^2`$ is the second-derivative matrix for Chebyshev polynomials. The Crank-Nicholson method is second-order accurate in the time and preserves the normalisation of the wave function.

Now we need an iterative method to find $`\varphi ^{n+1}`$. We write $`\varphi _k^{n+1}`$ and $`u_k^{n+1}`$ for the iterates and make the initial choice $`\varphi _0^{n+1}=\varphi ^n`$, then use the spectral method to solve (5) for $`u_0^{n+1}`$, use this in (4) to improve $`\varphi _0^{n+1}`$ to $`\varphi _1^{n+1}`$, and solve (5) to find $`u_1^{n+1}`$. This cycle can be repeated until the desired degree of convergence is reached.

The boundary conditions are that $`u`$ and $`r\varphi `$ should vanish at the origin and the outer edge $`r=L`$ of the grid. We expect there to be an outgoing flux of probability and we need to prevent it from reflecting back off of the boundary. For this we use the ‘sponges’ of Bernstein et al . The idea is to replace the $`i`$ on the left-hand-side in (3) by $`i+s(r)`$ where $`s(r)`$ is a positive function large (of order one) near the boundary but small close in. We typically use $`s(r)=ae^{b(rL)}`$ for constants $`a`$ and $`b`$. By this means the Schrödinger equation is converted to a heat equation near the boundary and the outgoing flux is absorbed.

We can perform various checks on the method outlined here. To check that the sponges are working as desired, we temporarily set $`\varphi =0`$ and consider the numerical evolution corresponding to the following explicit spherically-symmetric solution of the zero-potential Schrödinger equation:

$`r\psi =u={\displaystyle \frac{C\sqrt{\sigma }}{(\sigma ^2+2it)^{\frac{1}{2}}}}[\mathrm{exp}({\displaystyle \frac{(rvta)^2}{2(\sigma ^2+2it)}}+{\displaystyle \frac{ivr}{2}}{\displaystyle \frac{iv^2t}{4}})`$
$`\mathrm{exp}({\displaystyle \frac{(r+vt+a)^2}{2(\sigma ^2+2it)}}{\displaystyle \frac{ivr}{2}}{\displaystyle \frac{iv^2t}{4}})]`$ (6)
where $`C`$ is a normalisation constant (which can be found explicitly). This solution is a spherically-symmetric bump with a Gaussian profile centred initially at $`r=a`$ and moving radially with velocity $`v`$. It both disperses and moves off the grid. Computing it is a test for the Schrödinger evolution, and having it move smoothly off the grid is a test for the sponges. The initial data for the solution (6) will be evolved in the full evolution in Section 3.

Next, we can include a fixed $`\varphi `$ in (3), dropping (4) and jointly test the Schrödinger evolution and the sponges. Finally, we can test convergence of the method by varying the time step and/or the number of Chebyshev points chosen in the spectral method.

## 3 Results in the spherically-symmetric case

Having tested the method and the sponges satisfactorily, we evolve the ground state, that is we use as initial data the stationary-state of lowest energy which is known from  or . Since this is a stationary state it should evolve with a factor $`e^{iEt}`$ and this is just what we find (see figure 1: the phase is linear in time).

However if this evolution is continued for a long enough time, numerical errors accumulate and the solution drifts away from the ground state. By Fourier analysing the solution one finds lines in the power spectrum at (approximately) the normal frequencies for a perturbation about the ground state already found by linear perturbation theory in . Thus we really are getting the stable but perturbed ground state. We have also evolved from data which are the ground state with a small perturbation introduced by hand and the picture is the same.

Next we evolve the second spherically-symmetric stationary state taken from . We should expect this also to evolve as a stationary state at least for a while until numerical errors accumulate. What we obtain is shown in figure 2. This shows $`|r\psi |`$, so that the second state, which has a zero, appears as bimodal. The evolution proceeds as a stationary-state to about $`t=2000`$ and then there is a sudden change to a unimodal solution with some oscillation. This we claim is a perturbed and rescaled ground state solution, with total probability less than one.

To support this claim, we have again calculated the power spectrum and find lines at (approximately) the normal frequencies of perturbations about the ground state (to find agreement we need to rescale these frequencies with the factor corresponding to the rescaling of the residual ground state). We can also use the check suggested by equation (1). That is, we calculate both the probability $`p`$ and the action or conserved energy $``$ remaining on the grid at time $`t`$ and compute the bound on $`p`$ provided by (1). If the argument leading to this bound is correct then these two probabilities should converge on each other and as we see in figure 3 they do.

Again, we can evolve from the second state with a small perurbation introduced by hand rather than waiting for the numerical errors to accumulate when the collapse to the ground state is more rapid, or alternatively we can evolve from the second state but with a more stringent tolerance on the iteration determining $`\varphi `$, when the collapse takes longer.

Moving on, the decay of the third state is shown in figure 4. The initial state this time has three maxima and evolves for a while as a stationary state before collapsing to a rescaled ground state with the emission of some probability to infinity.

We have plotted the figure corresponding to figure 3 and we again find convergence.

The picture that emerges from these calculations is of nonlinear instability of all states after the ground state. Each higher state decays to the ground state either by the accumulation of numerical errors or because of an explicitly included perturbation. The end result is a noisy rescaled ground state, noisy because the eigenvalues for perturbation about the ground state are all imaginary.

The other spherically symmetric evolutions which we have calculated are with the initial data furnished by equation (6) with $`t=0`$. This is a Gaussian bump centered at $`r=a`$ with width $`\sigma `$ and moving with velocity $`v`$. The conserved energy of this data rises with $`v`$ (because the kinetic energy rises) and with $`\sigma `$ (because the particle is more localised). In all cases the solution disperses leaving a rescaled ground state which we can characterise by the probability residing in it. In figures 5, 6 and 7 we plot the residual probability against $`v`$, $`a`$ and $`\sigma `$ at different times.

From figure 5 we see that the residual probability is greatest if the lump is released from rest, but it is easier to trap probability if the initial velocity is ingoing than if it is outgoing. From figure 6 we see that for lumps released from rest, the residual probability is greater if the lump starts closer in, or in other words is more gravitationally bound. From figure 7 we see that reducing $`\sigma `$, which raises the energy, leads to more dispersion.

A key motivation for this set of calculations was the desire to check convergence of the method, and we have done this in a variety of ways. First for the time-step, with a Crank-Nicholson method we expect to have quadratic convergence, and this can be investigated with a Richardson quotient. We suppose that for some variable of interest the calculated value $`O_h`$ and actual values $`I`$ are related by

$$O_hI+Ah^k,$$

where the time step size is $`h`$, and $`k`$ is the order of the error, then we can calculate the Richardson quotient:

$$\frac{O_{h_1}O_{h_3}}{O_{h_2}O_{h_3}},$$

where $`h_1`$,$`h_2`$ and $`h_3`$ are three different values of the step size. With $`k`$ known, this will be a simple function of the $`h_i`$ and so will provide a check on $`k`$. Next for convergence in space we can repeat the calculation with different numbers $`N`$ of Chebyshev points. For both checks, the results converge as required.

With this calculation we confirm the picture of the SN evolution which we have been claiming: the Schrödinger equation tends to disperse probability but the Newtonian gravitational attraction holds it together; all states after the ground state are unstable and the evolution in general leads to a dispersion of probability to infinity, leaving some residual probability in a rescaled ground state.

## 4 The axially-symmetric SN equations

In this chapter, we solve the SN equations for an axially-symmetric system in 3 dimensions. The wave-function is now a function of polar coordinates $`r`$ and $`\theta `$ but is independent of polar angle. We shall first find stationary solutions. These typically have nonzero total angular momentum and include a dipole-like solution which appears to minimise the energy among wave functions which are odd in (the usual) $`z`$. We then consider the time-dependent problem. In particular we evolve the dipole-like state, and it turns out to be nonlinearly unstable - the two regions of probability density attract each other and fall together leaving a multiple of the ground state as the evolutions in Section 2 did.

With $`u=r\psi `$ the system of equations (2) becomes:

$`iu_t`$ $`=`$ $`u_{rr}{\displaystyle \frac{1}{r^2\mathrm{sin}\theta }}(\mathrm{sin}\theta u_\theta )_\theta +\varphi u,`$ (7)
$`{\displaystyle \frac{|u|^2}{r}}`$ $`=`$ $`(r\varphi )_{rr}+{\displaystyle \frac{1}{r\mathrm{sin}\theta }}(\mathrm{sin}\theta \varphi _\theta )_\theta .`$ (8)
For stationary solutions, the left-hand-side of (7) is replaced by $`Eu`$. Boundary conditions are that $`u=0`$ at $`r=0`$ and $`r=L`$, $`\varphi =0`$ at $`r=L`$ and finite at $`r=0`$, and $`u_\theta =0=\varphi _\theta `$ at $`\theta =0`$ and $`\theta =\pi .`$

To find stationary solutions we proceed as follows:

1. take as an initial guess for the potential $`\varphi ={\displaystyle \frac{1}{(1+r)}}`$;
2. using this potential solve the time-independent Schrödinger equation;
3. select an eigenfunction in such a way that the procedure will converge to a stationary state (this needs trial and error);
4. calculate the potential due to the chosen eigenfunction;
5. take the new potential to be the one obtained from step 4 above but symmetrised, since we require that $`\varphi `$ should be symmetric around the $`\theta ={\displaystyle \frac{\pi }{2}}`$ (otherwise numerical errors can cause the wave-function to move along the axis);
6. now provided that the new potential does not differ from the previous potential by some fixed tolerance in the norm of the difference, stop, otherwise continue from step 2.

In step 2 we solve the eigenvalue problem by a 2-dimensional spectral method, using Chebyshev differentiation in the directions of $`r`$ and $`\theta `$. In step 3, the eigenfunctions at the first iteration are labelled by the usual $`l,m`$ and $`n`$ quantum numbers, though with $`m=0`$ for axisymmetry. The idea is to choose one and run through the iteration, hoping for convergence. When the method converges, we do obtain a stationary state but we do not arrive at a one-to-one correspondence between solutions of the starting linear problem and the final nonlinear problem.

For the energy eigenvalue $`E`$ (which is not the conserved energy) we have the formula

$$E=|\psi |^2+\varphi |\psi |^2,$$
(9)
while the total angular momentum $`J^2`$, after integration by parts, is :

$$J^2=|\frac{\psi }{\theta }|^2$$
(10)
Both integrals are over $`𝐑^3`$.

We present in Table 1 the first few stationary states of the axially symmetric solution of the SN equations ordered by their energy and named for convenience axi1 to axi8.

In this table, axi1, axi3 and axi8 are spherically-symmetric solutions turning up again (as they should). Axi2, shown as a surface in figure 8 and as a contour plot in figure 9, is very much like a dipole solution and appears to be the solution minimising the energy among wavefunctions odd in $`z`$. If so it was found already in . Contour plots of axi4-axi7 are given in figures 10 to 13.

To solve the time-dependent axisymmetric SN-equations, we shall use an alternating direction implicit (or ADI) method (see e.g.,  or ) We split the Laplacian in (7) to write it as

$$\dot{u}=i(L_1+L_2\varphi )u$$
(11)
where

$`L_1`$ $`=`$ $`{\displaystyle \frac{1}{r}}{\displaystyle \frac{^2}{r^2}}`$ (12)
$`L_2`$ $`=`$ $`{\displaystyle \frac{1}{r}}[{\displaystyle \frac{^2}{\theta ^2}}+\mathrm{cot}\theta {\displaystyle \frac{}{\theta }}].`$ (13)
Next we introduce new variables $`S`$ and $`T`$ and write this in the formally equivalent form:

$`\mathrm{exp}({\displaystyle \frac{ih}{2}}L_1)S(t)`$ $`=`$ $`\mathrm{exp}({\displaystyle \frac{ih}{2}}L_2)u(t)`$
$`\mathrm{exp}({\displaystyle \frac{ih}{2}}L_2)T(t)`$ $`=`$ $`\mathrm{exp}({\displaystyle \frac{ih}{2}}L_1)S(t)`$ (14)
$`\mathrm{exp}({\displaystyle \frac{ih}{2}}\varphi )u(t+h)`$ $`=`$ $`\mathrm{exp}({\displaystyle \frac{ih}{2}}\varphi )T(t)`$
where we suppress the dependence on $`r`$ and $`\theta `$. To obtain a discrete form of (14) we linearise in the time-step $`h`$ to find

$`(1{\displaystyle \frac{ih}{2}}L_1)S^n`$ $`=`$ $`(1+{\displaystyle \frac{ih}{2}}L_2)u^n`$
$`(1{\displaystyle \frac{ih}{2}}L_2)T^n`$ $`=`$ $`(1+{\displaystyle \frac{ih}{2}}L_1)S^n`$ (15)
$`(1+{\displaystyle \frac{ih}{2}}\varphi ^{n+1})u^{n+1}`$ $`=`$ $`(1{\displaystyle \frac{ih}{2}}\varphi ^n)T^n`$
where the superscript indicates the value of discretised time.

For the Poisson equation (8) we have

$$(L_1+L_2)(r\varphi )=\frac{1}{r^2}|u|^2$$

which we solve at each instant by a Peaceman-Rachford ADI iteration . This means we discretise it as

$`(L_1+\rho )\varphi _{k+1}`$ $`=`$ $`(L_2\rho )\varphi _k+{\displaystyle \frac{1}{r^2}}|u|^2`$ (16)
$`(L_2+\rho )\varphi _{k+1}`$ $`=`$ $`(L_1\rho )\varphi _k+{\displaystyle \frac{1}{r^2}}|u|^2`$
and iterate with $`\varphi ^0=0`$. Here $`k`$ labels the iteration (at a fixed time) and $`\rho `$ is a small constant chosen so that the iteration converges.

The boundary conditions are that $`u`$ and $`\varphi `$ vanish at the outer boundary and so we need sponges as in Section 2 to prevent waves of probability reflecting back.

As an example, we take as initial data the dipole state of figure 8. The evolution of $`|\psi |`$ for this is shown in figure 14. (Note that $`\psi `$ is initially real but must become complex before again becoming approximately real in the remote future. The initial data for $`\psi `$ is odd as a function of $`z`$ but it ends up approximately even.)

Following the lessons learned in Section 3, we should expect the stationary state to be unstable. We can compute the evolution for short times to see that initially it remains a stationary state, that is that the only change is a phase growing linearly with time. However when numerical errors have built up the state becomes unstable and the two concentrations of probability density fall into each other. Probability leaves the grid, as does angular momentum, and we are left with a multiple of the ground state.

We can check convergence of the method by calculating a Richardson quotient with three different time-steps (to find that convergence is now linear in time), and by varying the number of Chebyshev points.

What we learn from this calculation is that, as well as dispersion which is what we mostly saw in Section 3, the solutions of the SN equation show gravitational attraction: lumps of probability density released from rest fall into each other. In the next section we shall see that, at least for the two-dimensional SN equations, lumps of probability can orbit each other.

## 5 The two-Dimensional SN equations

In this section, we shall consider the SN equations in a plane, that is in Cartesian coordinates $`x`$,$`y`$. We shall find a dipole-like stationary solution, and some solutions which are like rigidly rotating dipoles. These rigidly rotating solutions are unstable however and will merge, radiating angular momentum.

The SN equations in this case are

$`i\psi _t`$ $`=`$ $`\psi _{xx}\psi _{yy}+\varphi \psi `$ (17)
$`\varphi _{xx}+\varphi _{yy}`$ $`=`$ $`|\psi |^2`$
We use the ADI scheme as in equation (15) but with the understanding that now

$`L_1`$ $`=`$ $`{\displaystyle \frac{^2}{x^2}}`$
$`L_2`$ $`=`$ $`{\displaystyle \frac{^2}{y^2}}.`$
For the potential we use the counterpart of (16) with the same understanding. The boundary conditions are that $`\psi `$ and $`\varphi `$ vanish at the edges of the grid, which is now a large square. We still need sponges and we do not want them to have corners so we take functions like $`s(x,y)=\mathrm{min}[1,e^{0.5(\sqrt{(}x^2+y^2)20)}]`$. We can test the efficacy of the sponges by evolving moving two-dimensional Gaussians to see that they propagate off the grid, as they do.

Once we have the code, we can look for solutions corresponding to those found already in 3-dimensions. In particular we find a stable ground state and a dipole-like solution which is stable for a while before decaying to the ground state. We can seek solutions with no counterpart among those found already by making an ansatz of rigid rotation. That is we look for solutions which in polar coordinates $`r,\theta `$ take the form:

$$\psi (r,\theta ,t)=e^{iEt}\psi (r,\theta +\omega t)$$
(18)
for (real) constants $`E`$ and $`\omega `$. (Solutions of a similar nature in 3-dimensions would depend on all three spatial coordinates which is why we haven’t seen them so far.)

To separate the SN equations for a solution like (18), we go into the rotating frame with Cartesian coordinates $`X,Y`$ given by

$`X=x\mathrm{cos}\omega t+y\mathrm{sin}\omega t,`$
$`Y=x\mathrm{sin}\omega t+y\mathrm{cos}\omega t,`$
and the time-dependence separates off to leave the equations as

$`^2\psi +\varphi \psi i\omega (Y\psi _XX\psi _Y)`$ $`=`$ $`E\psi `$ (19)
$`^2\varphi `$ $`=`$ $`|\psi |^2,`$
To solve (19), we begin with a small value of $`\omega `$ and the dipole-like solution. Iteration leads to a solution, and we can study its change with increasing $`\omega `$. The wave-function $`\psi `$ is necessarily complex and we display the real and imaginary parts of such a solution, with $`\omega =0.005`$, in figures 15 and 16.

The solution just found is a stationary state but we must expect it to be unstable. If we use it as initial data for the time-dependent problem then we find as before that it evolves as a stationary state for some time (see figure 17) before numerical errors build up and it becomes unstable (figure 18). The two lumps of probability orbit each other about four times before collapsing into a single lump in about one orbital period. Probability and angular momentum are radiated off of the grid.

What we have found by this calculation is a confirmation of the picture of lumps of probability interacting gravitationally with each other. Here the lumps are in orbit around each other for a while before becoming unstable and collapsing into a single lump. By earlier work, we must then expect the probability to disperse leaving a rescaled ground state, as it does.

## Acknowledgement

The work described in this paper formed part of the D.Phil. thesis of the first author and he gratefully acknowledges the receipt of a grant from EPSRC.


format: text

text:
# Tau neutrino deep inelastic charged current interactions

## I Introduction

Results from the SuperKamiokande underground experiment measuring the atmospheric neutrino flux suggest that muon neutrinos oscillate into tau neutrinos with nearly maximal mixing . A test of the oscillation hypothesis is $`\nu _\tau `$ production of $`\tau `$ through charged current interactions, a process which will be studied in underground neutrino telescopes  as well as long-baseline experiments  measuring neutrino fluxes from accelerator sources. For precision measurements of oscillation mixing angles and eventually CP violation, neutrino cross sections will ideally be know to the level of a few percent. Eventually, measurements of neutrino-nucleon charged current interaction cross sections are expected to be at the 1% level at a neutrino factory .

The QCD theory of deep inelastic cross sections in the leading-twist approximation has proceeded to the point of evaluating next-to-next-to leading order (2 loop) perturbative QCD corrections to the coefficient functions for the structure functions $`F_1,F_2`$ and $`F_3`$  and approximations for the splitting functions . In the specific case of $`\nu _\mu `$-isoscalar nucleon ($`N`$) cross sections, target mass corrections and nuclear binding effects in a LO twist-2 approach of DIS, combined with the elastic peak and a modeling of higher twists in the continuum and resonance region of DIS, have also been investigated . Tau neutrino charged current interactions with nucleons have received less theoretical attention . Albright and Jarlskog, in Ref. , pointed out that there are two additional structure functions, $`F_4`$ and $`F_5`$ that contribute to the tau neutrino cross section. $`F_4`$ and $`F_5`$ are ignored in muon neutrino interactions because of a suppression factor depending on the square of the charged lepton mass ($`m_{\mathrm{}}`$) divided by the nucleon mass times neutrino energy, $`m_{\mathrm{}}^2/(M_NE_\nu )`$. At leading order, in the limit of massless quarks and target hadrons, $`F_4`$ and $`F_5`$ are

$`F_4`$ $`=0`$ (1)
$`2xF_5`$ $`=F_2,`$ (2)
where $`x`$ is the Bjorken-$`x`$ variable. These generalizations of the Callan-Gross relation $`F_2=2xF_1`$ are called the Albright-Jarlskog relations. As with the Callan-Gross relations, the Albright-Jarlskog relations are violated from kinematic mass corrections and at NLO***We will find below that Eq. (2) is not violated in massless NLO QCD. in QCD. We quantify the violation in Section III A.

Below $`E_\nu 10`$ GeV for muon neutrinos and higher energies for tau neutrinos, untangling the exclusive and inclusive contributions to the neutrino-nucleon cross section is difficult. The cross section in this energy range has quasi-elastic, resonant production (such as $`\mathrm{\Delta }`$ production), non-resonant pion production and other inelastic contributions. There are several phenomenological methods for avoiding double counting in this region. One method employs a cutoff in the hadronic invariant mass $`W>W_{\mathrm{min}}`$ such that the combined exclusive and $`W_{\mathrm{min}}`$ dependent inclusive cross section yields a total consistent with data . Using muon neutrino data, a common choice for $`W_{\mathrm{min}}`$ is 1.4 GeV. A second method involves looking at different final state multiplicities . By normalizing the calculated total cross section for a given multiplicity $`j`$, $`\sigma _{\mathrm{tot}}^j`$, to the data, one can determine the factors $`f_j`$ in $`\sigma _{\mathrm{tot}}^j=\sigma _{\mathrm{res}}^j+f_j\sigma _{\mathrm{DIS}}^j`$. The quantity $`\sigma _{\mathrm{DIS}}^j`$ is determined by the inelastic cross section and the hadronization scheme. The Monte Carlo for the Soudan experiment does the normalization at $`E_\nu =20`$ GeV .

We present here the deep-inelastic contribution to $`\nu _\tau N\tau X`$ incorporating next-to-leading order QCD corrections, collinear target mass and charmed quark mass corrections. The NLO structure functions $`F_4`$ and $`F_5`$ including charm quark production have been evaluated by Gottschalk in Ref. . Here, we correct a misprint and evaluate gluonic helicity states in $`D`$ dimensions rather than 4 dimensions. We show numerical results for the new structure functions and for the charged current cross sections, both with and without imposing a nonzero value for $`W_{\mathrm{min}}`$. To gauge the effects of perturbative delicacies such as higher-twist beyond the inclusion of target mass effects in the scaling variable, we also evaluate the effect of imposing a cutoff on $`Q^2`$. We find that the collinear cutoff effect of the tau mass $`m_\tau =1.78\mathrm{GeV}`$ is enough to guarantee that for $`E_{\nu _\tau }5\mathrm{GeV}`$, $`\sigma _{\mathrm{CC}}(\nu _\tau N)`$ is dominated by $`Q^2>1\mathrm{GeV}^2`$. While this is widely considered a perturbative scale in QCD, some caveats about the possible importance of higher twist are discussed.

In Section II, we show the formulae for the differential cross section in terms of structure functions, and the expressions for the structure functions at NLO. Charm quark mass corrections for the charged current process are displayed, together with the results for the $`m_c0`$ limit. In Section III, we exhibit our numerical results for the structure functions $`F_4`$ and $`F_5`$, and for $`\nu _\tau N`$ and $`\overline{\nu }_\tau N`$ interactions for neutrino energies up to 100 GeV. Cross sections are compared to those with incident muon neutrinos and uncertainties are discussed. We conclude in Section IV. In an Appendix we rewrite the Albright-Jarlskog relations in terms of helicity amplitudes and pinpoint the approximations in their derivation.

## II Formulary: $`\nu _\tau `$ DIS at NLO

Neglecting neither the target nucleon mass $`M_N`$ nor the final state lepton mass $`m_\tau `$, the charged current $`\nu _\tau `$ (anti-)neutrino differential cross section is represented by a standard set of 5 structure functions  Our normalization of $`F_4`$ differs from that in  by a factor of $`x`$.

$`{\displaystyle \frac{d^2\sigma ^{\nu (\overline{\nu })}}{dxdy}}`$ $`=`$ $`{\displaystyle \frac{G_F^2M_NE_\nu }{\pi (1+Q^2/M_W^2)^2}}\{(y^2x+{\displaystyle \frac{m_\tau ^2y}{2E_\nu M_N}})F_1^{W^\pm }`$ (3)
$`+`$ $`\left[(1{\displaystyle \frac{m_\tau ^2}{4E_\nu ^2}})(1+{\displaystyle \frac{M_Nx}{2E_\nu }})y\right]F_2^{W^\pm }`$ (4)
$`\pm `$ $`\left[xy(1{\displaystyle \frac{y}{2}}){\displaystyle \frac{m_\tau ^2y}{4E_\nu M_N}}\right]F_3^{W^\pm }`$ (5)
$`+`$ $`{\displaystyle \frac{m_\tau ^2(m_\tau ^2+Q^2)}{4E_\nu ^2M_N^2x}}F_4^{W^\pm }{\displaystyle \frac{m_\tau ^2}{E_\nu M_N}}F_5^{W^\pm }\}.`$ (6)
where $`\{x,y,Q^2\}`$ are the standard DIS kinematic variables related through $`Q^2=2M_NE_\nu xy`$ and where we have neglected factors of $`m_\tau ^2/2M_NE_\nu Q^2/M_W^2`$. These latter corrections from the $`q^\mu q^\nu /M_W^2`$ part of the massive boson propagator are negligible both at low and at high neutrino energies and will not enter our numerics. For completeness of this formulary, though, they can be included multiplicatively by replacing:

$$F_i^{W^\pm }F_i^{W^\pm }\times \left(1+ϵ_i\right)$$
(7)
with

$`ϵ_1`$ $`=`$ $`{\displaystyle \frac{m_\tau ^2(Q^2+2M_W^2)}{2M_W^4}}`$ (8)
$`ϵ_2`$ $`=`$ $`{\displaystyle \frac{E_\nu ^2m_\tau ^2y[4M_W^2+y(Q^2+m_\tau ^2)]}{M_W^4[4(y1)E_\nu ^2+m_\tau ^2+Q^2]}}`$ (9)
$`ϵ_3`$ $`=`$ $`0`$ (10)
$`ϵ_4`$ $`=`$ $`{\displaystyle \frac{Q^2(Q^2+2M_W^2)}{M_W^4}}`$ (11)
$`ϵ_5`$ $`=`$ $`{\displaystyle \frac{Q^2}{M_W^2}}+{\displaystyle \frac{(M_W^2+Q^2)(m_\tau ^2+Q^2)y}{2M_W^4}}`$ (12)
Kinematics determine the integration rangesA typographic error in Refs.  is corrected here in the lower limit for $`x`$. to be

$$\frac{m_\tau ^2}{2M_N(E_\nu m_\tau )}x1$$
(13)
and

$$abya+b$$
(14)
where

$`a`$ $`={\displaystyle \frac{1m_\tau ^2\left(\frac{1}{2M_NE_\nu x}+\frac{1}{2E_\nu ^2}\right)}{2\left(1+\frac{M_Nx}{2E_\nu }\right)}}`$ (15)
$`b`$ $`={\displaystyle \frac{\sqrt{\left(1\frac{m_\tau ^2}{2M_NE_\nu x}\right)^2\frac{m_\tau ^2}{E_\nu ^2}}}{2\left(1+\frac{M_Nx}{2E_\nu }\right)}}.`$ (16)
In the perturbative regime we can calculate the structure functions $`F_i^{W^\pm }`$ from parton dynamics by use of the factorization theorem

$$W^{\mu \nu }=\frac{d\xi }{\xi }f(\xi ,\mu ^2)\widehat{\omega }^{\mu \nu }|_{p^+=\xi P_N^+}.$$
(17)
relating the hadronic ($`W^{\mu \nu }`$) and partonic ($`\widehat{\omega }^{\mu \nu }`$) forward matrix element of the product of weak currents $`<J^\mu J^\nu >`$. In Eq. (17), $`f(\xi ,\mu ^2)`$ is a parton distribution function evaluated at factorization scale $`\mu `$ and the parton momentum fraction $`\xi `$ of the light cone momentum of the nucleon $`P_N^+(P_N^0+P_N^z)/\sqrt{2}`$. Intrinsic transverse momentum of the incoming parton is neglected throughout our discussion, i.e. $`p_{}=0`$ in (17).

For the neutrino energies of interest here, we can safely restrict the consideration to the first two quark generations. The light flavour contributions to Eq. (6) can be obtained from the $`m_c0`$ limit of the charm production component which we will, therefore, consider first. The charm production contribution to $`d\sigma `$ will be represented by the structure functions $`F_i^c`$. We introduce theoretical structure functions§§§To meet with an experimentalist’s point of view, we denote functional dependence on $`x`$ and $`Q`$ only in the LHS of Eq. (19); corresponding in theory to an exact knowledge of all mass parameters and an ideal validity of the renormalization group for observable cross sections. Residual scale ($`\mu `$) dependence will be investigated below.

$`_i^c(x,Q^2)=(1\delta _{i4})s^{}(\overline{\eta },\mu ^2)`$ $`+`$ $`{\displaystyle \frac{\alpha _s(\mu ^2)}{2\pi }}\{{\displaystyle _{\overline{\eta }}^1}{\displaystyle \frac{d\xi ^{}}{\xi ^{}}}[H_i^q(\xi ^{},{\displaystyle \frac{Q^2}{\mu ^2}},\lambda )s^{}({\displaystyle \frac{\overline{\eta }}{\xi ^{}}},\mu ^2)`$ (18)
$`+`$ $`H_i^g(\xi ^{},{\displaystyle \frac{Q^2}{\mu ^2}},\lambda )g^{}({\displaystyle \frac{\overline{\eta }}{\xi ^{}}},\mu ^2)]\}`$ (19)
for scattering off the CKM-rotated weak eigenstate

$$s^{}=\left|V_{s,c}\right|^2s+\left|V_{d,c}\right|^2d$$
(20)
and its QCD evolution partner

$$g^{}g\left(\left|V_{s,c}\right|^2+\left|V_{d,c}\right|^2\right)$$
(21)
i.e., $`ds^{}/d\mathrm{ln}Q^2=s^{}P_{qq}+g^{}P_{qg}`$.

In Eq. (19) we set the renormalization scale equal to the factorization scale and

$$\overline{\eta }=\frac{\eta }{\lambda },\frac{1}{\eta }=\frac{1}{2x}+\sqrt{\frac{1}{4x^2}+\frac{M^2}{Q^2}}$$
(22)
is the target mass corrected slow rescaling variable.The quantity $`\eta `$ is the Nachtmann variable . The charm mass dependence is included in the dimensionless $`\lambda Q^2/(Q^2+m_c^2)`$. In Eq. (19) the convolution variable $`\xi `$ in Eq. (17) has been traded for $`\xi ^{}=\overline{\eta }/\xi `$ which relates to the partonic CMS energy $`\widehat{s}=(p+q)^2`$ through $`1/\xi ^{}=\lambda (1+\widehat{s}/Q^2)`$.

The theoretical structure functions in Eq. (19) are obtained from a tensor projection on the partonic $`\widehat{\omega }^{\mu \nu }`$ and have been conveniently normalized to a simple $`𝒪(\alpha _s^0)`$ term. As indicated, the leading order contribution to $`_4^c`$ vanishes. In the presence of target mass, the physical structure functions $`F_i^{W^\pm }`$ in Eq. (6) \- relating to a tensor projection on the hadronic $`W^{\mu \nu }`$ \- are a mixtureSee Ref.  for details on the mixing of structure functions within the tensor basis. of $`_i`$. Explicitly in our case of CC charm production we have:

$`F_1^c`$ $`=`$ $`_1^c`$ (23)
$`F_2^c`$ $`=`$ $`2{\displaystyle \frac{x}{\lambda }}{\displaystyle \frac{_2^c}{\rho ^2}}`$ (24)
$`F_3^c`$ $`=`$ $`2{\displaystyle \frac{_3^c}{\rho }}`$ (25)
$`F_4^c`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}{\displaystyle \frac{(1\rho )^2}{2\rho ^2}}_2^c+_4^c+{\displaystyle \frac{1\rho }{\rho }}_5^c`$ (26)
$`F_5^c`$ $`=`$ $`{\displaystyle \frac{_5^c}{\rho }}{\displaystyle \frac{(\rho 1)}{\lambda \rho ^2}}_2^c`$ (27)
where

$$\rho ^21+\left(\frac{2M_Nx}{Q}\right)^2,$$
(28)
$`F_i^c=F_i^c(x,Q^2)`$ and $`_i^c=_i^c(x,Q^2)`$. These results are in agreement with Table V in Ref. . They rely only on the assignment of parton light cone momentum $`p^+`$ related to the nucleon light cone momentum $`p^+=\xi P_N^+`$, with $`p_{}=0`$ and $`\xi =\overline{\eta }`$ at leading order.

The NLO corrections $`H_{i=1,\mathrm{},5}^{q,g}`$ were first obtained in Ref.  and $`H_{1,2,3}^{q,g}`$ have been rederived in . We follow these references in notation and present the full set $`(i=1,\mathrm{},5)`$ including our independently rederived $`H_{4,5}^{q,g}`$. The fermionic NLO coefficient functions $`H_i^q`$ in Eq. (19) for charm production are calculated from the subprocess $`W^+sgc`$ and from the one-loop correction to $`W^+sc`$. They are given byFor the ease of notation, from here on $`\xi `$ will be a generic variable and no longer the light-cone momentum of Eq. (17). To strictly match Eq. (19) we should have $`\xi \xi ^{}`$ instead.

$$H_{i=1,2,3,5}^q(\xi ,\frac{Q^2}{\mu ^2},\lambda )=\left[P_{qq}^{(0)}(\xi )\mathrm{ln}\frac{Q^2}{\lambda \mu ^2}+h_i^q(\xi ,\lambda )\right]$$
(29)
where $`P_{qq}^{(0)}(\xi )={\displaystyle \frac{4}{3}}\left({\displaystyle \frac{1+\xi ^2}{1\xi }}\right)_+`$ and

$`h_i^q(\xi ,\lambda )={\displaystyle \frac{4}{3}}`$ $`\{`$ $`h^q+A_i\delta (1\xi )+B_{1,i}{\displaystyle \frac{1}{(1\xi )_+}}`$ (30)
$`+`$ $`B_{2,i}{\displaystyle \frac{1}{(1\lambda \xi )_+}}+B_{3,i}\left[{\displaystyle \frac{1\xi }{(1\lambda \xi )^2}}\right]_+\}`$ (31)
with

$`h^q=`$ $``$ $`\left(4+{\displaystyle \frac{1}{2\lambda }}+{\displaystyle \frac{\pi ^2}{3}}+{\displaystyle \frac{1+3\lambda }{2\lambda }}K_A\right)\delta (1\xi )`$ (32)
$``$ $`{\displaystyle \frac{(1+\xi ^2)\mathrm{ln}\xi }{1\xi }}+(1+\xi ^2)\left[{\displaystyle \frac{2\mathrm{ln}(1\xi )\mathrm{ln}(1\lambda \xi )}{1\xi }}\right]_+`$ (33)
and

$$K_A=\frac{1}{\lambda }(1\lambda )\mathrm{ln}(1\lambda ).$$
(34)
The coefficients in (31) for $`i=1,2,3,5`$ are given in Table I. A misprint in Ref.  concerning $`A_2`$ was already corrected in . Here, we correct a similar misprint concerning $`A_5`$.

The gluonic NLO coefficient functions $`H_i^g`$ in Eq. (19) for charm production, as calculated from the subprocess $`W^+gc\overline{s}`$, are given by

$$H_{i=\genfrac{}{}{0pt}{}{1,2,5}{3}}^g(\xi ,\frac{Q^2}{\mu ^2},\lambda )=\left[P_{qg}^{(0)}(\xi )\left(\pm L_\lambda +\mathrm{ln}\frac{Q^2}{\lambda \mu ^2}\right)+h_i^g(\xi ,\lambda )\right]$$
(35)
where $`P_{qg}^{(0)}(\xi )={\displaystyle \frac{1}{2}}\left[\xi ^2+(1\xi )^2\right],L_\lambda =\mathrm{ln}{\displaystyle \frac{1\lambda \xi }{(1\lambda )\xi }}`$ and

$$h_i^g(\xi ,\lambda )=C_0+C_{1,i}\xi (1\xi )+C_{2,i}+(1\lambda )\xi L_\lambda (C_{3,i}+\lambda \xi C_{4,i})$$
(36)
with

$$C_0=P_{qg}^{(0)}(\xi )\left[2\mathrm{ln}(1\xi )\mathrm{ln}(1\lambda \xi )\mathrm{ln}\xi \right].$$
(37)
The coefficients $`C_{k,i}`$ are given in Table II. As in Ref. , they differ from those in Ref.  by counting \- within dimensional regularization - the gluonic helicity states in $`D`$  rather than in $`4`$ dimensions.

The structure function $`F_4`$ is insensitive to any collinear physics up to $`𝒪(\alpha _s^1)`$ and the corresponding coefficients $`H_{i=4}^{q,g}`$ are, therefore, scheme and scale-independent functions. They are given by

$$H_{i=4}^q(\xi ,\lambda )=\frac{4}{3}\frac{\lambda (1\xi )\xi \left[1+(12\lambda )\xi \right]}{(1\lambda \xi )^2}$$
(38)
$$H_{i=4}^g(\xi ,\lambda )=2\lambda \xi \left[1\xi (1\lambda )\xi L_\lambda \right].$$
(39)
In the limit $`\lambda 1(m_c0)`$ and after an additional minimal subtraction of the collinear mass singularities in $`H_{i=1,2,3,5}^g`$ the $`H_i^{q,g}`$ reduce to the massless $`\overline{\mathrm{MS}}`$ coefficient functions:\**\**\**The choice $`\mu =Q`$ is made to match with the notation in . Retaining a general $`\mu Q`$ on the LHS of (40), (41) results in the massless coefficient functions for arbitrary scale. Obviously, this amounts to adding back the $`P^{(0)}\mathrm{ln}(Q^2/\mu ^2)`$ “splitting function times log” counter-terms.

$`\underset{\lambda 1}{lim}H_i^q(\xi ,{\displaystyle \frac{Q^2}{\mu }},\lambda )|_{\mu ^2=Q^2}`$ $`=`$ $`C_{F,i}^{(1)}(\xi )`$ (40)
$`\underset{\lambda 1}{lim}\left\{H_{i=\genfrac{}{}{0pt}{}{1,2,4,5}{3}}^g(\xi ,{\displaystyle \frac{Q^2}{\mu ^2}},\lambda )(1\delta _{i4})P_{qg}^{(0)}\mathrm{ln}{\displaystyle \frac{\mu ^2/Q^2}{1\lambda }}\right\}|_{\mu ^2=Q^2}`$ $`=`$ $`C_{G,i}^{(1)}(\xi )`$ (41)
where $`C_{F,i}^{(1)}`$ and $`C_{G,i}^{(1)}`$ for $`i=1,2,3`$ are the massless $`\overline{\mathrm{MS}}`$ coefficients in Appendix III of . Extending the results in Refs.  to include $`i=4,5`$ and within the notation of Ref.  we find

$`C_{F,4}^{(1)}(\xi )`$ $`=`$ $`{\displaystyle \frac{4}{3}}\xi `$ (42)
$`C_{F,5}^{(1)}(\xi )`$ $`=`$ $`C_{F,2}^{(1)}(\xi )`$ (43)
$`C_{G,4}^{(1)}(\xi )`$ $`=`$ $`2\xi (1\xi )`$ (44)
$`C_{G,5}^{(1)}(\xi )`$ $`=`$ $`C_{G,2}^{(1)}(\xi ).`$ (45)
To calculate the light quark contributions to CC $`\nu _\tau `$ DIS, the massless coefficient functions on the RHS of Eqs. (40)-(45) are used together with the PDFs which multiply the CKM matrix elements $`\left|V_{i,jc}\right|^2`$ in an obvious modification ($`\lambda 1`$, $`HC`$, …) of our Eq. (19). The result is added to the charm production component (19) to obtain the entire NLO structure function.

It is interesting to note that the equalities (43), (45) guarantee that the Albright-Jarlskog relation (2) is not violated in massless QCD at NLO. Eqs. (23)-(27) and the fact that $`H_2^{q,g}H_5^{q,g}`$ for charm production manifest, of course, a violation of Eq. (2) in the real world of massive target hadrons and of heavy quarks interacting through QCD. This observation must actually be expected to hold at all orders as we clarify in the Appendix.

## III Structure Functions and Cross Sections

### A Structure Functions

The Albright-Jarlskog relations (Eqs. (1) and (2)) are valid at leading order in the massless limit. Here, we show $`F_4`$ and $`2xF_5F_2`$ to demonstrate violations of these relations. We use the CTEQ6 parton distributions  which include estimates of the uncertainties in the distributions. Thus, in Section III B we will be able to quantify the error in the evaluation of $`\sigma _{\mathrm{CC}}(\nu _\tau N)`$ that is caused by our imperfect knowledge of the PDFs. The CTEQ6 fits are provided for $`\mu >\mu _0=1.3`$ GeV but can be extrapolated to lower values of $`\mu `$. When the factorization scale squared goes below $`\mu _{\mathrm{cut}}^2=0.5`$ GeV2, though, we choose to freeze it at $`\mu ^2=\mu _{\mathrm{cut}}^2`$. An alternate set of parton distribution functions used below is the Glück, Reya and Vogt GRV98 set  which evolves from $`\mu _0^2=0.4\mathrm{GeV}^2`$ with parametrizations provided above $`\mu _{\mathrm{cut}}^2=0.8`$ GeV2. In principle, CTEQ6 covers the parameter space of global PDF-related data within conservative errors. On the other hand, the evolution with three quark flavours ($`u,d,s`$) of GRV98 matches better with our evaluation of $`\sigma _{\mathrm{CC}}(\nu _\tau N)`$ at low neutrino energies where we assume light sea quarks only. In comparison, the CTEQ6 PDFs are evolved with a variable flavour number which is, strictly, not fully compatible with our approach. When we quantify the uncertainties from the PDF degrees of freedom we will, therefore, employ GRV98 for an independent comparison. Thereby, we can convince ourselves that the slight inconsistency mentioned above leads to no noticeable bias in practice. The CTEQ6 NLO $`\overline{\mathrm{MS}}`$ set of parton distribution functions is our default choice in the evaluations below. If not stated otherwise, the curves correspond to $`\mu =Q`$.

In Fig. 1, we show $`F_4(x,Q^2=2\mathrm{GeV}^2)`$. The solid lines in the figure include target mass corrections of Eq. (26), while the dashed lines are with the target mass set to zero, so that $`\eta x`$ and $`\rho =1`$. The leading order (LO) curve with $`M_N=0`$ shows $`F_4=0`$, even with charm mass corrections included. A leading order violation would occur if the initial quark masses were set to non-zero values. Including the target mass corrections shows the effect of the mixing of $`_4`$ with $`_2`$ and $`_5`$ in Eq. (26). The NLO corrections have an effect primarily at small-$`x`$.

Fig. 2 shows the violation of the second Albright-Jarskog relation, at fixed $`Q^2=2`$ GeV2. Even at leading order, $`2xF_5F_20`$. This is due to charm quark mass corrections in $`W^+s^{}c`$. The magnitude at small-$`x`$ reflects the impact of the $`s^{}`$ sea distribution. Target mass effects incorporated by Eq. (27) are not significant. Including the NLO corrections makes small changes to the curve at small-$`x`$. Both of the figures show that in evaluations of the total charged current cross section, the naive Albright-Jarlskog relations are good approximations to the NLO results. This is true at low energies, where $`\sigma _{\mathrm{CC}}(\nu _\tau N)`$ does not probe small-$`x`$ and at high energies where $`F_{4,5}`$ are suppressed, anyway.

### B Cross Sections

The cross sections for neutrino and antineutrino charged current interactions with isoscalar nucleons are shown in the first two panels of Fig. 3. We use our default set of CTEQ6 parton distribution functions with the factorization scale equal to the renormalization scale $`\mu ^2`$, and $`\mu ^2=Q^2`$. Below $`\mu ^2=\mu _{\mathrm{cut}}^2`$, we set $`\mu ^2=\mu _{\mathrm{cut}}^2`$ in the parton distribution functions and in $`\alpha _s`$, but we keep the explicit $`Q^2`$ dependence in the differential cross section of Eq. (6) as well as in the counter-logs $`\mathrm{ln}Q^2/\mu _{\mathrm{cut}}^2`$ of the NLO coefficient functions. For $`Q^2<\mu _{\mathrm{cut}}^2`$, a noticable impact of these logs is indicative of long-distance strong interaction. The evaluation can, therefore, not be trusted perturbatively, whenever it becomes sensitive to the technical choice $`\mu >\mu _{\mathrm{cut}}`$. Most of our results below are, however, completely insensitive to it. In Fig. 3 muon (anti-)neutrino cross sections appear with dashed lines, while the solid lines show the tau (anti-)neutrino cross sections. The upper curves show no cuts while the lower curves have $`W^2=Q^2(x^1+1)+M_N^2>(1.4\mathrm{GeV})^2`$ and $`Q^2>1`$ GeV2. As we show in Fig. 4, for $`W_{\mathrm{min}}=1.4`$ GeV, the tau neutrino CC cross section is fairly insensitive to the $`Q^2`$ cut of 1 GeV2. At, for example, $`E_\nu =20`$ GeV, $`\sigma _{\mathrm{CC}}(Q_{\mathrm{min}}^2=1\mathrm{GeV}^2)/\sigma _{\mathrm{CC}}(Q_{\mathrm{min}}^2=0)=0.93`$ for $`\nu _\tau N`$ CC interactions. The $`Q^2`$ cut has a larger impact on $`\sigma _{\mathrm{CC}}(\nu _\mu N)`$, where $`\sigma _{\mathrm{CC}}(Q_{\mathrm{min}}^2=1\mathrm{GeV}^2)/\sigma _{\mathrm{CC}}(Q_{\mathrm{min}}^2=0)=0.85`$ for $`E_\nu =20`$ GeV, due to the fact that for nearly massless leptons, $`Q^2>Q_{\mathrm{min}}^2`$ cuts out a larger fraction of the available phase space.

The small changes in the CC cross sections with $`Q_{\mathrm{min}}^2=1\mathrm{GeV}^2`$ lead one to expect that non-perturbative effects at low $`Q^2`$ are unlikely to be large when one also applies the $`W_{\mathrm{min}}`$ cut of 1.4 GeV. At low energies, without the $`W_{\mathrm{min}}`$ cut, a substantial contribution to the cross sections comes from $`Q^2<1`$ GeV2, however, and this is precisely where the DIS cross section is only a rough approximation to the true cross sections with quasi-elastic and resonant as well as non-resonant contributions.

The bottom panel of Fig. 3 shows the ratio of the $`\nu _\tau N`$ to $`\nu _\mu N`$ CC cross sections (solid lines) and the same ratio for antineutrinos. Shown are the uncut results, but the results with the $`W_{\mathrm{min}}`$ and $`Q^2`$ cuts agree to within 3% for $`E_\nu >20`$ GeV. Of note is the fact that even at $`E_\nu =10^3`$ GeV, the $`\nu _\tau N`$ to $`\nu _\mu N`$ CC cross section ratio is still 5% below unity. At 100 GeV, the ratio is 0.76. There are two reasons for the deficit in the $`\nu _\tau `$ CC cross section: the reduced phase space and the contribution of $`F_5`$. The reduced phase space is reflected in the integration limits for $`x`$ and $`y`$ (Eqs. (13) and (14)) . This is responsible for about half of the suppression of the $`\nu _\tau `$ CC cross section relative to $`\sigma _{\mathrm{CC}}(\nu _\mu N)`$. In Eq. (6), the $`F_5`$ term appears with a minus sign, and no factor of $`x`$. Since $`F_5F_1q(x,Q^2)`$, there is a small-$`x`$ enhancement of its contribution to the cross section at high energies. The $`F_5`$ term accounts for the rest of the suppression of the $`\nu _\tau `$ cross section at high energies. The tau mass corrections to the prefactors of $`F_1,F_2`$ and $`F_3`$ become negligible at high energies because the low-$`x`$ rise of $`q(x)`$ is tempered by factors of $`x`$ or $`y`$ for these structure function.

Since mass effects for tau neutrino interactions persist to 1 TeV, it is interesting to compare to the case of muon neutrino CC interactions at low energies where the muon mass is important. Since the factor of lepton mass comes into the equations via $`m_\tau ^2/E_\nu `$, the energy for muon neutrino interactions equivalent to 1 TeV for $`\nu _\tau `$ is $`E_\nu =3.5`$ GeV. At this energy, $`\sigma _{\mathrm{CC}}(\nu _\mu N)`$ including $`m_\mu =0.106`$ GeV (without cuts) is 2% lower than the cross section with $`m_\mu =0`$. This smaller suppression is due to the fact that at low energies, one does not get significant small-$`x`$ contributions to $`F_5`$. Furthermore, for $`\nu _\mu N`$ scattering, this energy is in the range where the DIS approximation to the total cross section is not reliable.

Charm production in neutrino interactions is a small contribution. In Fig. 5, we show the total charged current cross section and the separate cross section for $`\nu NcX`$ as a function of incident neutrino energy. The dashed curves are for incident $`\nu _\mu `$, the solid curves for incident $`\nu _\tau `$. At $`E_\nu =100`$ GeV, charm production contributes about 7% of the cross section for both $`\nu _\mu `$ and $`\nu _\tau `$.

The K-factor, a comparison of the NLO to LO††††††The LO result is obtained by neglecting all $`𝒪\left(\alpha _s^1\right)`$ terms in Section II and employing the CTEQ6L PDFs. charged current cross sections for incident tau neutrinos is shown in Fig. 6. As to be expected, NLO corrections are most significant near threshold where virtualities $`Q^2`$ are lower and $`\alpha _s`$ is larger than at higher energies. At $`E_\nu =10`$ GeV, $`K=`$NLO/LO=1.12, reducing to $`K=1`$ at $`E50`$ GeV for the evaluation with no cuts. The three curves show that with $`W_{\mathrm{min}}>1.4`$ GeV and $`Q^2>2`$ GeV2, the K-factor is the same as for the uncut cross section except for $`E_\nu <10`$ GeV where the cuts improve the perturbative convergence. The K-factor with only the $`W_{\mathrm{min}}`$ cut is even slightly lower, confirming the expectation based on Fig. 4 that perturbation theory is well behaved and its convergence does not have to be improved by any further $`Q_{\mathrm{min}}`$ cut on top of $`W_{\mathrm{min}}`$. Overall, the $`K`$ factor is reasonably close to one to indicate trustworthy perturbative NLO predictions, even more so when the $`W_{\mathrm{min}}`$ cut is imposed.

### C Uncertainties

In this Section we investigate a few factors that cause theoretical errors in the evaluation of $`\sigma _{\mathrm{CC}}`$. A full assessment of uncertainties will have to be based on the combination of the quasi-elastic and resonant channels with our DIS results. At higher energies, though, where DIS becomes dominant, we already can provide a good guide towards error estimates.

An uncertainty in our evaluation of the cross section is due to the factorization scale dependence of the cross section at fixed order in perturbation theory. To evaluate the uncertainty due to the scale dependence, we have varied the factorization and renormalization scale ($`\mu =\mu _F=\mu _R`$) over a very wide range of $`\mu ^2=0.110Q^2`$ for energies between 5 GeV and 100 GeV. A selection of energies is shown in Fig. 7. We show the ratio of $`\sigma _{\mathrm{CC}}(\nu _\tau N)`$ as a function of $`\mu ^2/Q^2`$ to the cross section at $`\mu ^2=Q^2`$. The flat $`\mu `$ dependence is a reassuring feature of the NLO calculation as opposed to the monotonic decrease with $`\mu `$ observed in LO, see Fig. 8. While we discuss here the full range of $`\mu ^2`$ as plotted in Fig. 7, the perturbative stability observed in Fig. 6 suggests that one would very likely overestimate the uncertainty from higher orders by such a wide scale variation. This seems even more so if we look at the scale dependent K-factor in Fig. 9 which prefers scales close to the canonical DIS choice $`\mu =Q`$ where $`K1`$. Still, scale choices are arbitrary to a large extent and we prefer to discuss the full picture instead of narrowing it down to some window for $`\mu `$ around $`\mu =Q`$. At $`E_\nu =5`$ GeV, the variation is the largest for the cross section evaluated with no cuts, ranging between 1.16 and 0.66. As explained in Section III A, for $`\mu ^2<\mu _{\mathrm{cut}}^2`$, we always set $`\mu =\mu _{\mathrm{cut}}`$‡‡‡‡‡‡ I.e. strictly the plot should be labeled $`\sigma _{\mathrm{CC}}\left(\mathrm{max}\{\mu ,\mu _{\mathrm{cut}}\}\right)/\sigma _{\mathrm{CC}}\left(\mathrm{max}\{Q,\mu _{\mathrm{cut}}\}\right)`$. . For small $`\mu /Q`$, the ratio without cuts in $`\{W,Q\}`$ is nearly constant as a function of $`\mu ^2/Q^2`$, showing the degree to which the uncut cross section comes from small $`\mu `$ values where $`\mu _{\mathrm{cut}}`$ takes over and where the perturbative treatment is unreliable. With $`\{W,Q\}`$ cuts applied, the variation over the range of scales is slightly less. The scale uncertainty at low energies underlies a larger uncertainty associated with the DIS approximation near threshold where the quasi-elastic and resonant contributions are significant.

At higher energies, the variation in the ratio is between $`0.851`$. For $`E_\nu =20`$ GeV, the ratio in Fig. 7 is between 0.95 and 1 for $`\mu ^2/Q^20.24`$.

The uncertainty in the cross section due to uncertainties in the parametrization of the parton distribution functions (PDF) is harder to quantify than the scale dependence. One approach is to use different parton distribution functions than the default CTEQ6 set. As a comparison, the GRV98 PDF set, at $`E_\nu =20`$ GeV with $`\mu =Q`$, yields a cross section only about one percent smaller than the CTEQ6 set. At $`E_\nu =5`$ and 10 GeV, the GRV98 cross sections are 10% and 4% lower than the CTEQ6 cross sections when no cuts are applied. The larger deviation at lower energies is due to the different high-$`x`$ distributions in the two PDF sets. Including $`W>1.4`$ GeV results in a deviation of 3.6% at $`E_\nu =5`$ GeV, even less for higher energies. Table III represents results for a few selected energies with a cut in $`W`$ applied. As commented above in Section III A, the PDFs also differ in their treatment of the number of active flavors and in the value of the strong coupling constant.

The CTEQ collaboration has also provided distributions in addition to their best fit set . The 20 dimensional parameter space to which the PDFs are sensitive yields 40 PDF sets with plus/minus variations on the eigenvector directions in that space. The resulting error estimate on $`\sigma _{\mathrm{CC}}(\nu _\tau N)`$ from evaluating the 40 sets is 3% at $`E_\nu =20`$ GeV, see Table III for other values of $`E_\nu `$. Overall, the GRV98 results lie within the uncertainty estimate suggested by the CTEQ6 eigenvector PDFs. We can, therefore, be confident of the absence of a systematic effect from the number of flavours and that the statistical uncertainties as encoded in the CTEQ6 sets provide a realistic PDF error estimate for $`\sigma _{\mathrm{CC}}(\nu _\tau N)`$.

We have incorporated kinematic corrections due to including the target hadron mass $`M_N`$ by employing the parton light cone fraction $`\xi `$, which equals the Nachtmann variable $`\eta `$  at leading order in the massless quark limit. One finds that $`\eta `$ is much different than Bjorken-$`x`$ at large $`x`$. For example, for $`Q^2=2(10)`$ GeV2, $`\eta =0.45(0.49)`$ at $`x=0.5`$ and $`\eta =0.75(0.92)`$ at $`x=1`$. The use of $`\eta `$ rather than $`x`$ in the structure functions has the largest impact at high $`x`$ and low $`Q^2`$.

Target mass corrections are also included via Eqs. (23)-(27), in which the $`_i`$ are mixed with target mass dependent prefactors for a given $`F_i`$. These formulae are based on the assignment of $`p_+=\xi P_+`$ to the light-cone momentum of the massless incident parton ($`p^2=0`$) given $`P^2=M_N^2`$. The parton and nucleon are assumed to have collinear momenta, $`p_{}=0`$. The formalism is discussed in detail in Ref. , however, other choices for the model of including target mass effects are possible, for example, including parton transverse momentum  $`p_{}0`$. The model of Ref.  reproduces the kinematic corrections for the leading twist operator product expansion result . As a consequence, our results here serve only as a guide to the magnitude of target mass corrections to the charged current cross sections. For the transverse structure function $`F_T`$ of Ref. , which equals $`_2`$ in the collinear parton approximation, the difference between the $`p_{}0`$ approximation giving $`F_T`$ and our order $`p_{}^0`$ approximation is less than 10% over all $`x`$ and $`Q^2`$ as low as 1 GeV2, rapidly falling below 5% at $`Q^2=`$2 GeV2. Low energy $`𝒪(p_{}^2)`$ effects and the target mass treatment in general border onto dynamical higher twists at large-$`x`$ . Those concepts will be investigated in more detail when future work will combine the DIS cross-sections with the non-DIS channels.

The effects of the target mass corrections as implemented here on $`F_4`$ and $`2xF_5F_2`$, as shown in Figs. 1 and 2, are small. In terms of the total cross section, the lowest energies are most affected. In the absence of kinematic cuts, $`\sigma _{\mathrm{CC}}(\nu _\tau N)(E_\nu =10\mathrm{GeV})`$ for $`M_N=0`$ is 8% larger than the cross section including the target mass via $`\eta `$ and Eqs. (23)-(27). For $`E>20`$ GeV, the deviation is less than 5%, down to the 2% level at 50 GeV incident neutrino energy. When one includes $`W>W_{\mathrm{min}}=1.4`$ GeV, the effect of target mass corrections on the $`\nu _\tau N`$ CC cross section is less than 2% for $`E_\nu >8`$ GeV. The reduced effect is due to the fact that the $`W_{\mathrm{min}}`$ value reduces the region of integration for large $`x`$ since $`2M_NE_\nu y(1x)W_{\mathrm{min}}^2M_N^2`$.

## IV Conclusions

The NLO corrections for the $`\nu _\tau N`$ charged current cross section have a relatively large impact at low neutrino energies near threshold and less of an impact at high energies. We have already shown the K-factor in Fig. 6. In Fig. 10, we show $`\sigma _{\mathrm{CC}}(\nu _\tau N)/E_\nu `$ versus neutrino energy at NLO including target mass corrections as implemented via Eqs. (23)-(27) and charm mass correction. This is compared to the naive evaluation, neglecting masses and NLO corrections, where the Albright-Jarlskog relations are correct. At high energies, the target mass corrections are negligible. The main effect is due to the charm mass threshold. At lower energies, the larger QCD K-factor is compensated by the reduction in the cross section due to target mass and charm mass effects.

In the theoretical evaluation of the $`\nu _\tau N`$ charged current cross section, the scale dependence in the PDFs and $`\alpha _s`$ is a large uncertainty at low energies. The parametrization of the PDFs does not seem to be a large uncertainty in the evaluation of the total cross section, especially when one applies DIS cuts on $`W^2`$ and $`Q^2`$.

Target mass corrections and the importance of $`Q^2<1`$ GeV2 are small at high energies, but they are significant in the 10-20 GeV energy range and lower, especially for the cross section without kinematic cuts. Our implementation of target mass corrections neglects parton transverse momentum ($`p_{}^2`$). This approach makes the inclusion of NLO QCD corrections straightforward, however, it neglects corrections of order $`M_N^2/Q^2`$ induced by nonzero $`p_{}^2`$ . As a consequence, one should view our target mass corrections as approximate.

In any case, the DIS cross section in the energy region below $`E_\nu 1020`$ GeV is difficult to interpret. At best, the uncut cross section is a crude approximation to the total cross section including resonant hadron production, e.g., $`\mathrm{\Delta }`$ production. The cross section with cuts is likely a better representation of the non-resonant neutrino-nucleon interactions, however, the issue of avoiding double-counting in combining resonant and non-resonant interactions is not solved theoretically. Phenomenological approaches are being explored  for applications to $`\nu _\mu N`$ interactions in the few GeV region that may also be applied to $`\nu _\tau N`$ interactions at higher energies. The universality of $`M_N^2/Q^2`$ corrections, carried over from copious electromagnetic interaction data to the context of weak interactions, is not completely clear.

The $`\nu N`$ cross section is an important ingredient in current and future atmospheric and neutrino factory experiments. Our evaluation of the NLO corrections for $`\nu _\tau `$N CC interactions including charm mass corrections and an estimate of target mass effects is part of a larger theoretical program to understand the inelastic $`\nu N`$ cross section over the full energy range relevant to experiments.

###### Acknowledgements.

We thank K. Ellis, M. Goodman, R. Jakob, F. Olness, W.K. Tung and U.K. Yang for discussions. This research was supported in part by the National Science Foundation under Grants PHY-0070443 and PHY-9802403 and DOE Contract No. FG02-91ER40664.

## Appendix A: The Albright-Jarlskog Relations

In this Appendix we pinpoint the approximations for which the Albright-Jarlskog relations (AJRs) (1) and (2) hold beyond the naive parton model (where they are exact). Below, we will rewrite the AJRs in terms of helicity amplitudes. First, however, we can immediately tell from Eqs. (23)–(27) that the mixing of $`\{_2^c,_4^c,_5^c\}`$ for $`M_N0`$ will violate both (1) and (2). We will, therefore, restrict the following discussion to $`M_N=0`$. The charm mass will be retained to trace its (somewhat less obvious) impact on the AJRs. In terms of helicity projections

$`W_0`$ $`=`$ $`\epsilon _0^\mu \epsilon _0^\nu W_{\mu \nu }`$ (A1)
$`W_s`$ $`=`$ $`\left(\epsilon _0^\mu \epsilon _q^\nu +\epsilon _q^\mu \epsilon _0^\nu \right)W_{\mu \nu }`$ (A2)
$`W_q`$ $`=`$ $`\epsilon _q^\mu \epsilon _q^\nu W_{\mu \nu }`$ (A3)
$`W_\pm `$ $`=`$ $`\epsilon _\pm ^\mu \epsilon _{\pm }^{\nu }{}_{}{}^{}W_{\mu \nu }`$ (A4)
with polarization vectors in terms of virtual W momentum $`q`$ and longitudinal reference vector $`k`$

$`\epsilon _\pm ^\mu `$ $`={\displaystyle \frac{1}{\sqrt{2}}}(0,1,i,0)`$ $`(\mathrm{𝑡𝑟𝑎𝑛𝑠𝑣𝑒𝑟𝑠𝑒})`$ (A5)
$`\epsilon _q^\mu `$ $`={\displaystyle \frac{q^\mu }{\sqrt{q^2}}}`$ $`(\mathrm{𝑠𝑐𝑎𝑙𝑎𝑟})`$ (A6)
$`\epsilon _0^\mu `$ $`={\displaystyle \frac{(q^2)k^\mu +(kq)q^\mu }{\sqrt{(q^2)[(kq)^2]}}}`$ $`(\mathrm{𝑙𝑜𝑛𝑔𝑖𝑡𝑢𝑑𝑖𝑛𝑎𝑙})`$ (A7)
the tensor basis $`_{i=1,\mathrm{},5}^c`$ can be written as:

$`_1^c`$ $`={\displaystyle \frac{1}{2}}(W_++W_{})`$ (A8)
$`_2^c`$ $`={\displaystyle \frac{\lambda }{2}}(W_++W_{}+2W_0)`$ (A9)
$`_3^c`$ $`={\displaystyle \frac{1}{2}}(W_+W_{})`$ (A10)
$`_4^c`$ $`={\displaystyle \frac{1}{2}}(W_0+W_qW_s)`$ (A11)
$`_5^c`$ $`={\displaystyle \frac{1}{2}}(W_++W_{}+2W_0W_s).`$ (A12)
We then see that

$$F_2^c2xF_5^c=2\frac{x}{\lambda }_2^c2x_5^c=xW_s$$
(A13)
singles out the interference term $`W_s`$ between scalar and longitudinal polarization. The latter involves a contraction with $`\epsilon _q^\mu q^\mu `$. Now, as long as we make the single $`W`$ boson exchange approximation, the DIS process is equivalent to an effectively abelian (electroweak) interaction if all quarks are massless (or of the same mass). Then, $`F_22xF_5=0`$ is guaranteed by naive gauge invariance under $`\epsilon _q^\mu \epsilon _q^\mu q^\mu /\sqrt{q^2}=0^\mu `$. This is a stronger statement than helicity conservation for massless spin-1/2 quarks because gauge invariance holds to any order in $`\alpha _s`$ while helicity conservation breaks down when non-collinear NLO radiation generates angular momentum. A transition $`W^+sc`$ with $`0=m_sm_c`$, however, is necessarily non-diagonal in flavour space and the interaction cannot be abelianized. Thus, naive gauge invariance does not hold and the second AJR (2) is not protected from charm mass corrections. It will, however, hold at any order in $`\alpha _s`$ for massless quarks.

Worth mentioning is also that $`F_4^c=0`$ at $`𝒪\left(\alpha _s^0\right)`$ even for the charm production process, indicating that $`W_0+W_qW_s=0`$ at LO. This should be compared to the longitudinal structure function $`F_22xF_1W_0`$ which does not vanish by helicity conservation in LO when a massive charm quark is produced in the final state. Now, $`F_4`$ is obtained from

$$W_0+W_qW_s=\frac{k^\mu k^\nu }{Q^2}W_{\mu \nu }.$$
(A14)
At the parton level, $`k`$ is the incoming parton momentum $`p`$. $`F_4`$ does, therefore, not receive corrections from the final state charm mass in LO as long as initial state down and strange masses are zero ($`p^2=0`$).