text:
# Stochastic semiclassical gravity and fluctuations during inflation

## I Introduction

Let us briefly summarize the road to stochastic semiclassical gravity. It starts with quantum field theory in a curved spacetime which is now a well understood and well defined theory at least for free fields . In this theory the gravitational field is the classical field of general relativity, that is the metric of the spacetime, and the quantum fields propagate in such a spacetime. Since the spacetime is now dynamical it is not always possible to define a physically meaningful vacuum state for the quantum field and when this is possible in some “initial” times it is usually unstable, in the sense that it may differ from the vacuum state at latter times, and spontaneous creation of particles occurs. Applications of this in cosmology, such as particle production in expanding Friedmann-Robertson-Walker models , and black hole physics, such as Hawking radiation , are well known. This is one aspect of the interaction of gravity with quantum matter fields

Another aspect of this interaction is the back-reaction of the quantum fields on the spacetime. Since the gravitational field couples to the stress tensor of matter fields, the key object here is the expectation value in a given quantum state of the stress energy tensor of the quantum field, which is a classical observable. However, since this object is quadratic in the field operator, which is only well defined as a distribution on spacetime, it involves ill defined quantities which translate into ultraviolet divergencies. To be able to define a physically meaningful quantity a regularization and a renormalization procedure is required. The ultraviolet divergencies associated to the expectation value of the stress energy tensor are also present in Minkowski spacetime, but in a curved background the renormalization procedure is more sophisticated as it needs to preserve general covariance. A regularization procedure which is specially adapted to the curved background is the so called point-splitting method . The final expectation value of the stress energy tensor using point splitting or any other reasonable regularization technique is essentially unique, modulo some terms which depend on the spacetime curvature and which are independent of the quantum state. This uniqueness was essentially proved by Wald  who investigated the criteria that a physically meaningful expectation value of the stress energy tensor ought to satisfy.

Now the back-reaction problem may be formulated in terms of the so called semiclassical Einstein equations. These are Einstein equations which have the expectation value in some quantum state of the stress energy tensor as a matter source. The back-reaction problem was investigated in cosmology, in particular to see whether cosmological anisotropies could be damped by back-reaction . This was an earlier attempt  previous to inflation to explain why the universe is so isotropic at present.

Another step forward in the back-reaction problem was the use of effective action methods  so familiar in quantum field theory. These methods were of great help in the study of cosmological anisotropies since they allowed the introduction of familiar perturbative treatments into the subject. The most common effective action method, however, led to equations of motion which were not real because they were taylored to compute transition elements of quantum operators rather than expectation values. Fortunately the appropriate technique had already been developped by Schwinger and Keldysh  in the so called Closed Time Path (CTP) or in-in effective action method, and was soon adapted to the gravitational context . These techniques were then applied to different problems in the cosmological context incuding the effects of arbitrary perturbations on homogeneous backgrounds . As a result one was now able to deduce the semiclassical Einstein equations by the CTP functional method: starting with an action for the interaction of gravity with matter fields, treating the matter fields as quantum fields and the gravitational field at tree level only.

The semiclassical Einstein equations have limitations, it is clear that even outside Planck scales if fluctuations on the expectation value of matter fields are large, and that depends on the quantum state, the semiclassical equations should break down . One expects, in fact, that a better approximation would describe the gravitational field in a probabilistic way. In other words, that the semiclassical equations should be substituted by some Langevin-type equations with a stochastic source that describes the quantum fluctuations. A significant step in this direction was made by Hu  who proposed to view the back-reaction problem in the framework of an open quantum system, where the quantum fields are seen as the “environment” and the gravitational field is seen as the “system”. Following this proposal a systematic study of the connection between semiclassical gravity and open quantum systems resulted in the development of a framework were semiclassical Einstein-Langevin equations could be derived . The key technical factor to most of these results was the use of the influence functional method of Feynman and Vernon  for the description of the system-environment interaction when only the state of the system is of interest. The CTP method for open systems involves, in fact, the influence functional.

However although several Einstein-Langevin equations were derived, the results were somewhat formal and some concern could be raised on the physical reality of the solutions of the stochastic equations for the gravitational field. This is related to the issue of the environment induced quantum to classical transition. In fact, for the existence of a semiclassical regime for the dynamics of the system one needs two requeriments, in the language of the consistent histories formulation of quantum mechanics . The first is decoherence, which guarantees that probabilities can be consistently assigned to histories describing the evolution of the system, and the second is that these probabilities should peak near histories which correspond to solutions of classical equations of motion. The effect of the environment is crucial on the one hand to provide decoherence  and on the other hand to produce both dissipation and noise to the system through back-reaction, thus inducing a semiclassical stochastic dynamics on the system. As shown by Gell-Mann and Hartle  in an open quantum system stochastic semiclassical equations are obtained after a coarse graining of the environmental degrees of freedom and a further coarse graining in the system variables. That this mechanism could also work for decoherence and classicalization of the metric field was not so clear lacking a quantum description of the gravitational field, and the analogy could be made only formaly .

Thus, an axiomatic approach to the Einstein-Langevin equations which was independent of the open system analogy was suggested: it was based on the formulation of consistent dynamical equations for a perturbative correction to semiclassical gravity able to account for the lowest order stress energy fluctuations of matter fields . It was later shown that these same equations could be derived, in this general, case from the influence functional of Feynman and Vernon in which, the gravitational field is treated at tree level and and quantum fields are quantized, the first being, in fact, the “system” and the seconds the “environment” .

Here we review some of these developments in semiclassical stochastic gravity and also some of its applications. In section II a brief sketch of semiclassical gravity is given. In section III the axiomatic approach to the Einstein-Langevin equations is discussed. To illustrate the relation between the semiclassical, stochastic semiclassical and quantum theories, a simplified model of linear gravity is used. Throughout this section we use a simplified notation, avoiding tensorial indices when possible and emphasizing the conceptual aspects. In section IV an important application of stochastic gravity is discussed in some detail, it concerns the computation of the two-point correlations of the metric perturbations induced by the fluctuation in the stress energy tensor of the inflaton field during inflation. The results agree with the standard results but the present method in which the gravitational field and the matter fields are treated separately may have some advantages over the methods were both the metric and matter perturbations are treated on the same footing from the start. Finally, in section V we summarize our result and briefly discuss other applications. We should mention that an extensive and stimulating review of stochastic semiclassical gravity is given in ref. , it includes the suggestive idea of considering classical general relativity it terms of collective variables and of viewing semiclassical gravity as mesocopic physics.

## II Semiclassical gravity

Semiclassical gravity is a theory which describes the interaction of the gravitational field assumed to be a classical field with matter fields which are quantum. It is supposed to be some limit of the still unknown theory describing the interaction of quatum gravity with quantum fields. Due to the lack of the quantum theory, the semiclassical limit cannot be rigurously derived. However, it can be formally derived in several ways. One of them is the leading-order $`1/N`$ approximation of quantum gravity , where $`N`$ is the number of independent free quantum fields which interact with gravity only, and where one keeps finite the value of $`NG`$, where $`G`$ is Newton’s gravitational constant. In this limit, after path integration one arrives at a theory in which formally the gravitational field can be treated as a c-number (i.e. is quantized at tree level) and the quantum fields are fully quantized. If we call $`g`$ the metric tensor and $`\widehat{\varphi }`$ the scalar field (for simplicity we consider just one scalar field) one arrives at the semiclassical Einstein equation as the dynamical equation for the metric $`g`$:

$$G_g=8\pi G\widehat{T}^R_g,$$
(1)
where $`\widehat{T}=T[\widehat{\varphi }^2]`$ is the stress energy tensor in a simplifyied notation, which is quadratic in the field operator $`\widehat{\varphi }`$. This operator, being the product of distribution valued operators, is ill defined and needs to be regularized and renormalized, the $`R`$ in $`\widehat{T}^R`$ means that the operator has been renormalized. The angle brackets on the right hand side mean that the expectation value of the stress tensor operator is computed in some quantum state, say $`|\psi `$, compatible with the geometry described by the metric $`g`$. On the left hand side $`G_g`$ stands for the Einstein tensor for the metric $`g`$ together with the cosmological constant term and other terms quadratic in the curvature which are generally needed to renormalize the stress energy tensor operator. The quantum field operator $`\widehat{\varphi }`$ propagates in the background defined by the metric $`g`$, it thus satisfies a Klein-Gordon equation, which we write also schematically as

$$\mathrm{}_g\widehat{\varphi }=0,$$
(2)
where $`\mathrm{}_g`$ stands for the D’Alambert operator in the background of $`g`$. Equation (1) is the semiclassical Einstein equation, it is the dynamical equation for the metric tensor $`g`$ and describes the back-reaction of the quantum matter fields on the geometry. A solution of semiclassical gravity consists of the set $`(g,\widehat{\varphi },|\psi )`$ where $`g`$ is a solution of (1), $`\widehat{\varphi }`$ is a solution of (2) and $`|\psi `$ is the quantum state in which the expectation value of the stress energy tensor in (1) is computed.

For a free quantum field this thery is robust in the sense that it is consistent and fairly well understood. Note that it is in some sense unique as a theory where the gravitational field is classical. In fact, the (classical) gravitational field interacts with matter fields through the stress energy tensor, and the only reasonable c-number stress energy tensor that one may construct  with the operator $`\widehat{T}`$ is just the right hand side of (1). However the scope and limits of the theory are not so well understood as a consequence of the lack of the full quantum theory. It is assumed that the semiclassical theory should break down at Planck scales, which is when simple order of magnitude estimes suggest that the quantum effects of gravity cannot be ignored: the gravitational energy of a quantum fluctuation of energy in a Planck size region, determined by the Heisenberg uncertainty principle, are of the same order of magnitude.

There is also another situation when the semiclassical theory should break down, namely, when the fluctuations of the stress energy tensor are large. This has been emphasized by Ford and collaborators. It is less clear how to quatify what a large fluctuation here means and some criteria have been proposed . Generally this depends on the quantum state and may be illustrated by the example used in ref.  as follows.

Let us assume a quantum state formed by an isolated system which consists of a superposition with equal amplitude of one configuration with mass $`M_1`$ and another with mass $`M_2`$. Semiclassical theory as described in (1) predicts that the gravitational field of this system is produced by the average mass $`(M_1+M_2)/2`$, that is a test particle will move on the background spacetime produced by such a source. However one would expect that if we send a succesion of test particles to probe the gravitational field of the above system half of the time they would react to the field of a mass $`M_1`$ and the other half to the field of a mass $`M_2`$. If the two masses differ substantially the two predictions are clearly different, note that the fluctuations in mass of the quantum state is of the order of $`(M_1M_2)^2`$. Although the previous example is suggestive a word of caution should be said in order not to take it too literaly. In fact, if the previous masses are macroscopic the quantum system decoheres very quickly  and instead of a pure quantum state it is described by a density matrix which diagonalizes in a certain pointer basis. Thus for observables associated to this pointer basis the matrix density description is equivalent to that provided by a statistical ensemble. In any case, however, from the point of view of the test particles the predictions differ from that of the semiclassical theory.

## III Einstein-Langevin equation

The purpose of semiclassical stochastic gravity is to be able to deal with the situation of the previous example in which the predictions of the semiclassical theory may be inacurate. Consequently, our first point is to characterize the quantum fluctuations of the stress energy tensor.

The physical observable that measures these fluctuations is $`\widehat{T}^2\widehat{T}^2`$. To make this more precise let us introduce the tensor operator $`\widehat{t}=\widehat{T}\widehat{T}\widehat{I}`$, where $`\widehat{I}`$ is the identity operator, then we introduce the noise kernel as the four index bi-tensor defined as the expectation value of the anticommutator of the operator $`\widehat{t}`$:

$$N(x,y)=\frac{1}{2}\{\widehat{t}(x),\widehat{t}(y)\}_g.$$
(3)
The subindex $`g`$ here means that this expectation value in taken in the background metric which is a solution of the semiclassical equation (1). An important property of the symmetric bi-tensor $`N(x,y)`$ is that it is finite because the tensor operator $`\widehat{t}`$ is finite since the ultraviolet divergencies of $`\widehat{T}`$ are cancelled by the substraction of $`\widehat{T}`$. Since the operator $`\widehat{T}`$ is selfadjoint $`N(x,y)`$, which is the expectation value of an anticommutator, is real and positive semi-definite. This last property allows for the introduction of a classical Gaussian stochastic tensor $`\xi `$ defined by

$$\xi (x)_c=0,\xi (x)\xi (y)_c=N(x,y).$$
(4)
This stochastic tensor is symmetric $`\xi _{\mu \nu }=\xi _{\nu \mu }`$ and divergenceless, $`^\mu \xi _{\mu \nu }=0`$, as a consequence of the fact that the stress tensor operator is divergenceless. The subindex $`c`$ means that the expectation value is just a classical average. Note that we assume that $`\xi `$ is Gaussian just for simplicity in order to include the main effect. The idea now is simple we want to modify the semiclassical Einstein equation (1) by introducing a linear correction to the metric tensor $`g`$, such as $`g+h`$, which accounts consistently for the fluctuations of the stress energy tensor. The simplest equation is,

$$G_{g+h}=8\pi G(\widehat{T}^R_{g+h}+\xi ),$$
(5)
where $`g`$ is assumed to be a solution of equation (1). This stochastic equation must be thought of as a linear equation for the metric perturbation $`h`$ which will behave consequently as a stochastic field tensor. Note that the tensor $`\xi `$ is not a dynamical source, since it has been defined in the background metric $`g`$ which is solution of the semiclassical equation. Note also that this source is divergenceless with respect to the metric, and it is thus consistent to write it on the right hand side of the Einstein equation. This equation is gauge invariant with respect to diffeomorphisms defined by any field on the background spacetime . If we take the statistical average equation (5) becomes just the semiclassical equation for the metric $`g+h`$ where now the expectation value of $`\widehat{T}`$ is taken in the perturbed spacetime.

The stochastic equation (5) is known as the Einstein-Langevin equation. It is linear in $`h`$, thus its solutions can be written as the sum of a solution of the homogeneous equation plus a stochastic part: $`h=h^h+h^s`$. The equation predicts that the gravitational field has stochastic fluctuations over the background $`g`$. The correlation function for the gravitational field is simply given by $`h^s(x)h^s(y)_c`$. This is the physically most relevant observable, to find it requires to solve the Einsten-Langevin equation and to know the noise kernel $`N(x,y)`$. Note that the noise kernel should be thought of as a distribution function, the limit of coincidence points has meaning only in the sense of distributions. Explicit expressions of this kernel in terms of the two point Wightman functions is given in , expression based on point-splitting methods have also been given in .

This stochastic theory goes beyond semiclassical gravity in the following sense. The semiclassical theory, which is based on the expectation value of the stress energy tensor, carries information on the field two-point correlations only, since $`\widehat{T}`$ is quadratic in the field operator $`\widehat{\varphi }`$. The stochastic semiclassical theory on the other hand, is based on the noise kernel (3) which is quartic in the field. It is thus clear that it carries information beyond the semiclassical theory. In some sense the theory represents a middle road between the semiclassical and the quantum theory. It does not carry information, however, on the graviton-graviton interaction.

### A A simplified model

To clarify this point it is useful to introduce a toy model which would describe exactly a linear theory such as electromagnetism, and captures some essential features of linearized gravity. Let us assume that the gravitational equations are described by the linear equations for the field $`h`$, with a source $`T[\varphi ^2]`$. The semiclassical equations now correspond to

$$\mathrm{}h=\widehat{T},$$
(6)
where now $`\widehat{T}`$ depends on the field operator $`\widehat{\varphi }`$. Note that this theory does not correspond to linearized semiclassical gravity around the Minkowski background since in that case $`\widehat{T}=0`$ assuming that $`T`$ does not depend on $`h`$ (in the linearized semiclassical theory $`\widehat{T}`$ is, in fact, linear $`h`$ ). The model, however, can be extended to linearized gravity .

The solutions of this equation may be written in terms of the retarded propagator $`G_{xy}`$ as

$$h_x=h_x^0+G_{xx^{}}\widehat{T}_x^{},$$
(7)
where $`h_x^0`$ is the homogeneous solution which is determined by the initial conditions.

Let us now consider the quantum theory which, in Heisenberg representation, may be written as

$$\mathrm{}\widehat{h}=\widehat{T}.$$
(8)
The solutions of this equation may be written as

$$\widehat{h}_x=\widehat{h}_x^0+G_{xx^{}}\widehat{T}_x^{},$$
(9)
and one may compute the two point quatum correlation function as

$$\widehat{h}_x\widehat{h}_y=\widehat{h}_x^0\widehat{h}_y^0+G_{xx^{}}G_{yy^{}}\widehat{T}_x^{}\widehat{T}_y^{},$$
(10)
where the expectation value is taken in the quantum state in which both fields $`\varphi `$ and $`h`$ have been quantized, and we have used that for the free field $`\widehat{h}^0=0`$.

At this point we may now compare with the stochastic theory as described by equation (5) which, in this simplified model, may be written as

$$\mathrm{}h=\widehat{T}+\xi ,$$
(11)
where $`\xi `$ is a Gaussian stochastic source defined by $`\xi _x_c=0`$ and $`\xi _x\xi _y_c=\widehat{T}_x\widehat{T}_y\widehat{T}_x\widehat{T}_y`$, where the expectation values on the right hand side are defined in a given state of the field $`\widehat{\varphi }`$ and the subscript $`c`$ on the left hand sides means statistical average, see (3) and (4). Now the solution of this equation may be written in terms of the retarded propagator as,

$$h_x=h_x^0+G_{xx^{}}\left(\widehat{T}_x^{}+\xi _x^{}\right),$$
(12)
from where the two point correlation function for the classical field $`h`$, after using the definition of $`\xi `$ and that $`h^0_c=0`$, is given by

$$h_xh_y_c=h_x^0h_y^0_c+G_{xx^{}}G_{yy^{}}\widehat{T}_x^{}\widehat{T}_y^{}.$$
(13)
Comparing (10) with (13) we see that the respective second terms on the right hand side are identical provided the expectation values are computed in the same quantum state for the field $`\widehat{\varphi }`$, note that we have assumed that $`T`$ does not depend on $`h`$. The fact that the field $`h`$ is also quantized in (10) does not change the previous statement. In the real theory of gravity $`T`$, in fact, depends also on $`h`$ and then the previous statement is only true approximately, i.e perturbatively in $`h`$. The nature of the first terms on the right hand sides of equations (10) and (13) is different: in the first case it is the two point expectation value of the free quatum field $`\widehat{h}_0`$ whereas in the second case it is the average of the two point classical average of the homogeneous field $`h_0`$, which depends on the initial contions. Now we can still make these terms to be equal to each other if we assume for the homogeneous field $`h`$ a distribution of initial conditions such that $`h_x^0h_y^0_c=\widehat{h}_x^0\widehat{h}_y^0`$. Thus, under this assumption on initial conditions for the field $`h`$ the two point correlation function of (13) equal the quantum expectation value of (10) exactly. Thus in a linear theory as in the model just described one may just use the statistical description given by (11) to compute the quantum two point function of equation (10). Of course, the statistical description is not able to account for graviton-graviton effects which go beyond the linear approximation in $`\widehat{h}`$.

### B Functional approach

To end this section we should mention that the Einstein-Langevin equation (5) may also be formaly derived using the CTP functional method . As remarked in the introduction the CTP functional was introduced by Schwinger  to compute expectation values. One just considers the interaction of the gravitational field $`g`$ at tree level and of the quantum field $`\varphi `$ fully quantum. Then the effective action for the gravitational field is derived after integrating out the degrees of freedom of the quantum field, and the CTP influence action reduces basically to the Feynman and Vernon influence functional  used in quantum open systems. Here the system is the gravitational field and the environment is the quantum field. The stochastic terms for the gravitational field are found by suitably interpreting some pure imaginary term which appear in the influence action. These terms are closely connected to Gell-Mann and Hartle decoherence functional  used to study decoherence and classicalization in open quantum systems. The net result of these analogies is that the interaction with the environmennt induces fluctuations in the system dynamics. It is precisely the noise kernel introduced in (3) that accounts for this effect.

## IV Gravitational fluctuations during inflation

An important application of stochastic semiclassical gravity is the derivation of the cosmological perturbations generated during inflation . Let us consider the Lagrangian density for the inflaton field $`\varphi `$ of mass $`m`$

$$(\varphi )=\frac{1}{2}(\varphi )^2+\frac{1}{2}m^2\varphi ^2,$$
(14)
which is the basis of the simplest chaotic inflationary model . The conditions for the existence of an inflationary period, which is characterized by an accelerated expansion of the spacetime, is that the value of the field averaged over a region with the typical size of the Hubble radius is higher than the Planck mass $`m_P`$. This is because in order to solve the cosmological horizon and flatness problem more than 60 e-folds of expansion are needed, to achieve this the scalar field should begin with a value higher than $`3m_P`$. Furthermore, as we will see, the large scale anisotropies measured  restrict the inflaton mass to be of the order of $`10^6m_P`$.

We want to study the metric perturbations produced by the stress tensor fluctuations of the inflaton field on the homogeneous background of a flat Friedmann-Robertson-Walker model, described by the cosmological scale factor $`a(\eta )`$, where $`\eta `$ is the conformal time, which is driven by the homogenous inflaton field $`\varphi (\eta )=\widehat{\varphi }`$. Thus we write the inflaton field in the following form

$$\widehat{\varphi }=\varphi (\eta )+\widehat{\phi }(x),$$
(15)
where $`\widehat{\phi }(x)`$ corresponds to a free massive quantum scalar field with zero expectation value on the homogeneous background metric: $`\widehat{\phi }_g=0`$. Restrincting ourselves to scalar-type perturbations the perturbed metric $`\stackrel{~}{g}=g+h`$ can be written in the longitudinal gauge as,

$$ds^2=a^2(\eta )[(1+2\mathrm{\Phi }(x))d\eta ^2+(12\mathrm{\Psi }(x))\delta _{ij}dx^idx^j],$$
(16)
where the metric perturbations $`\mathrm{\Phi }(x)`$ and $`\mathrm{\Psi }(x)`$ correspond to Bardeen’s gauge invariant variables .

The Einstein-Langevin equation as described in the previous section is gauge invariant, and thus we can work in a desired gauge and then extract the gauge invariant quantities. The Einstein-Langevin equation (5) reads now:

$$G^{(0)}8\pi G\widehat{T}^{(0)}_g+G^{(1)}(h)8\pi G\widehat{T}^{(1)}(h)_g=8\pi G\xi ,$$
(17)
where the two first terms cancel, that is $`G^{(0)}8\pi G\widehat{T}^{(0)}_g=0`$, as the background metric satisfies the semiclassical Einstein equations. Here the subscripts $`(0)`$ and $`(1)`$ refer to functions in the background metric $`g`$ and linear in the metric perturbation $`h`$, respectively. The stress tensor operator $`\widehat{T}`$ for the minimally coupled inflaton field in the perturbed metric $`\stackrel{~}{g}=g+h`$ is:

$$\widehat{T}_{\mu \nu }=_\mu \widehat{\varphi }_\nu \widehat{\varphi }+\frac{1}{2}\stackrel{~}{g}_{\mu \nu }(_\rho \widehat{\varphi }^\rho \widehat{\varphi }+m^2\widehat{\varphi }^2).$$
(18)
Now using the decomposition of the scalar field into its homogeneous and inhomogeneous part, see (15), and the metric $`\stackrel{~}{g}`$ into its homogenous background $`g`$ and its perturbation $`h`$, the renormalized expectation value for the stress tensor can be written as

$$\widehat{T}^R[\stackrel{~}{g}]=\widehat{T}[\stackrel{~}{g}]_{\varphi \varphi }+\widehat{T}[\stackrel{~}{g}]_{\varphi \phi }+\widehat{T}^R[\stackrel{~}{g}]_{\phi \phi },$$
(19)
where only the homogeneous solution for the scalar field contributes to the first term. The second term is proportional to $`\widehat{\phi }[\stackrel{~}{g}]`$ which is not zero because the field dynamics is considered on the perturbed spacetime, i.e. this term includes the coupling of the field with $`h`$. The last term corresponds to the expectation value to the stress tensor for a free scalar field on the spacetime of the perturbed metric.

We can now compute the noise kernel $`N(x,y)`$ defined in equation (3), which after using the previous decomposition can be written as

$$\{\widehat{t},\widehat{t}\}[g]=\{\widehat{t},\widehat{t}\}_{\varphi \phi }[g]+\{\widehat{t},\widehat{t}\}_{\phi \phi }[g],$$
(20)
where we have used the fact that $`\widehat{\phi }_g=0=\widehat{\phi }\widehat{\phi }\widehat{\phi }_g`$ for Gaussian states on the background geometry. We have considered the vacuum state to be the Euclidean vacuum which is preferred in the de Sitter background, and this state is Gaussian. In the above equation the first term is quadratic in $`\widehat{\phi }`$ whereas the second one is quartic, both contributions to the noise kernel are separately conserved since both $`\varphi (\eta )`$ and $`\widehat{\phi }`$ satisfy the Klein-Gordon field equations on the background spacetime. Consequently, the two terms can be considered separately. On the other hand if one treats $`\widehat{\phi }`$ as a small perturbation the second term in (20) is of lower order than the first and may be neglected consistently, this corresponds to neglecting the last term of (19). The stress tensor fluctuations due to a term of that kind were considered in ref. .

We can now write down the Einstein-Langevin equations (17). It is easy to check that the space-space components coming from the stress tensor expectation value terms and the stochastic tensor are diagonal, i.e. $`\widehat{T}_{ij}=0=\xi _{ij}`$ for $`ij`$. This, in turn, implies that the two functions characterizing the scalar metric perturbations are equal: $`\mathrm{\Phi }=\mathrm{\Psi }`$ in agreement with ref. . The equation for $`\mathrm{\Phi }`$ can be obtained from the $`0i`$-component of the Einstein-Langevin equation, which in Fourier space reads

$$2ik_i(H\mathrm{\Phi }_k+\mathrm{\Phi }_k^{})=\frac{8\pi }{m_P^2}\xi _{k\mathrm{\hspace{0.17em}0}i},$$
(21)
where $`k_i`$ is the comoving momentum component associated to the comoving coordinate $`x^i`$. Here primes denote derivatives with respect to the conformal time $`\eta `$ and $`H=a^{}/a`$. A non-local term of dissipative character which comes from the second term in (19) should also appear on the left hand side of equation (21), but since we are mainly interested in the fluctuating part we have ignored this term. To solve this equation, whose left hand side comes from the linearized Einstein tensor for the perturbed metric , we need the retarded propagator for the gravitational potential $`\mathrm{\Phi }_k`$,

$$G_k(\eta ,\eta ^{})=i\frac{4\pi }{k_im_P^2}\left(\theta (\eta \eta ^{})\frac{a(\eta ^{})}{a(\eta )}+f(\eta ,\eta ^{})\right),$$
(22)
where $`f`$ is a homogeneous solution of (21) related to the initial conditions chosen. For intance, if we take $`f(\eta ,\eta ^{})=\theta (\eta _0\eta ^{})a(\eta ^{})/a(\eta )`$ the solution would correspont to “turning on” the stochastic source at $`\eta _0`$.

The correlation function for the metric perturbations is now given by

$$\mathrm{\Phi }_k(\eta )\mathrm{\Phi }_k^{}(\eta ^{})_c=(2\pi )^2\delta (\stackrel{}{k}+\stackrel{}{k}^{})^\eta 𝑑\eta _1^\eta ^{}𝑑\eta _2G_k(\eta ,\eta _1)G_k^{}(\eta ^{},\eta _2)\xi _{k0i}(\eta _1)\xi _{k^{}0i}(\eta _2)_c.$$
(23)
The correlation function for the stochastic source , which is connected to the stress tensor fluctuations through the noise kernel is given by,

$$\xi _{k0i}(\eta _1)\xi _{k0i}(\eta _2)_c=\frac{1}{2}\{\widehat{t}_{0i}^k(\eta _1,\widehat{t}_{0i}^k(\eta _2)\}_{\varphi \phi }=\frac{1}{2}k_ik_i\varphi ^{}(\eta _1)\varphi ^{}(\eta _2)G_k^{(1)}(\eta _1,\eta _2),$$
(24)
where $`G_k^{(1)}(\eta _1,\eta _2)=\{\widehat{\phi }_k(\eta _1),\widehat{\phi }_k(\eta _2)\}`$ is the $`k`$-mode Hadamard function for a free minimally coupled scalar field which is in the Euclidean vacuum on the de Sitter background.

It is useful to compute the Hadamard function for a massless field and consider a perturbative expansion in terms of the dimensionless parameter $`m/m_P`$. Thus we consider $`\overline{G}_k^{(1)}(\eta _1,\eta _2)=a(\eta _1)a(\eta _2)G_k^{(1)}(\eta _1,\eta _2)=0|\{\widehat{y}_k(\eta _1),\widehat{y}_k(\eta _2)\}|0=2\left(u_k(\eta _1)u_k^{}(\eta _2)\right)`$ with $`\widehat{y}_k(\eta )=a(\eta )\widehat{\phi }_k(\eta )=\widehat{a}_ku_k(\eta )+\widehat{a}_k^{}u_k^{}(\eta )`$ and where $`u_k=(2k)^{1/2}e^{ik\eta }(1i/\eta )`$ are the positive frequency $`k`$-mode for a massless minimally coupled scalar field on a de Sitter background, which define the Euclidean vacuum state: $`\widehat{a}_k|0=0`$ .

The background geometry, however, is not exactly that of de Sitter spacetime, for which $`a(\eta )=(H\eta )^1`$ with $`\mathrm{}<\eta <0`$. One can expand in terms of the “slow-roll” parameters and assume that to first order $`\dot{\varphi }(t)m_P^2(m/m_P)`$, where $`t`$ is the physical time. The correlation function for the metric perturbation (23) is the computed, see ref.  for details. The final result, however, is very weakly dependent on the initial conditions as one may understand from the fact that the accelerated expansion of de quasi-de Sitter spacetime during inflation erases the information about the initial conditions. Thus one may take the initial time to be $`\eta _0=\mathrm{}`$ and obtain to lowest order in $`m/m_P`$ the expression

$$\mathrm{\Phi }_k(\eta )\mathrm{\Phi }_k^{}(\eta ^{})_c8\pi ^2\left(\frac{m}{m_P}\right)^2k^3(2\pi )^3\delta (\stackrel{}{k}+\stackrel{}{k}^{})\mathrm{cos}k(\eta \eta ^{}).$$
(25)
From this result two main conclusions are derived. First, the prediction of an almost Harrison-Zel’dovich scale-invariant spectrum for large scales, i.e. small values of $`k`$. Second, since the correlation function is of order of $`(m/m_P)^2`$ a severe bound to the mass $`m`$ is imposed by the gravitational fluctuations derived from the small values of the Cosmic Microwave Background (CMB) anisotropies detected by COBE. This bound is of the order of $`(m/m_P)10^6`$ . One possible advantge of the Einstein-Langevin approach to the gravitational fluctuations in inflaton over the approach based on the quantization of the linear perturbations of both the metric and the inflaton field , is that an exact treatment of the inflaton quantum fluctuations is in principle possible, keeping the metric perturbations to linear order. On the other hand although the gravitational fluctuations are here assumed to be classical, the correlation functions obtained correspond to the quantum expectation values of the quantum metric perturbations , at least in the linear regime. This means that even in the absence of decoherence the fluctuations predicted by the Einstein-Langevin equation, whose solutions do not describe the actual dynamics of the gravitational field any longer, still give the correct quantum two-point functions.

## V Summary and outlook

We have reviewed the semiclassical theory of gravity as the theory of the interaction of classical gravity with quantum matter fields. The most important equations in this theory are the semiclassical Einstein equations (1) which describe the back-reaction of the gravitational fluctuations in its interaction with the quantum fields. We noticed that the theory may seriously fail when the fluctuations on the stress energy tensor of the quantum fields are important. We have then sought an axiomatic approach by which the semiclassical equations can be corrected in order to take into account those fluctuations. These equations turn out to be uniquely defined and are the Einstein-Langevin equations (5) which are linear in the metric perturbations $`h`$ over the semiclassical background. These equations predict stochastic fluctuations in the metric perturbations induced by the stress tensor fluctuations described by the noise kernel (3).

We have also noticed that the Einstein-Langevin equations can be formally derived from the open quantum system paradigm, in which the gravitational fluctuations $`h`$ and the quantum field interact when the interest is in the dynamics of the gravitational field. Thus treating the quantum field as the “environment” and the gravitational field as the “system”. The mathematical tools to carry out this approach are the CTP functional method and, in this context, its closely related Feynman and Vernon influence functional.

We have finally used the stochastic theory in the inflationary cosmological context. We have computed the two-point correlation functions of the metric fluctuations during a quasi-de Sitter expansion induced by the stress tensor fluctuations of the inflaton field. The results are in agreement with other approaches to the same problem , an approximate Harrison-Zel’dovich spectrum is predicted. We noticed that in our approach the quantum fields and the gravitational fields are treated separately, and this may have some advantages to go one step further and consider the quantum field fully, not just to linear order.

In is worth mentioning that other applications of the stochastic theory have been carried out and others are in progress. Thus, the fluctuations on the ground state of semiclassical gravity, which consists of the Minkowski metric and the quantum state in its vacuum state, have been considered . The computation of the two-point correlations of the linearized Einstein tensor indicate that a typical correlation length is present.

Other applications of stochastic semiclassical gravity to semiclassical cosmology have been performed , some including thermal fields . It has been shown that noise produced by a quantum field on the cosmological scale factor of an isotropic closed Friedmann-Robertson-Walker, in the presence of a cosmological constant, may take the scale factor from a region where it is nearly zero to a region where it describes a de Sitter inflationary era . Thus jumping over the barrier by activation, this is the semiclassical analogue of the tunneling from nothing in quantum cosmology  and gives yet another mechanism to produce inflation.

An important application of stochastic gravity which is now beginning is in the physics of black holes . In particular in black hole thermodynamics, the stress tensor fluctuations near the black hole horizon may induce fluctuations in the horizon area. The relevance of this back-reaction effect in Hawking radiation has not been yet explored, although preliminary investigations seem to indicate that Hawking result should not be substantially different . The contribution of the horizon fluctuations to the black hole entropy  is another tantalizing issue that may deserve some attention in the present context.

###### Acknowledgements.

I am very grateful to Professors Gr. Tsagas and D. Papadopoulos for giving me the oportunity to participate at the conference on Applied and differential geometry, Lie algebras and general relativity, and for their kind and genereous hospitality. I am also grateful to Rosario Martín and Albert Roura for their unvaluable contribution to the work described here. I also thank Albert Roura for a critical reading of the manuscript and many fruitful discussions. This work has been partially supported by the CICYT Research Project No. AEN98-0431.


format: text

text:
# Precise Measurement of Dimuon Production Cross-Sections in 𝜈_𝜇Fe and 𝜈̄_𝜇Fe Deep Inelastic Scattering at the Tevatron

## I Introduction

Oppositely-charged dimuon production in neutrino-nucleon deep inelastic scattering (DIS) provide an excellent source of information on the structure of the nucleon, the dynamics of heavy quark production, and the values of several fundamental parameters of the Standard Model (SM) of particle physics. These events are produced most commonly in charged-current (CC) neutrino DIS interactions when the incoming neutrino scatters off a strange or down quark to produce a charm quark in the final state, which subsequently fragments into a charmed hadron that decays semi-muonically. This distinct signature is easy to identify and measure in massive detectors, which allows for the collection of high statistics data samples. Consequently, dimuon events have played a significant role in the last 20 years in understanding charm production in DIS.

Since the contribution of the down quark to charm production is Cabibbo suppressed, scattering off a strange quark is responsible for a significant fraction of the total dimuon rate: $`90\%`$ in $`\overline{\nu }_\mu N`$ and $`50\%`$ in $`\nu _\mu N`$ scattering. Therefore the measurement of this process provides an excellent source for the determination of the strange quark parton distribution function (PDF) of the nucleon. Other methods involving measurements of the difference in parity violating structure functions from $`\nu `$ and $`\overline{\nu }`$ scattering off an isoscalar target and the ratio of parity conserving structure functions from charged lepton and neutrino DIS suffer from theoretical and experimental systematic uncertainties. Its prominent role in DIS charm production makes the $`s`$ quark PDF an important ingredient in tests of two-scale quantum chromodynamics (QCD), the two scales being the squared invariant momentum transfer $`Q^2`$ and the charm mass $`m_c`$. The $`s`$ quark PDF also enters into background calculations for new physics searches at the Tevatron and the LHC. For example, in some scenarios the super-symmetric top quark $`\stackrel{~}{t}`$ is best searched for via its loop decay to a charm quark and a neutralino, $`\stackrel{~}{t}c\stackrel{~}{\chi }^0`$. The dominant background to this mode is from the gluon-strange quark process $`gsW^{}c`$, with the $`W^{}`$ decaying leptonically. Accurate separate measurement of $`s`$ and $`\overline{s}`$ PDFs may also shed light on the interplay between perturbative and non-perturbative QCD effects in the nucleon.

In the following we present a new measurement of dimuon production from scattering of $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ beams from an iron target in the NuTeV experiment at the Fermilab Tevatron. This measurement exploits the high purity $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ beams produced by the Fermilab Sign-Selected Quadrupole Train (SSQT) to extend acceptance for muons down to $`5`$ GeV$`/c`$ momentum without any ambiguity in the determination of the muon from the primary interaction. Our analysis proceeds as follows: We first perform a fit to NuTeV $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ dimuon data using a leading order (LO) QCD cross-section model to obtain values for effective charm mass $`m_c`$, parameters describing the size and shape of the nucleon strange and anti-strange sea, and parameters describing the fragmentation of charm quarks into hadrons followed by their subsequent semi-muonic decay. The LO parameters allow us to make contact with previous measurements and provide an accurate description of the dimuon data. We then go on to use the LO model for acceptance and resolution corrections in extracting, for the first time, the cross-sections $`d\sigma \left(\nu _\mu /\overline{\nu }_\mu Fe\mu ^{}\mu ^\pm X\right)/dxdy`$ for forward secondary muons tabulated in bins of neutrino energy $`E_\nu `$, Bjorken scaling variable $`x`$, and inelasticity $`y`$. These cross-section tables provide the most model-independent convenient representation of $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ dimuon data. They may be used to test any cross-section calculation for dimuon production from iron, provided that the model is augmented with a fragmentation/decay package such as that provided by PYTHIA . We then repeat the cross-section extraction on an older CCFR data set from the same detector; and, after demonstrating consistency between the CCFR and NuTeV data, we combine the cross-section tables to obtain the most statistically precise high energy neutrino dimuon production collected to date. We finish by performing a LO QCD fit to the combined cross-section tables and extract precise determinations of $`m_c`$, the CKM parameter $`V_{cd}`$, and size and shape parameters for the strange sea.

## II Leading Order Charm Production

Dimuon production from charm depends on three different components: the charm production cross-section, the fragmentation of the charm quark to a charmed hadron, and the semileptonic decay of the charmed hadron. In LO QCD, charm production arises from scattering off a strange or down quark: $`\nu _\mu +\mathrm{N}\mu ^{}+c+\mathrm{X}`$, where the second muon is produced from the semileptonic decay of the charmed hadron (Fig. 1). In this case, the dimuon production cross-section factorizes into the form:

$$\frac{d^3\sigma _{\mu \mu }^{LO}(\nu _\mu N\mu ^{}\mu ^+X)}{dxdydz}=\frac{d^2\sigma _c^{LO}(\nu _\mu NcX)}{d\xi dy}\text{ }D(z)B_c.$$
(1)
Here, $`y`$ is the inelasticity, $`z`$ is the fraction of charm quark momentum carried by the charm quark hadron, and $`\xi =x\left(1+m_c^2/Q^2\right)\left(1x^2M^2/Q^2\right)`$ is the fraction of the nucleon’s momentum carried by the struck quark, with $`x`$ the Bjorken scaling variable and $`M`$ the nucleon mass. $`D(z)`$ is the fragmentation function for the charm quark; $`B_c`$, the semi-muonic branching fraction for charmed hadrons; both averaged over all charmed hadrons produced in the final state. The charm production cross-section $`d\sigma (\nu _\mu NcX)/d\xi dy`$ for an isoscalar target $`N`$ can be expressed as

$`{\displaystyle \frac{d^2\sigma _c^{LO}\left(\nu _\mu N\mu ^{}cX\right)}{d\xi dy}}`$ $`=`$ $`{\displaystyle \frac{2G_F^2ME/\pi }{\left(1+Q^2/M_W^2\right)^2}}\left(1{\displaystyle \frac{m_c^2}{2ME_\nu \xi }}\right)`$ (3)
$`\times \left[\left|V_{cs}\right|^2s(\xi ,Q^2)+\left|V_{cd}\right|^2\left(u(\xi ,Q^2)+d(\xi ,Q^2)\right)/2\right],`$
where $`G_F=1.17\times 10^5`$ GeV2, $`M_W=80.42`$ GeV$`/c^2`$, and $`V_{cs}`$/$`V_{cd}`$ are the $`cs`$/$`cd`$ CKM matrix elements. The corresponding $`\overline{\nu }_\mu `$ process has the quarks replaced by their anti-quark partners. One observes that the non-strange contributions to the cross-section are large in neutrino mode, where they dominate at high $`\xi `$; and small in anti-neutrino mode, where $`\overline{s}(\xi ,Q^2)`$ dominates at all $`\xi `$. In the case where $`s(\xi ,Q^2)=\overline{s}(\xi ,Q^2)`$, the neutrino $`\xi `$ distribution determines the relative size of the $`s(\xi ,Q^2)`$ and $`\left(u(\xi ,Q^2)+d(\xi ,Q^2)\right)/2`$ contributions, the anti-neutrino $`\xi `$ distributions determines the shape of $`s(\xi ,Q^2)`$, the energy dependence of the cross-section determines $`m_c`$, the ratio of the dimuon cross-section to the single muon inclusive cross-section sets $`B_c`$, and the energy distribution of the charm decay muon constrains $`D\left(z\right)`$.

The non-strange PDFs are determined from inclusive $`\nu N`$ and $`\overline{\nu }N`$ CC scattering, and from charged lepton scattering; they contribute negligible uncertainty to $`\sigma _c`$. The final state charmed hadron admixture is obtained from neutrino data from an emulsion experiment that has been corrected for improved knowledge of charmed hadron lifetimes. Extraction of CKM matrix elements from the data requires an independent determination of $`B_c`$; this analysis uses the Particle Data Group charm semi-leptonic branching fractions convolved with the species production cross-sections just mentioned. The result is $`B_c^{ext}=0.093\pm 0.009`$.

In next-to-leading-order (NLO) QCD, additional diagrams complicate the expression for $`\sigma _c`$ and also spoil the factorization of the various components of $`\sigma _{\mu \mu }`$. The charm production cross-section also becomes dependent on the QCD factorization and re-normalization schemes and their respective scales.

## III Experimental measurement and analysis technique

In an ideal situation one would like to present direct measurements of the differential charm production cross-sections $`d\sigma _c^{\nu \left(\overline{\nu }\right)N}/dxdy`$ at several different neutrino energies. Neither NuTeV nor its predecessor CCFR measures charm, but rather dimuons. The charm cross-section is thus related to the data by model-dependent corrections for charm fragmentation and decay and by experimental effects of resolution, acceptance, and neutrino flux. One way of handling these issues is to fit a parametric model directly to the data and extract parameters from the model. This approach was used in the past for LO QCD and NLO QCD in the variable flavor ACOT scheme.

The approach taken here begins with the same idea, a LO QCD parametric fit based on Eq. 1. Events passing selection criteria detailed in the next section are binned separately in $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ mode in the quantities

$`E_{VIS}`$ $`=`$ $`E_{\mu 1}+E_{\mu 2}+E_{HAD},`$ (4)
$`x_{VIS}`$ $`=`$ $`{\displaystyle \frac{4E_{\mu 1}\left(E_{\mu 1}+E_{\mu 2}+E_{HAD}\right)\mathrm{sin}^2\theta _{\mu 1}/2}{2M\left(E_{\mu 2}+E_{HAD}\right)}},\mathrm{and}`$ (5)
$`z_{VIS}`$ $`=`$ $`{\displaystyle \frac{E_{\mu 2}}{E_{\mu 2}+E_{HAD}}},`$ (6)
where $`E_{\mu 1}`$ is the energy of the primary muon with the same lepton number as the beam, $`E_{\mu 2}`$ is the energy of the other muon, $`E_{HAD}`$ is the observed hadronic energy in the calorimeter, and $`\theta _{\mu 1}`$ is the scattering angle of the primary muon. The “$`VIS`$” subscript indicates that these quantities differ from the true values of $`x`$, $`E`$, and $`z`$, due to the energy carried away by the neutrino from charm decay and due to detector smearing. Other quantities of interest for comparison purposes are

$`y_{VIS}`$ $`=`$ $`{\displaystyle \frac{E_{\mu 2}+E_{HAD}}{E_{\mu 1}+E_{\mu 2}+E_{HAD}}}\mathrm{and}`$ (7)
$`Q_{VIS}^2`$ $`=`$ $`2ME_{VIS}x_{VIS}y_{VIS}.`$ (8)
A binned likelihood fit is performed which compares the data to a model composed of a charm source described by Eqs. 1 and 3. The charm events are augmented with a contribution from dimuon production through $`\pi /K`$ decay in the charged-current neutrino interaction’s hadron shower and then processed through a detailed Monte Carlo (MC) simulation of the detector and the same event reconstruction software used for the data. The MC dimuon sample is normalized to the data through use of the inclusive single muon event rates in $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ mode. The fit varies a common charm mass $`m_c`$, branching fraction $`B_c`$, and fragmentation parameter $`ϵ`$ for both modes, and two parameters for each mode, $`(\kappa _\nu ,\alpha _\nu )`$ and $`(\kappa _{\overline{\nu }},\alpha _{\overline{\nu }})`$, that describe the magnitude and shape of the $`s`$ and $`\overline{s}`$ quark PDFs. The strange sea parameters are defined by

$`s(x,Q^2)`$ $`=`$ $`\kappa _\nu {\displaystyle \frac{\overline{u}(x,Q^2)+\overline{d}(x,Q^2)}{2}}\left(1x\right)^{\alpha _\nu }\mathrm{and}`$ (9)
$`\overline{s}(x,Q^2)`$ $`=`$ $`\kappa _{\overline{\nu }}{\displaystyle \frac{\overline{u}(x,Q^2)+\overline{d}(x,Q^2)}{2}}\left(1x\right)^{\alpha _{\overline{\nu }}}.`$ (10)
In these parameterizations***This parameterization differs slightly from previous LO analyses in the definition of $`\kappa _\nu `$ and $`\kappa _{\overline{\nu }}`$ in the general case of $`\alpha _\nu `$ and $`\alpha _{\overline{\nu }}0`$. The motivation for not using the older definitions (of the form $`s(x,Q^2)=\frac{\kappa _\nu }{2}\left(\overline{u}(x,Q^2)+\overline{d}(x,Q^2)\right)\left(1x\right)^{\alpha _\nu }`$ $`\times `$ $`\frac{_0^1𝑑x\left(\overline{u}(x,Q^2)+\overline{d}(x,Q^2)\right)}{_0^1𝑑x\left(\overline{u}(x,Q^2)+\overline{d}(x,Q^2)\right)\left(1x\right)^{\alpha _\nu }}`$) is to avoid a procedure that requires information about the PDF outside the experimentally accessible $`x`$ range of the experiment., values of $`\kappa _\nu =\kappa _{\overline{\nu }}=1`$ and $`\alpha _\nu =\alpha _{\overline{\nu }}=0`$ would imply an SU(3)-flavor symmetric sea; previous measurements have yielded $`\kappa `$ values around $`0.4`$, and $`\alpha `$ values consistent with zero (within large errors). The fragmentation process is described using the Collins-Spiller fragmentation function :

$$D(z,ϵ)=N\left[(1z)/z+ϵ_C(2z)/(1z)\right](1+z)^2[1(1/z)ϵ_C/(1z)]^2.$$
(11)
The Peterson function  was also tried but produced worse agreement between MC and data.

The analysis proceeds based on the observation that the dimuon data are well described by the LO fit. The information on charm production from the LO fit is, in fact, not of great importance; our goal is to use the fit parameters to help construct a cross-section. This task is performed by forming a new grid in $`x_{VIS}`$, $`y_{VIS}`$, and $`E_{VIS}`$ in the data and MC, and a corresponding grid of $`x,y,`$ and $`E`$ in the MC. The dimuon cross-section is computed at the weighted center of each $`(x,y,E)`$ bin $`i`$. The MC also predicts, with the result of the LO fit, the number of events in each $`(x_{VIS},y_{VIS},E_{VIS})`$ bin $`j`$. The MC can further be used to establish a correspondence between $`(x,y,E)`$ and $`(x_{VIS},y_{VIS},E_{VIS})`$ bins; this is accomplished by finding the bin $`j_i`$ in $`(x_{VIS},y_{VIS},E_{VIS})`$ space that receives the largest fraction of events produced in $`(x,y,E)`$ bin $`i`$. After this procedure the cross-section for $`(x,y,E)`$ bin $`i`$ is then determined by

$$\left\{\frac{d^2\sigma _{\mu \mu }(\nu _\mu N\mu ^{}\mu ^+X)}{dxdy}\right\}_{(x,y,E)_i}=\frac{D_i}{𝒩_{i,\text{LO fit}}^{}}_{E_{\mu 2}>E_{\mu 2\mathrm{min}}}𝑑z𝑑\mathrm{\Omega }\left\{\frac{d^3\sigma _{\mu \mu }(\nu _\mu N\mu ^{}\mu ^+X)}{dxdydz}\right\}_{(x,y,E)_i\text{LO fit}}.$$
(12)
In this expression the left hand side represents the measured cross-section for dimuon production as a function of $`x`$, $`y`$, and $`E`$ in bin $`i`$ with the requirement that the second muon in the event exceeds the threshold used in the experiment, $`E_{\mu 2}>E_{\mu 2\mathrm{min}}`$. On the right hand side, $`D_i`$ is the number of data corresponding to bin $`i`$, and $`𝒩_{i,\text{LO fit}}^{}`$ is the corresponding number of events predicted by the LO fit. In the integrand, $`d^3\sigma _{\mu \mu }(\nu _\mu N\mu ^{}\mu ^+X)/dxdydz`$ is taken from Eq. 1, and the integral over the fragmentation variable $`z`$ and charmed hadron decay variable $`\mathrm{\Omega }`$ maintains the condition $`E_{\mu 2}>E_{\mu 2\mathrm{min}}`$. The procedure for defining $`D_i`$ and $`𝒩_{i,\text{LO fit}}^{}`$ is detailed further in Appendix A.

The end result is two tables, one each for $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ mode, of the “forward” dimuon cross-section (which is closest to what the experiment actually measures). It will be shown later that these tables can be used to re-extract the LO fit parameters and can be combined with similar tables from the CCFR experiment. More details on the apparatus, event selection, and analysis procedure will be given first.

## IV NuTeV Experimental Details

### A Detector and beamline

The NuTeV (Fermilab-E815) neutrino experiment collected data during 1996-97 with the refurbished Lab E neutrino detector and a newly installed Sign-Selected Quadrupole Train(SSQT) neutrino beamline. Figure 2 illustrates the sign-selection optics employed by the SSQT to pick the charge of secondary pions and kaons which determine whether $`\nu _\mu `$ or $`\overline{\nu }_\mu `$ are predominantly produced. The SSQT produced beam impurities of $`\overline{\nu }_\mu `$ ($`\nu _\mu `$) events in $`\nu _\mu `$ mode ($`\overline{\nu }_\mu `$ mode) at the $`10^3`$ level. During NuTeV’s run the primary production target received $`1.13\times 10^{18}`$ and $`1.41\times 10^{18}`$ protons on target in neutrino and anti-neutrino modes, respectively, resulting in inclusive CC samples of $`1.3\times 10^6`$ events in neutrino mode and $`0.46\times 10^6`$ $`\overline{\nu }`$ inclusive CC events. The very low “wrong-flavor” backgrounds imply that only one muon charge measurement is needed to make the correct assignments for $`E_{\mu 1}`$ and $`E_{\mu 2}`$ described above.

The Lab E detector, described in detail elsewhere, consists of two major parts, a target calorimeter and an iron toroid spectrometer. The target calorimeter contains 690 tons of steel sampled at 10 cm(Fe) intervals by 84 3m $`\times `$ 3m scintillator counters and at 20 cm(Fe) intervals by 42 3m $`\times `$ 3m drift chambers. The toroid spectrometer consists of four stations of drift chambers separated by three iron toroid magnets that provide a $`p_T`$ kick of $`2.4`$ GeV$`/c`$. The toroid magnets were set to always focus the muon with same lepton number as the beam neutrino. Precision hadron and muon calibration beams monitored the calorimeter and spectrometer performance throughout the course of data taking. The calorimeter achieves a sampling-dominated resolution of $`\sigma _E/E=2.4\%87\%/\sqrt{E}`$ and an absolute scale uncertainty of $`\delta E/E=0.4\%`$. The spectrometer’s multiple-Coulomb-scattering-dominated muon momentum resolution is $`\sigma _p/p=`$ $`11\%`$, and the muon momentum scale is known to $`\delta p/p=1.0\%`$.

### B Data selection

A typical dimuon event has the characteristics shown in Fig. 3. In this figure, the toroid can be seen to focus the leading muon originating from the leptonic vertex and to de-focus the secondary muon, which originates most probably from charm decay. In the event shown both muons pass through the toroid, and both their signs are measured. In events where the sign of one muon is not measured, it is assumed to be opposite the one measured. Since the sign of the primary muon is known because of the sign selection of the SSQT, the measurement of the sign of only one muon is sufficient to identify the primary and secondary muon in the event. The rate of the same sign dimuon events with both muons toroid analyzed is very low.

Candidate opposite sign dimuon events, in both data and Monte Carlo, were selected using the following criteria:

* The event must occur in coincidence with the beam and fire the penetration trigger (charged-current interaction trigger).
* The incident neutrino energy, $`E_\nu `$, must be greater than $`20`$ GeV and the energy of the hadronic shower, $`E_{had}`$, greater than $`10`$ GeV.
* In order to ensure event containment, only events occurring within an active fiducial volume are accepted: the transverse vertex positions $`(V_x,V_y)`$ must be satisfy $`127`$ cm$`<V_{x,y}<127`$ cm and $`\sqrt{V_x^2+V_y^2}<152.4`$ cm, and the longitudinal vertex position must lie between counters 15 and 80, which corresponds to $`2.7`$ and $`13.3`$ hadronic interaction lengths from the upstream and downstream ends of the calorimeter, respectively.
* Two muons must be identified in the event and satisfy the following criteria:

+ The energy of each of the two muons, $`E_{\mu 1,\mu 2}`$, must be greater than 5 GeV.
    + The time obtained from fitting each track must be within 36 ns of the trigger time.
    + One of the muons must be toroid-analyzed and the energy of the toroid-analyzed muon at the entrance of the toroid must be greater than $`3`$ GeV.
    + At least one toroid-analyzed muon must pass through at least $`2/3`$ of the toroid.
    + The toroid-analyzed muons must hit the front face of the toroid inside a circle of radius $`R_{FF}<152`$ cm, and more than $`80\%`$ of the path length of the muon must be in the toroid steel.
* Finally, in order to remove mis-reconstructed events, a requirement is imposed on the reconstructed $`x_{VIS}`$ kinematic variable: $`0x_{VIS}1`$.

The final event sample contains $`5102`$ $`\nu _\mu `$-induced and $`1458`$ $`\overline{\nu }_\mu `$-induced events. Of these, $`2280/655`$ in $`\nu _\mu /\overline{\nu }_\mu `$ mode have both muons reconstructed in the toroid spectrometer. All other events have only one. The mean $`E_{VIS}`$ of the events is $`157.8`$ GeV, the mean $`Q_{VIS}^2=21.1`$ GeV2, and the mean $`x_{VIS}=0.14`$. The overall reconstruction efficiency, including the detector acceptance, is $`60\%`$ for events with $`E_{\mu 2}`$ 5 GeV, and $`80\%`$ when $`E_{\mu 2}`$ is above 30 GeV.

### C Detector Simulation

A hit-level Monte Carlo simulation of the detector based on the GEANT package  was used to model the detector response and provide an accurate representation of the detector geometry. The Monte Carlo events were analyzed using the same reconstruction software used in the data analysis. The detector response in the simulation was tuned to both hadron and muon test beam data at various energies. To ensure the accurate modeling of the muon reconstruction efficiency the drift chamber efficiencies implemented in the simulation were measured in the data as a function of time.

The primary neutrino interactions were generated using the LO QCD model and fragmentation function described in Sec. III. Electroweak radiative corrections based on the model by Bardin  were applied to this cross-section. The main background source from ordinary CC interactions in which a pion or kaon produced in the hadronic shower decays muonically was generated following a parameterization of hadron test beam muoproduction data for simulating secondary decays, and the LEPTO  package for the decays of primary hadrons . The total probability to produce such muons with momentum greater than 4 GeV/c is $`2\times 10^4`$ for events with $`E_{had}`$ 30 GeV, and $`10^3`$ for $`E_{had}`$ 100 GeV. The contribution from the primary hadron decays is roughly two times larger than that from the secondary decays.

### D Neutrino Flux and Normalization

The total flux, energy spectra, and composition for both neutrino and anti-neutrino beams are calculated using a Monte Carlo simulation of the beamline based on the DECAY TURTLE program  and production data from Atherton  as parameterized by Malensek  for thick targets. This flux is used to generate an inclusive charged-current interaction Monte Carlo sample using the GEANT based hit-level detector Monte Carlo described in section IV.C.

The predicted flux is then tuned so that the inclusive charged-current interaction spectra in the Monte Carlo match the data. Selection criteria for this sample of inclusive charged-current interactions are exactly the same as those used to select the dimuon sample with the requirement for two muons removed. Flux corrections of up to $`15\%`$ are applied in bins of neutrino energy and transverse vertex position to force the single muon data and Monte Carlo to agree. In addition, an overall factor is determined for each beam (neutrino or anti-neutrino) that absolutely normalizes the single muon Monte Carlo to the data. The dimuon Monte Carlo uses the flux determined with the above procedure; and it is absolutely normalized to the inclusive single muon charged-current data through the flux tuning procedure.

The procedure used to tune the flux to the observed single muon rate is iterative since the event rate observed in the detector depends on the convolution of the cross-section with neutrino flux, and the result of the charm measurement has a small effect on the total cross-section. Corrections found from the single muon Monte Carlo/data comparisons are used in the dimuon Monte Carlo that is used to determine the dimuon cross-section, and thus the charm production cross-section within our LO model. The charm cross-section results are then used in the single muon Monte Carlo and the procedure is repeated until the flux parameters do not change (in practice the convergence is very fast). Figure 4 shows a comparison between data and Monte Carlo for the single muon (flux) sample; the level of agreement is very good. Note that muon and hadron energies are not adjusted separately in the flux tuning procedure.

## V Results

### A NuTeV and CCFR Leading Order QCD fits

The LO QCD fits were performed using three different parton distribution function (PDF) sets with their corresponding QCD evolution kernels: GRV94LO  and CTEQ4LO , as implemented in the PDF compilation PDFLIB , and a Buras-Gaemers parameterization (BGPAR) that has been used extensively in this experiment and its CCFR predecessor. In the BGPAR case, an explicit Callan-Gross relation violation is implemented by replacing the term $`1\frac{m_c^2}{2ME_\nu \xi }`$ in Eq. 3 with $`(1+R_L)(1+(\frac{2M\xi }{Q})^2)^1(1y\frac{Mxy}{2E})+\frac{xy}{\xi }`$, where $`R_L`$, the ratio of longitudinal to transverse $`W^\pm N`$ cross-sections, is taken from a fit to electroproduction data. Table I lists fit results with the rightmost column showing the combined $`\chi ^2`$ for $`\nu `$ and $`\overline{\nu }`$ modes. All three models have the same good level of agreement with the dimuon data. Figure 6 illustrates the quality of the BGPAR fit by comparing it to the data for the kinematic variables used directly in the fit; figure 7 shows a comparison for variables not used directly in the fit. One observes that differences in the choice of PDF parameterization result in different charm production parameters, indicating significant model dependence at LO. Because of this dependence, the most relevant quantity to extract is the dimuon production cross-section.

### B CCFR Leading Order QCD fits

The CCFR dimuon data set used in the CCFR LO  and NLO  analysis is used to extract LO QCD parameters using the same procedure described in the previous section. Results, shown in Table II, are consistent with the NuTeV fits except for the fragmentation function shape parameter $`ϵ`$. This difference is caused by the fact that in NuTeV a much higher percentage of low energy muons is used, thus the $`z_{VIS}`$ distribution has significant shape differences to that of CCFR. In addition, the uncertainty in the $`\pi /K`$ background parameterization is larger for low muon energies; this effect is reflected in the size of the error on the determination of $`ϵ`$.

Table II shows results of a combined fit to the data from both NuTeV and CCFR performed using the same procedure and the BGPAR PDF set. Since the fragmentation shape parameter is different for the two experiments, the Monte Carlo sample for each is reweighted to the appropriate $`ϵ`$ and then kept fixed in the combined fit to simplify the fitting procedure.

### C Comparing CCFR with NuTeV

Different BGPAR PDF sets are used to analyze NuTeV and CCFR data. Thus, since we have shown that different PDF sets can produce different LO parameters, the numbers presented in Table II should not be compared directly with each other. Furthermore, the combined fit uses the NuTeV BGPAR PDF set, an approach that is not completely rigorous. Parameters in Table II should be compared to their corresponding values in Table III, which gives the result of the fit to the dimuon cross-section table discussed in the next section. In this case, the use of any PDF set or model is equally appropriate, since the purpose is to compare the results of the fit to the dimuon data to those of the fit to the extracted cross-section tables.

We caution readers against taking $`\kappa `$ and $`\alpha `$ parameters for strange seas extracted using the BGPAR sets and using them to construct $`s(x,Q^2)/\overline{s}(x,Q^2)`$ PDFs according to Eqs. 9-10 from a different PDF set. Since the NuTeV/CCFR PDFs are not publicly available, a safer course would be to take the strange sea parameters for CTEQ or GRV. Even then, it is important to match $`\kappa `$ and $`\alpha `$ from our fits with the appropriate PDF set.

### D Systematic uncertainties

Although the main thrust of the analysis is the extraction of the dimuon production cross-section, the various sources of systematic uncertainty are presented by listing their contributions to the LO fit parameters. This is done for reasons of clarity, since the individual systematic uncertainty contributions add too many entries in the cross-section tables. These methods are completely equivalent since the systematic uncertainty on the parameters of the LO fits propagates directly to the cross-section measurement.

The main sources of the systematic uncertainties arise from modeling uncertainties in the Monte Carlo. The most significant are the $`\pi /K`$ decay background simulation, the detector calibration from the analysis of test beam hadron and muon data as a functions of energy and position, and the overall normalization. In addition, in the case of the BGPAR fits, the uncertainty on the longitudinal structure function is important.

The systematic error sources are given in Table IV. For the combined NuTeV+CCFR fit, systematic errors begin to dominate the uncertainties for several fit parameters.

### E NuTeV Cross-Section Tables

As was shown in the previous section, charm cross-section parameters depend on the details of the charm production model used in the fit. In addition, a fragmentation and decay model must be used to extract the charm production cross-section from the observed dimuon rate, introducing further model dependence. These model dependences are exacerbated by the experimental smearing correction, due to the missing neutrino energy of the charmed hadron decay and the detector resolution, and the substantial acceptance correction for low energy muons. In order to minimize model dependencies, we choose to present our result in the form of a dimuon production cross-section.

We have shown in the previous section that we can obtain very good description of our dimuon data, independent of the charm production model assumptions in our Monte Carlo. We thus limit the use of the Monte Carlo to correct experimental effects to the measured dimuon rate, and we limit the measurement of this rate to regions of phase space where the acceptance of the experiment is high. In this case, the only model dependence comes from potential smearing effects close to our acceptance cuts, i.e., the uncertainty on production of events produced with kinematic variables outside our cuts that smeared to reconstructed values within the cuts. The prediction for this kind of smearing depends on the underlying physics model, which is not well constrained by our data, since it involves phase space not accessed by our data. However, this model dependence is a second order effect.

The dimuon cross-section extraction procedure depends on the ability of our Monte Carlo to describe the data. The “true” three-dimensional phase space $`(x,y,E)`$ is divided into a number of bins. The grid is set up so that there is the same number of Monte Carlo events in each bin for any projection into one dimension. The cross-section for dimuon production with $`E_{\mu 2}5`$ GeV is calculated according to Eq. 12 at the center of gravity of each bin. The size of the bins in “visible” phase space $`(x_{VIS},y_{VIS},E_{VIS})`$ is defined in such a way so that there is a correspondence between a visible bin $`j_i`$ and a generated $`i`$ established by the condition that bin $`j_i`$ contains at least $`60\%`$ of the events from bin $`i`$. The relative error on the cross-section is taken to be equal to $`\delta _{\sigma _2}/\sigma _{2\mu }=\sigma _{2\mu }/\sqrt{D\left(j_i\right)}`$. Cross-section tables were produced using all three BGPAR, GRV, and CTEQ models; figure 8-11 demonstrate the insensitivity of the cross-section to PDF choice. Further details of the binning procedure are defined in Appendix A.

### F Charm Production Fits Using the Dimuon Cross-Section Tables

Using the cross-section tables involves similar steps as used in the direct data analysis, but with all detector and flux dependent effects removed. One must provide a model for charm production, fragmentation, and decay; construct the dimuon cross-section number for each entry in the table; and perform a $`\chi ^2`$ fit. The $`\chi ^2`$ function should employ the statistical and systematic errors in each bin added in quadrature.

One must also account for correlations between the various table entries. These correlations derive from our use of the LO fit to parameterize the data and from our method of binning; they are an inherent consequence of the incomplete kinematic reconstruction of the dimuon final state. We have adopted a pragmatic approach towards handling this issue. Rather than compute a large correlation matrix, we inflate the cross-section errors in each bin by a factor that is typically $`1.4.`$ This factor is chosen so that an uncorrelated $`\chi ^2`$ fit to the tables returns the same parameter errors as a direct fit to the data. A consequence of this error inflation is an apparent best fit $`\chi ^2`$ per bin $`\left(\chi ^2/\text{bin}\right)`$ that is approximately $`1/4`$ rather than $`1`$. This low value results from overcounting the number of degrees of freedom (DOF) present in the table and is compensated for by using an effective DOF per bin determined via MC calculation and given along with cross section data in the tables. Further details of the determination of the effective DOF are presented in Appendix A.

We tested this fitting procedure on the tables using the same BGPAR, fragmentation, and decay models used to obtain the table; Table III summarizes this study. Both the parameter values and their uncertainties obtained from fitting to the cross-section table agree with the corresponding values obtained by fitting directly to the data. It has also been verified that GRV94 parameters, for example, can be extracted from a cross-section table constructed with either the CTEQ or BGPAR model so as to agree with parameters obtained by fitting directly to the data.

While our cross-checks in fitting the cross-section tables entail using the same physics model used in generating the tables, we emphasize that the table presents a set of physical observables which may be used to test any dimuon production model. For the most interesting case of dimuon production through charm, a typical model test would consist of the following stepsA simple PYTHIA implementation of the first three steps may be obtained at www-e815.fnal.gov, or by contacting the authors.:

* For each bin, generate a (large) ensemble of $`N_{GEN}`$ events with $`x`$, $`y`$, and $`\stackrel{}{p}_c`$, where $`\stackrel{}{p}_c`$ denotes the lab momentum of the produced charm quark, according to the model charm production differential cross-section $`d\sigma _{cmod}/dxdyd\stackrel{}{p}_c`$
* Fragment the charm quarks into hadrons and decay the charmed hadrons (using, for example, PYTHIA), and determine $`N_{PASS}`$ the number of events which have a charm decay muon with $`E_{\mu 2}5`$ GeV.
* The cross-section to compare to the table value is then $`d\sigma _{\mu \mu mod}/dxdy=\left(N_{PASS}/N_{GEN}\right)\times 𝑑\stackrel{}{p}_c𝑑\sigma _{cmod}/𝑑x𝑑y𝑑\stackrel{}{p}_c`$.
* A fit should then be performed to minimize

$$\chi ^2=\underset{tablebins}{}\frac{\left(d\sigma _{\mu \mu mod}/dxdyd\sigma _{\mu \mu data}^+/dxdy\right)^2}{\sigma _{stat}^2+\sigma _{syst}^2}$$
(13)
in each beam mode with respect to the desired parameters in $`d\sigma _{cmod}/dxdyd\stackrel{}{p}_c.`$
* The confidence level for the fit may be obtained by comparing the $`\chi ^2`$ to the sum of effective DOF for table bins used in the fit.

To provide further experimental information for model testing, the following kinematic quantities are given along with the cross-section in each bin: $`E_{HAD}`$, the mean visible hadronic energy; $`E_{\mu 2}`$, the mean energy of the secondary muon; $`p_{T2in}^2`$, the mean square of the secondary muon’s transverse momentum in the event scattering plane; and$`p_{T2out}^2`$, the mean square of the secondary muons transverse momentum perpendicular to the event scattering plane. These quantities are computed from the dimuon data, with the LO fit used only for acceptance and smearing corrections. Tables V-XVI contain the measurementsThese data may also be obtained electronically at www-e815.fnal.gov.. Appendix B contains a supplementary discussion of the cross sections at high $`x`$.

## VI Summary

We present a measurement of the dimuon production cross-section from an analysis of the data of the NuTeV neutrino DIS experiment at the Tevatron. NuTeV data are combined with an earlier measurement from the CCFR experiment that used the same detector but a different beamline. A leading order QCD analysis of charm production performed on the combined data yields the smallest errors to date on model parameters describing the charm mass, the size and shape of the strange sea, and the mean semi-muonic branching fraction of charm. The leading order QCD model describes NuTeV and CCFR data very well, but cross-section model parameters extracted vary depending on the particular choice of model. The extracted dimuon production cross-section, by contrast, is insensitive to the choice of leading order QCD model and should be of most use to the community of phenomenologists.

## A Cross-Section Table Binning Procedure

In this analysis we report the differential cross-sections $`d\sigma \left(\nu _\mu /\overline{\nu }_\mu Fe\mu ^{}\mu ^\pm X\right)/dxdy`$ for forward secondary muons tabulated in bins of neutrino energy $`E`$, Bjorken scaling variable $`x`$, and inelasticity $`y`$. The measurement is obtained by using a LO Monte Carlo fit to the data to find the correspondence (mapping) between the “true” (unsmeared) and reconstructed (smeared) phase space, as discussed in Sec. III. This procedure maps the statistical fluctuations of the observed event rate to the “true” phase space bins where the cross-section is reported. The consistency of the procedure is checked by comparing fits of various models to the extracted cross-section tables to fits of the same models performed directly to the data. The criteria for the check is that the obtained central values and the errors on the model parameters are the same in both cases. In order to meet these criteria the binning of the smeared and the unsmeared phase space has to be appropriately selected.

Usually in cross-section measurements, the procedure followed is to bin both the “true” and the reconstructed variables using the same grid, and select the bin size empirically so that for each bin the purity is maximized and the smearing contribution from other bins is minimized§§§Purity is defined here as the fraction of events which have both their unsmeared and smeared variables belonging to the same bin.. In such a method the correspondence between smeared and unsmeared bins is trivial since the same grid of bins is used. In the case of the dimuon cross-section measurement, a complication arises from the large smearing effects due to the missing neutrino energy in the reconstructed dimuon final state. Unlike detector resolution effects, this smearing is not a symmetric function of the “true,” unsmeared variables, but rather an asymmetric mapping similar to electroweak radiative corrections. Because of this effect, we have followed a more elaborate procedure to map the visible phase space bins to those of the “true” phase space. This procedure allows us to obtain the desired high purity for each bin and also achieve stability of the result independently of the binning choice. In addition to mapping, our procedure allows us to take into account the significant bin-to-bin correlations which arise from the large smearing corrections without having to construct a full error matrix. This correlation matrix can be calculated, but it is too unwieldy to be useful; and it is difficult to incorporate effects of correlated systematic errors in a meaningful way. In our treatment, we estimate an effective number of degrees of freedom which allows us to obtain from an uncorrelated fit to the cross-section tables the same fit parameter errors as the ones obtained by directly fitting to the data, using the usual $`\mathrm{\Delta }\chi ^2=1`$ definition.

We begin the description of the technique by defining more precisely the factors $`D_i`$ and $`𝒩_{i,\text{LO fit}}^{}`$ from Eq. 12. In this equation, the expression on the left hand side represents the measured cross-section for dimuon production as a function of $`x`$, $`y`$, and $`E`$ in bin $`i`$ of the true phase space (with the requirement $`E_{\mu 2}>E_{\mu 2\mathrm{min}}`$). On the right hand side, $`D_i`$ is the number of data events in bin $`i`$ (which has to be determined using the Monte Carlo mapping procedure), and $`𝒩_{i,\text{LO fit}}^{}`$ is the number of events predicted in bin $`i`$ by the LO fit. To estimate the number of data events associated with bin $`i`$, we start by selecting bins in generated phase space $`(x,y,E)_i`$, requiring equal number of Monte Carlo events in each projection of the true (generated) phase space, so the number of events in each $`(x,y,E)_i`$ bin is approximately the same. The visible phase space is divided using the same algorithm but with a much finer grid than the generated space. The mapping matrix $`_{ij}`$ makes the correspondence between visible and generated phase spaces:

$`_{ij}=𝒩_{ij}/𝒩_j,`$

where $`𝒩_{ij}`$ is the number of Monte Carlo events generated in bin $`i`$ which end up in visible bin $`j`$ and $`𝒩_j`$ is the total number of Monte Carlo events in visible bin j. The coverage fraction $`𝒞`$ in visible space is defined as

$$𝒞=\underset{j}{\overset{N(𝒞)}{}}𝒩_{ij}/𝒩_i,$$
(A1)
where $`𝒩_i`$ is the number of Monte Carlo events in generated bin $`i`$. In the case of $`𝒞=1`$ the sum is performed over all visible bins; otherwise, the summation goes over the bins with the highest $`𝒩_{ij}/𝒩_i`$ ratios until the desired fractional coverage is obtained. Using the above definitions, and for a given coverage $`𝒞`$, we can define the number of data events which belong to a given true phase space bin $`i`$ (where the cross-section is reported) as

$$D_i=\underset{j𝒞}{}_{ij}D_j.$$
(A2)
The number of Monte Carlo events in this generated bin, for a given coverage $`𝒞`$, is redefined as

$$𝒩_i^{}=\underset{j𝒞}{}_{ij}𝒩_j.$$
(A3)
The cross-section error for each bin $`i`$ should be proportional to the visible events “mapped” in that bin, so the following expression is used to assign it:

$$\delta _i=\frac{d\sigma _{\mu \mu }(\nu _\mu N\mu ^{}\mu ^+X)_i}{dxdy}/\sqrt{𝒩_i}.$$
(A4)
Note that the Monte Carlo statistics contribution comes from the total number of Monte Carlo events generated in bin $`i`$. The multiplicative factor $`D_i/𝒩_i^{}`$ in Eq. 12 cancels out to first order the model dependence of the extracted cross-section and approximately transfers the statistical fluctuation in the visible phase space to the true phase space.

The procedure as described above is incomplete, since the true bins $`i`$ are not statistically independent. As we stated in the introduction to this section, we do not calculate the full error matrix, but rather estimate an effective (independent) number of degrees of freedom. It is possible to estimate the number of independent degrees of freedom by calculating the contribution to the total number of degrees of freedom from each bin

$$DOF_i=\frac{_{i𝒞}_{ij}𝒩_i}{_{j𝒞}𝒩_j}.$$
(A5)
It is obvious that the effective number of degrees of freedom depends on the selected coverage fraction $`𝒞`$, and should decrease as $`𝒞`$ increases. Figure 12 shows the number of effective degrees of freedom (DOF) as a function of the coverage area $`𝒞`$ (solid curve). The dotted curve in Fig. 12 shows the $`\chi ^2`$ obtained as a result of the fit to the table. We conclude that our method produces the correct number for $`\chi ^2/DOF`$, if the effective number of degrees of freedom is used, for a coverage fraction between $`55\%`$ and $`90\%`$.

The coverage area percentage used for our reported cross-section result is based on a Monte Carlo study. In this study a Monte Carlo sample is used to produce “cross-section” tables and then fits are performed to the tables and directly to the sample. Using as guidelines the criteria that there should be no pull on fit parameters in a fit to the cross-section table versus a direct fit, and that an uncorrelated fit to the tables should yield the same fit parameter errors as the ones obtained by a direct fit, we selected a coverage area $`𝒞=60\%`$. The effective number of degrees of freedom which corresponds to this value should be used with all fits performed on the cross-section tables we present in this article.

## B The high $`x`$ region

This section presents a supplementary investigation of the high $`x`$ region ($`x>0.5`$). The objective of this study is to ascertain whether there is any indication of an enhancement in the cross-section that we may be missing due to our use of wide $`x`$ bins (dictated by the low observed event rate at high $`x`$). Such an enhancement could be caused by an unusually large strange sea, particularly in neutrino mode, which has been advocated to resolve certain discrepancies between inclusive charged lepton and neutrino scattering . Previous dimuon analyses  may have missed this effect due to dependence on the particular model used to parameterize the strange sea distribution.

In order to minimize model-dependent corrections, we report our high x cross-section measurements as fractions of the total dimuon cross-section. Similarly, to quote a limit for the $`x>0.5`$ cross-section we use the observed data rate for $`x_{VIS}>0.5`$. This is a conservative way to set a limit, since by the kinematic effect of the missing decay neutrino energy, the contribution to a given $`x_{VIS}`$ bin always comes from $`x<x_{VIS}`$.

The cross-section ratio of the dimuon cross-section for $`x>0.5`$ to the total dimuon cross-section in a given energy bin can be expressed as

$$\frac{\sigma _{x>0.5}^{2\mu }}{\sigma ^{2\mu }}=\frac{N_{x>0.5}}{N_{tot}}\frac{M}{M_{E>5}},$$
(B1)
where $`N_{x>0.5}`$ is the number of observed events for $`x_{VIS}>0.5`$, $`N_{tot}`$ is the total number of observed dimuon events, $`M`$ is the Monte Carlo prediction with all experimental cuts applied, and $`M_{E>5}`$ is the Monte Carlo prediction without the $`\pi /K`$ decay contribution and with only the $`E_{\mu 2}>5`$ GeV cut applied. For this study, we use the same data selection criteria described in section IV B, exept for the $`x_{VIS}`$ selection, which is changed to: $`0x_{VIS}2`$.

In the anti-neutrino dimuon data sample we define three energy bins and we record the number of the observed dimuon events with $`x_{VIS}>0.5`$ in the data, together with the Monte Carlo prediction and all the relevant information for Eq. B1. Here the Monte Carlo prediction is very well constrained by our dimuon data in the full $`x`$ range, since the observed rate is mostly due to scattering on $`\overline{s}`$ quarks. The results are presented in Table XVII. The systematic and statistical error on the Monte Carlo prediction is negligible for this discussion. Treating the Monte Carlo prediction as a “background,” and using Eq. B1, we set cross-section ratio upper limits at 90% CL, for any additional source of $`x>0.5`$ dimuons, of 0.0012, 0.007, and 0.009, respectively, in each one of the energy bins defined in Table XVII (counting from lower to higher energy bin).

We follow the same procedure in the neutrino data sample. Here, there is an additional complication in the interpretation of the result since a large contribution from valence quark events is expected. The Monte Carlo prediction for the rates of the non-strange sea component is not directly constrained by our dimuon data, but rather by inclusive structure function measurements. The results are presented in Table XVIII. It is worth noticing that within the model and PDF sets used in this analysis (BGPAR) the Monte Carlo prediction for the contribution of the strange sea is only on the order of 3.7% of the total rate for $`x_{VIS}>0.5`$; most of the contribution (69%) in this model comes from the valence quarks.

Using Eq. B1 and treating the Monte Carlo prediction as a “background” source, we set 90% CL limits for an additional cross-section source at $`x>0.5`$. We find that for the first and last energy bins in table XVIII this additional source cannot be larger than 0.006 and 0.013 of the total dimuon cross-section, while for the 153.9-214.1 bin there is less than 5% probability that there is an additional source consistent with our data (note that we have a 1.75$`\sigma `$ negative yield compared to our background prediction).