text:
# KEK-TH-742 On the Ambiguity of the Solution of the Muskhelishvili-Omnes Integral Equation

(January 2001)

The solution of the Muskhelishvili-Omnes Integral Equation is ambiguous by a real polynomial. The coefficients of this polynomial can be fixed either by the knowledge of the low energy parameters or by the asymptotic behavior of the form factor. The role of the contact terms of the low energy effective Lagrangian is explicitly analysed.

There have been discussions of contact terms in the phenomenological models describing low energy pion phenomena . They seem to play an important role in analysing the pion form factor and the process $`\gamma \pi ^0\pi ^+\pi ^{}`$. The main difference between this model and the usual vector meson model is the existence of a contact term where the photon interacts directly with the charged pions and not by the intermediate of the vector $`\rho `$ meson. Phenomenological analysis of this contact term in the form factor has recently been given .

The purpose of this note is to analyze the existence of the contact term within the spirit of the more general approach to the form factor problem using the Muskhelishvili-Omnes (MO) integral equation . It is based on the assumption of the elastic unitarity and dispersion relation for the pion form factor. The main conclusion of our analysis is that the existence of the contact term is in conflict with unitarity in the form of the phase theorem. Its existence can, however, be simulated by a polynomial ambiguity contribution to the solution of the MO equation. This polynomial is interpreted physically as either the correction due to the inelastic contribution to the unitarity relation or the modification of the asymptotic behavior of the form factor.

Let us first discuss the solution of this problem from a more general viewpoint. The form factor $`V(s)`$ is an analytic function with a cut from $`4m_\pi ^2`$ to $`\mathrm{}`$.

Assuming that $`V(s)`$ is polynomially bounded and that $`V(s)s^{(n+1)}0`$ as $`s\mathrm{}`$, $`n>0`$, we can write the following subtracted dispersion relation:

$$V(s)=1+a_1s+\mathrm{}a_ns^n+\frac{s^{n+1}}{\pi }_{4m_\pi ^2}^{\mathrm{}}\frac{ImV(z)dz}{(z)^{n+1}(zsiϵ)}$$
(1)
where we have used the Ward identity $`V(0)=1`$. Around $`s=0`$, the dispersion integral is of the order $`(s)^{n+1}`$ and can be neglected, hence the low energy theorem is recovered. Needless to say, $`a_n`$ are, apart from a factorial factor $`n!`$, the derivatives of $`V(s)`$ evaluated at $`s=0`$. These could be taken from experiments or could be determined by sum rules using some assumptions on the asymptotic behavior of the form factor.

For example, if we assume that as $`s\mathrm{}`$, $`V(s)/s0`$, a once subtraction dispersion relation for $`V(s)`$ can be written; the first and higher derivatives of $`V(s)`$ can be written as:

$$V^{(n)}(0)=\frac{n!}{\pi }_{4m_\pi ^2}^{\mathrm{}}\frac{ImV(z)dz}{z^{n+1}}$$
(2)
The mathematical problem is now clear: Find the solution of the integral equation of the MO type for $`V(s)`$ with its imaginary part given by the elastic unitarity:

$$ImV(s)=V(s)e^{i\delta (s)}\mathrm{sin}\delta (s)$$
(3)
with the boundary conditions around $`s=0`$ given by Eq. (1) and where $`\delta `$ is the the strong P-wave $`I=1`$ $`\pi \pi `$ phase shifts.

To solve this integral equation Eq. (1) which is of MO type , let us define the function $`\mathrm{\Omega }(s)`$ normalized to unity at $`s=0`$:

$$\mathrm{\Omega }(s)=\mathrm{exp}(\frac{s}{\pi }_{4m_\pi ^2}^{\mathrm{}}\frac{\delta (z)dz}{z(zsiϵ)})$$
(4)
The solution for our integral equation is:

$$V(s)=P_n(s)\mathrm{\Omega }(s)$$
(5)
where $`P_n(s)`$ is a polynomial in $`s`$ of order n with real coefficients, $`P_n(0)=1`$. From eq. (5), the phase of the form factor $`V(s)`$ is $`\delta (s)`$ which is also a direct consequence of the unitarity relation . It will be called in the following as the phase theorem.

The ambiguity of the solution of the MO equation due to the polyomial is not surprising because its solution depends on the knowledge of the asymptotic behavior of the form factor. The polynomial could be interpreted as extra parameters which are introduced in the problem without violating the phase theorem; in many cases they can be determined using experimental data at low energy.

This is indeed true by expanding the function $`\mathrm{\Omega }(s)`$ in a power series in $`s`$, and compare it with the derivatives given in Eq. (5) and Eq. (1). The expansion of $`\mathrm{\Omega }(s)`$ in the Taylor’s series around $`s=0`$ is possible because $`\mathrm{\Omega }(s)`$ is an analytic function with a cut from $`4m_\pi ^2`$ to $`\mathrm{}`$.

In the special case where only two terms in the series are known, i.e. the pion charge by Ward’s indentity and its first derivative $`a_1`$ which is proportional to the r.m.s. pion radius by experiment, the solution of our integral equation is given by:

$$V(s)=\{1+s(a_1\mathrm{\Omega }^{^{}}(0))\}\mathrm{\Omega }(s)$$
(6)
where $`\mathrm{\Omega }^{^{}}(0)`$ denotes the first derivative of $`\mathrm{\Omega }(s)`$ evaluated at $`s=0`$ and is given by:

$$\mathrm{\Omega }{}_{}{}^{}(0)=\frac{1}{\pi }_{4m_\pi ^2}^{\mathrm{}}\frac{\delta (z)dz}{z^2}$$
(7)
The presence of the term $`\mathrm{\Omega }^{^{}}(0)`$ is to ensure the boundary condition for $`V(s)`$ is satisfied. If $`a_1\mathrm{\Omega }^{^{}}(0)=0`$ there is no need to introduce a momnomial term in the form factor and hence the asymptotic behavior of the form factor is solely determined by $`\mathrm{\Omega }(s)`$. (Using partial fraction for the dispersion integral together with the sum rules Eq. (2) the asymptotic behavior of $`V(s)`$ can be determined).

It is straightforward to generalise the solution of Eq. (1) for other values of $`n`$. For example when $`n=2`$, the solution for the integral equation Eq.( 1) is obtained by adding to the curly bracket on the righthand side of Eq. (6) a term:

$$s^2(a_2a_1\mathrm{\Omega }^{^{}}(0)a_0\frac{\mathrm{\Omega }^{^{\prime \prime }}(0)}{2}+\mathrm{\Omega }^{}_{}{}^{}2(0))$$
(8)
Let us now look at the experimental data. For a given set of P-wave $`\pi \pi `$ phase shifts $`\delta `$, one can, in principle, construct the function $`\mathrm{\Omega }(s)`$ as given by Eq. (4) by the numerical method. We prefer instead to construct this function using the inverse amplitude method or the  Padé approximant method for the one loop Chiral Perturbation Theory for the pion form factor and verify that the resulting phases are in agreement with experiments . The function $`\mathrm{\Omega }(s)`$ can be parametrised as follows:

$$\mathrm{\Omega }(s)=\frac{1}{1s/s_R\frac{1}{96\pi ^2f_\pi ^2}\{(s4m_\pi ^2)H_{\pi \pi }(s)+2s/3\}}$$
(9)
where $`f_\pi =0.093GeV`$, and $`s_R`$ is related to the $`\rho `$ mass squared $`m_\rho ^2`$ by requiring that the real part of the denominator of Eq. (9) vanishes at the $`\rho `$ mass; $`H_{\pi \pi }(s)`$ is a well-known integral over the phase space factor:

$$H_{\pi \pi }(s)=(22\sqrt{14m_\pi ^2/s})\mathrm{log}\frac{\sqrt{s}+\sqrt{s4m_\pi ^2}}{2m_\pi }+i\pi \sqrt{14m_\pi ^2/s}$$
(10)
for $`s>4m_\pi ^2`$; for other values of s, $`H_{\pi \pi }(s)`$ can be obtained by analytic continuation. $`\mathrm{\Omega }(s)`$ as given by Eq. (9) has a ghost at a very large negative value of $`s`$ with a small residue and is irrelevant for our low energy calculation.

Although we do not try the best fit to the experimental data on the P-wave $`\pi \pi `$ phase shifts, an adequate agreement is obtained between theory and experiment as shown in Fig. (1) where $`m_\rho =0.775GeV`$.

Eq. (9) can also be obtained by the vector meson dominace model  with the strong $`\rho \pi \pi `$ coupling constant given by the KSRF relation .

If it was assumed that there is no polynomial ambiguity $`P_n(s)=1`$, we would have:

$$V(s)=\mathrm{\Omega }(s)$$
(11)
Fitting this expression with the physical value of the $`\rho `$ mass, the pion rms radius calculated by this expression is too low by more than $`10\%`$ and the maximum value of the form factor squared at the $`\rho `$ peak (the $`\rho `$ leptonic width) is too low by $`30\%`$, Fig.(2). This implies that the assumed elastic unitarity relation is not good enough. One can either rewrite the integral equation taking into account of the inelastic effect which was previously studied  or to make more subtractions in the integral equation and keeping the elastic unitarity relation, in order to minimize the high energy contribution. We choose first the latter method. (The two-loop inverse amplitude or the  Padé approximant method are given recently and are found to give excellent agreement with the pion form factor data ).

Let us write a twice subtracted dispersion relation, including the inputs as the pion charge and the experimental value of the pion radius. The integral equation now reads:

$$V(s)=1+a_1s+\frac{s^2}{\pi }_{4m_\pi ^2}^{\mathrm{}}\frac{V(z)e^{i\delta (z)}\mathrm{sin}\delta (z)dz}{z^2(zsiϵ)}$$
(12)
where $`a_1=1/6(<r^2>)`$ with $`<r^2>`$ is the rms squared radius of the pion. The factor $`z^2`$ in the denominator of the dispersion integral suppresses the high energy contribution.

Using the experimental value $`<r^2>=0.439\pm .008fm^2`$  , the solution of the integral solution is given by Eq. (6) and is numerically equal to

$$V(s)=(1+0.15s/m_\rho ^2)\mathrm{\Omega }(s)$$
(13)
The $`\rho `$ leptonic width and the square of the modulus of the pion form factor are now in agreement with the experimental data, Fig.(2).

The other possibility to get agreement with experimental data is to introduce the inelastic effect which could occur at higher energy. The integral equation now reads:

$$V(s)=1+\frac{s}{\pi }_{4m_\pi ^2}^{\mathrm{}}𝑑z\frac{f^{}(z)V(z)+\sigma (z)}{z(zsiϵ)}$$
(14)
where $`f(s)=(2i)^1(\eta e^{2i\delta (s)}1)`$, $`\eta `$ is the inelastic factor ($`\eta =1`$ for elastic $`\pi \pi `$ scattering) and $`\sigma `$ is the inelastic spectral function due to the contribution of higher intermediate states in the unitarity relation. The solution of this integral equation is known:

$$V(s)=\mathrm{\Omega }(s)\{1+\frac{s}{\pi }_{s_i}^{\mathrm{}}𝑑z\frac{2Re(\sigma (z)e^{i\delta (z)})}{(1+\eta (z))e^{i\delta (z)}\mathrm{\Omega }(z)}\frac{1}{z(zsiϵ)}\}$$
(15)
$`s_i`$ is the inelastic threshold and for the pion form factor problem one can take it to be around $`11.3GeV^2`$. In reality, it could be as high as 2 $`GeV^2`$ in the $`\rho ^{}`$ resonance where the inelastic effect would become important.

The dispersion integral in Eq. (15) is seen to be an analytic function with a cut starting at $`s_i`$ to $`\mathrm{}`$ and hence it can be expressed as a Taylor’s series around $`s=0`$. The radius of convergence of this series is $`s_i`$. Hence the RHS of Eq. (15) is simply the solution of the MO equation as given by Eq. (5) and the fit to the experimental data is given by Eq. (13).

For a feel for the inelastic effect, let us see how does the $`\rho {}_{}{}^{}(1500)`$ , which dominantly decays into $`4\pi `$, influence the pion form factor. The inelastic spectral function can then be approximated by a delta function in the dispersion integral which yields a term proportional to $`\mathrm{\Omega }(s)`$ multiplying with the $`\rho ^{}`$ propagator. This assures the validity of the phase theorem below the inelastic threshold and that neither the $`\rho ^{}`$ propagator alone nor a pure contact term can be written. For $`s`$ in the $`\rho `$ region, one can expand the $`\rho ^{^{}}`$ propagator in a series in $`s`$ and fit to the experimental data to get the desired result Eq. (13).

Let us now compare our result with models with a contact term  and, in particular, the hidden local symmetry model . At the tree level, the hidden local symmetry model gives:

$$V^{HLS}(s)=\frac{a}{2}+1+\frac{a}{2}\frac{m_\rho ^2}{m_\rho ^2s}$$
(16)
which can be rewritten as:

$$V^{HLS}(s)=\frac{m_\rho ^2+s(a/21)}{m_\rho ^2s}$$
(17)
In order to make the $`\rho `$ unstable we could try to replace the factor $`m_\rho ^2/(m_\rho ^2s)`$ in Eq. (16) by $`\mathrm{\Omega }(s)`$ leaving the contact term unchanged. This would violate the phase theorem unless the constant term cancels each other to give a net result which is equal to $`\mathrm{\Omega }(s)`$; this solution does not fit with the experimental data as explained above.

We should use instead Eq. (17) and interprete it as the limit large $`s`$, $`V^{HLS}(1a/2)`$ (in the limit of zero width). Then the hidden symmetry solution for the vector form factor must contain a polynomial $`P_1(s)`$ or:

$$V^{HLS}=(1+\frac{a/21}{m_\rho ^2}s)\mathrm{\Omega }(s)$$
(18)
Compare this equation with Eq. (13) we have $`a=2.30`$. Within the validity of the above interpretation of the contact term, this result is in agreement with that obtained by the authors of the reference  who begin with Eq. (16), then make the $`\rho `$ unstable. The unitarisation scheme is done with the $`N/D`$ method for the P-wave $`\pi \pi `$ scattering.

The limit of large value of $`s`$ of $`V^{HLS}`$ should be taken with some reservation, because the inelastic effect interpretation of the polynomial ambiguity is only valid for $`s<s_i`$, Eq. (15) and $`s_i`$ is experimentally not very large.

Another example of the contact term  is discussed in the $`\gamma \pi ^0\pi ^+\pi ^{}`$ process, in particular, the hidden symmetry model . In the VMD model, the $`\gamma \pi ^0\pi ^+\pi ^{}`$ is written as:

$$A^{3\pi }(s,t,u)=\frac{\lambda }{3}(\frac{m_\rho ^2}{m_\rho ^2s}+\frac{m_\rho ^2}{m_\rho ^2t}+\frac{m_\rho ^2}{m_\rho ^2u})$$
(19)
where $`s,t,u`$ are the standard invariant kinematics. With a contact term it can be written as:

$$A^{3\pi }(s,t,u)=\frac{\lambda }{3c}\{\frac{m_\rho ^2}{m_\rho ^2s}+\frac{m_\rho ^2}{m_\rho ^2t}+\frac{m_\rho ^2}{m_\rho ^2u}c\}$$
(20)
where $`c`$ is proportional to the strength of the contact term and $`\lambda `$ is given by the anomaly:

$$\lambda =\frac{e}{4\pi ^2f_\pi ^3}=9.7GeV^3$$
(21)
Eq. (20) can be rearranged to give:

$$A^{3\pi }(s,t,u)=\frac{\lambda }{3}\{[\frac{m_\rho ^2}{m_\rho ^2s}(1+\frac{c}{3c}\frac{s}{m_\rho ^2})]+[st]+[su]\}$$
(22)
Eq. (19) yields a decay width $`\mathrm{\Gamma }(\rho \pi \gamma )=36KeV`$ which is too small compared with the experimental data. Eq. (20) with $`c=1`$ yield, on the other hand, a width which is 9/4 times larger and appears to be in better agreement with the data .

To make the $`\rho `$ unstable, we could make a substitution $`m_\rho ^2/(m_\rho ^2s)`$ with $`\mathrm{\Omega }(s)`$ in Eq. (19) and Eq. (20), then both resulting equations would violate the phase theorem. Following the spirit of the form factor calculation, one could try to make this replacement in Eq. (22) but the resulting equation still violates the phase theorem when the P-wave amplitude is projected out from the full amplitude. How to get a correct result which satisfies the elastic unitarity relation or the phase theorem , is a lenghthy problem and will be dealt in a forthcoming publication . We would like to point out however the ambiguity of the corresponding integral equation, which is a generalisation of the MO equation, also exists but its nature is more complicated than that given in Eq. (6).

We have shown in this note two methods to get agreement with experimental data. The first one consists in making a twice subtractions in the MO integral equation in order to suppress the high energy contribution in the dispersion integral and hence it allow us to use the approximation of the elastic unitarity for the pion form factor. The second method is to take into account of the inelastic effect which, around the $`\rho `$ region, can be simulated as the polynomial ambiguity of the solution of the MO equation. In both methods, the extra parameter is determined by the (low energy) rms radius of the pion and hence we can only calculate the $`\rho `$ leptonic width and also its total width Eq. (9). (It should be noticed that the rms pion radius can be calculated by a rapidly converging sum rule in terms of the magnitude of the pion form factor factor and the P-wave $`\pi \pi `$ phase shifts using Eq. (2) with $`n=1`$ ).

Within the frame work of dispersion relation and unitarity which gives rise to the MO equation, there is no room for a constant contact term. Its existence could, however, be introduced as a large energy behavior of the form factor or some higher energy inelastic effect as explained in this note.

This work is completed at the Institute for Nuclear and Particle Studies at the KEK High Energy Accelerator Research Organization. The author would like to thank Professors M. Kobayashi and H. Sugawara for hospitality.

Figure Captions

Fig.1 : The phases of the function $`\mathrm{\Omega }(s)`$ or $`V(s)`$ in degrees are given as a function of $`s(GeV^2)`$. The experimental data are taken from references .

Fig.2 :The square of the modulus of the function $`\mathrm{\Omega }(s)`$ is given as a function of $`s`$ in $`GeV^2`$ (dashed line). The square of the modulus of the pion form factor $`V(s)`$ as given by Eq. (13) (solid line). Experimental data are taken from references .


format: text

text:
# Implementation of the Five Qubit Error Correction Benchmark

## Abstract

The smallest quantum code that can correct all one-qubit errors is based on five qubits. We experimentally implemented the encoding, decoding and error-correction quantum networks using nuclear magnetic resonance on a five spin subsystem of labeled crotonic acid. The ability to correct each error was verified by tomography of the process. The use of error-correction for benchmarking quantum networks is discussed, and we infer that the fidelity achieved in our experiment is sufficient for preserving entanglement.

Robust quantum computation requires that information be encoded to enable removal of errors unavoidably introduced by noise Shor (1995); Steane (1996); Knill and Laflamme (1997); Shor (1996); Aharonov and Ben-Or (1996); Kitaev (1997); Knill et al. (1998); Preskill (1998). Thus, every currently envisaged scalable quantum computer has encoding, decoding and error-correction procedures among its most frequently used subroutines. It is therefore critical to verify the ability to implement these procedures with sufficient fidelity. The experimental fidelities achieved serve as useful benchmarks to compare different device technologies and to determine to what extent scalability can be claimed.

Liquid state nuclear magnetic resonance (NMR) is currently the only technology that can be used to investigate the dynamics of more than four qubits Cory et al. (1997); Gershenfeld and Chuang (1997). Although it is not practical to apply it to more than about ten qubits Warren (1997), it can be used to investigate the behavior of quantum networks on representative physical systems to learn more about and begin to solve the problems that will be encountered in future, more scalable device technologies. In this Letter, we describe an experimental implementation using NMR of a procedure for benchmarking the one-error correcting five-qubit code. This is the shortest code that can protect against depolarizing one-qubit errors Bennett et al. (1996); Laflamme et al. (1996). The experiment is one of the most complex quantum computations implemented so far and the first to examine the behavior of a quantum error-correcting code which protects against all one-qubit quantum noise. We discuss the principles underlying error-correction benchmarks and offer a sequence of specific challenges to be met by this and future experiments. The ultimate challenge is to demonstrate the ability to reduce the destructive effects of independent depolarizing errors. This seems to be out of reach of liquid state NMR experiments. Our experiment shows an average polarization preservation of $`67\%`$ corresponding to an entanglement fidelity of $`.75`$. This easily achieves the goal of demonstrating the preservation of entanglement in principle.

The five-qubit code. A quantum error-correcting code for encoding a qubit is a two-dimensional subspace of the state space of a quantum system. In the case of interest, the quantum system consists of five qubits. The code, $`C_5`$, can be specified as one of the 16 two-dimensional joint eigenspaces of the four commuting operators

$`\sigma _z^{(2)}\sigma _y^{(3)}\sigma _y^{(4)}\sigma _x^{(5)}`$ , $`\sigma _z^{(1)}\sigma _y^{(2)}\sigma _y^{(3)}\sigma _x^{(4)},`$
$`\sigma _y^{(2)}\sigma _z^{(3)}\sigma _z^{(4)}\sigma _z^{(5)}`$ , $`\sigma _x^{(1)}\sigma _z^{(2)}\sigma _x^{(3)}\sigma _z^{(4)}.`$ (1)
Here, $`\sigma _x^{(k)},\sigma _y^{(k)},\sigma _z^{(k)}`$ are the Pauli spin operators acting on qubit $`k`$. This is an instance of a stabilizer code Gottesman (1996); Calderbank et al. (1997), and the fact that it can be used to correct any one-qubit error is due to the property that every product of one or two Pauli operators acting on different qubits anticommutes with at least one of the operators in Eq. (1).

A typical application of a quantum error-correcting code is to protect a qubit’s state in a noisy quantum memory. As implemented in our experiment, the procedure begins with qubit $`2`$ containing the state to be protected and *syndrome* qubits $`1,3,4,5`$ in the initial state $`|1`$. A unitary *encoding* transformation is applied to map the two-dimensional input state space to the subspace of the code $`C_5`$. In the application, the five qubits are then stored in the noisy memory. In our experiment, we explicitly applied one of the correctable error. The information is retrieved by first applying the inverse of the encoding transformation to *decode* the state. The properties of the code guarantee that if the error was a Pauli operator acting on one of the qubits, which one occurred is reflected in the state of the syndrome qubits. To complete the process, conditional on the syndrome qubits’ state, it is necessary to *correct* the state of qubit $`2`$ by applying a Pauli operator. Quantum networks for encoding, decoding and error-correction are shown in Fig.1.

Benchmarking quantum codes.

The purpose of a benchmark is to compare the performance of different devices on the same task. Since quantum codes will be used to maintain information in future quantum computing devices, they are excellent candidates for benchmarking the reliability of proposed quantum processors. A basic quantum code benchmark consists of a sequence of steps that implement encoding, evolution, decoding and error-correction networks. In the simplest cases, the steps are applied to one qubit’s state using several ancilla qubits for the intermediate steps. An experimental implementation measures the reliability (see below) with which the qubit’s state is processed. It is necessary to include a means for verifying that a code with the desired properties was indeed implemented. For error-correcting and for stabilizer codes, this can be done by inserting $`180^{}`$ pulses applied to individual qubits in the evolution step and observing their effect on the output. In this case, the verification relies on the assumption that such pulses can be reliably implemented, a property that needs to be independently verified.

To allow for unbiased comparison of devices, the reliability measurement and the verification steps of the benchmark need to be standardized. There are many different ways of quantifying reliability. The best known such measure, *fidelity*, is based on the geometry of the state space. If the input state is $`|\psi `$ and the output density matrix is $`\rho `$, then the fidelity of the output is given by $`F(|\psi ,\rho )=\psi |\rho |\psi `$. This can be seen to be the probability of measuring $`|\psi `$ in a measurement which distinguishes this state from the orthogonal states. In our case we are interested in an arbitrary state of one qubit. One quantity of interest would be the worst case pure state fidelity, which minimizes the fidelity of the implemented benchmark over pure input states. However, an easier to use quantity is the *entanglement fidelity* $`F_e`$ Schumacher (1996), the fidelity with which a Bell state on the qubit and a perfect reference qubit is preserved. Entanglement fidelity does not depend on the choice of Bell state and has the property that $`F_e=1`$ if and only if the process perfectly preserves every input state. To avoid experimentally implementing the reference qubit, $`F_e`$ can be determined from the fidelities of pure states. Define $`|\pm =(|0\pm |1)/\sqrt{2}`$ and $`|\pm i=(|0\pm i|1)/\sqrt{2}`$ (the eigenstates of $`\sigma _x`$ and $`\sigma _y`$, respectively). Let $`F_s`$ be the fidelity of the process for input $`|s`$. Then

$$F_e=(F_0+F_1+F_++F_{}+F_{+i}+F_i)/41/2.$$
(2)
We advocate the use of entanglement fidelity as the standard reliability measure to be given when describing the results of a quantum benchmark involving processing of a quantum bit.

The standard verification procedure for a quantum code benchmark needs to be such that sufficiently high fidelity cannot be achieved without having implemented a code with the desired properties. For codes defined as the common eigenspace of a commuting set of products of Pauli operators (stabilizer codes), it is in principle enough to verify that applying a product $`P`$ of Pauli operators during the evolution has the expected effect on the output, namely that it differs from the input by the application of a Pauli operator $`\sigma (P)`$ determined by the code and the applied product. A single fidelity measure may be obtained by applying $`\sigma (P)`$ to the output and by averaging the resulting entanglement fidelities over all $`P`$. To make this procedure experimentally feasible, one may randomize the choice of $`P`$ and use statistical methods to estimate the desired average.

For benchmarks involving a quantum error-correcting code, the emphasis is on having corrected the set of errors $``$ for which it was designed, and verification involves applying the errors in $``$ during the evolution and observing the extent to which they are indeed corrected. Ideally, the errors occur naturally in the course of evolution, and one would like to see that information is preserved better by encoding it. In order to investigate the code in a controlled way, it is easier to apply different errors explicitly and observe the effect on the reliability of the process. The experiment described here involves measuring the entanglement fidelity for each of the one-qubit Pauli operators applied during the evolution.

When implementing a benchmark, it is useful to have some goals in mind. Each goal implies the demonstration of a non-trivial result. For benchmarks involving codes designed to correct independent errors on qubits, we offer a sequence of four such goals depending on how well the implementation succeeds at protecting against various error models. Most involve comparing the entanglement fidelity for two situations involving a specific error process $`_i`$. In the first, the information is stored in any one of the qubits, giving an optimum $`F_{e,1}(_i)`$. In the second, the information is stored by using the implemented code, giving an experimentally determined $`F_{e,C}(_i)`$. Numerical goals for $`F_{e,C}(_i)`$ are given under the assumption that the implementation induces depolarizing noise. The goals are: 1. Improvement where $`_1`$ is depolarization of each qubit with some probability $`p`$. For $`C_5`$ under the assumption that the error behavior is identically depolarizing (quantified by $`F_{e,C_5}`$) for each possible Pauli-product error during the evolution, this requires $`F_{e,C_5}>0.97`$, giving an improvement when $`p=0.08713`$. See Fig. 2. 2. Improvement where $`_2`$ is the process that first randomly chooses a qubit and then depolarizes it. For our code this requires $`F_{e,C_5}(_2)>.85`$. 3. Preservation of some entanglement for $`_2`$. This requires $`F_{e,C}(_2)>.5`$ Bennett et al. (1996). 4. Improvement in the presence of the *demonic* error process $`_4`$ that, knowing the method for storing the qubit, chooses the worst possible one-qubit depolarizing error and applies it. In this case we need $`F_{e,C}(_4)>.25`$. One of the ultimate goals might be to demonstrate that the code can be implemented sufficiently well to permit preservation of information by means of concatenation.

Experimental implementation. We used 13C labeled trans-crotonic acid (Fig. 3) synthesized as described in Knill et al. (2000), but with deuterated acetone as a solvent. A standard 500 MHz NMR spectrometer (DRX-500 Bruker Instruments) with a triple resonance probe was used to run the experiments. The five spin $`1/2`$ systems used for the code are the methyl group (M), C1, C2, C3 and C4. The methods of Knill et al. (2000) were used to prepare the methyl group as an effective spin $`1/2`$ system and to initialize the labeled pseudo-pure state $`\mathrm{𝟏}\sigma _z\mathrm{𝟏𝟏𝟏𝟏𝟏}`$ on all active nuclei. Here, $`\mathrm{𝟏}=|11|`$ and the last two nuclei are H1 and H2. The selection method was based on gradients, and the pseudo-pure state was subjected to a “crusher” gradient. To absolutely guarantee the pseudo-pure state, more randomization is required (see Knill et al. (2000)) but we did not implement this. H1 and H2 were not used and were only affected by some hard pulse refocusings on the protons. The state of H1 and H2 (up or down) induces an effective frequency shift on the other nuclei depending on the coupling constants and was compensated for in phase calculations. To greatly reduce the effect of radio frequency (RF) inhomogeneity, we used the nutation based selection scheme of Knill et al. (2000), applied to both the proton and the carbon transmitters. The quantum networks of Fig. 1 were directly translated to pulse sequences, again using the methods described in Knill et al. (2000). The only significant use of manual intervention was to place the refocusing pulses. The evolution period between encoding and decoding was carefully isolated from both the preceding and the following pulses: It implemented the identity unitary operator, or one of the one qubit $`180^{}`$ rotations by refocusing the molecule’s internal Hamiltonian and applying an extra inversion or by shifting the phase by $`180^{}`$. The qubit’s output state appeared on C1 at the end of the experiment. The peak group associated with C1 was observed in each experiment. Spectra were analyzed by comparing the spectrum of the pseudo-pure state $`\mathrm{𝟏}\sigma _x\mathrm{𝟏𝟏𝟏𝟏𝟏}`$ to the output, using the knowledge of the peak positions and shapes to compute relative intensities and phases. *No* phase adjustment was made after phasing the pseudo-pure state spectrum. This was possible since the relative phase is precomputed by the pulse compiler and integrated into the acquisition.

We performed one experiment for each of the $`16`$ possible evolutions with one-qubit or no Pauli error, each of the three initial states $`\sigma _x`$, $`\sigma _y`$ or $`\sigma _z`$ on C1, and each of three observations (no pulse, $`90^{}`$ $`X`$ pulse, or $`90^{}`$ $`Y`$ pulse) on C1. This resulted in a total of 112 experiments, each of which was repeated sufficiently often to get better than $`8\%`$ error in the inferred state of C1. For each evolution $`E`$ and input $`\sigma _u`$, we determined the amount of signal $`P(e,u)`$ in the correct direction in the output relative to the input signal. This requires “tracing out” the other spins, which was done by adding the intensities of each peak in the C1 spectrum that is associated with the $`\mathrm{𝟏𝟏}`$ state on $`H1`$ and $`H2`$. (The spectrum of C1 resolves all its couplings.) Thus, except for noise, $`1P(E,u)1`$. Under the assumption that input $`\mathrm{𝟏}𝟙\mathrm{𝟏𝟏𝟏𝟏𝟏}`$ results in no observable signal, the entanglement fidelity for a given evolution $`e`$ is given by $`F_e(E)=(P(E,x)+P(E,y)+P(E,z)+1)/4`$. We did not verify the assumption in this experiment, but note that it has been verified in related experiments Cory et al. (1998), and could be enforced by modifying the process with random pairs of canceling $`180^{}`$ pulses before and after the implemented pulse sequence. This would also enforce the depolarizing noise model for the implementation while preserving the observed polarization.

Results. Typical spectra compared to the spectrum of the input pseudo-pure state are shown in Fig. 4. The relative polarization after the error-correction procedure in the correct output state varies between $`48\%`$ and $`87\%`$. The distribution is shown in Fig. 5. The inferred entanglement fidelity for goals $`2`$ and $`3`$ is $`F_{e,C}(_2)=.75`$, with an estimated error of less than $`.02`$. Thus we clearly met goal $`3`$.

The reduction in polarization is due to thermal relaxation, incompletely refocused couplings (part of the pulse compiler’s optimization trade-offs), pulse errors due to non-ideal implementation of $`180^{}`$ and $`90^{}`$ pulses, and RF inhomogeneity (less than $`2\%`$ after our selection procedure). Most of the error is explained by the known relaxation times, suggesting that this is what limits the fidelities that can be attained using liquid state NMR. We note that the estimated phase relaxation times of well over a second are high when compared to those of nuclei in other molecules used in NMR quantum information processing experiments to date.

Discussion. Benchmarking quantum devices for quantum information processing is crucial both for comparing different device technologies and for determining how much control over a device is achievable and how to best achieve it. Given the need for and difficulty of achieving robust quantum control, we advocate the use of quantum coding benchmarks to determine the fidelity of the implementation of standard, verifiable processes. Unlike the experimentally implemented versions (up to 5 qubits) of the popular quantum searching and order-finding algorithms Jones and Mosca (1998); Jones et al. (1998); Chuang et al. (1998a, b); Marx et al. (1999); Vandersypen et al. (2000), quantum codes offer a rich source of complex and verifiable quantum procedures required in currently envisioned quantum computer architectures. Liquid state NMR has been used to implement several interesting, small quantum codes Cory et al. (1998); Leung et al. (1999), and a seven qubit cat-state benchmark Knill et al. (2000). In this letter, we have given specific goals for benchmarks involving error-correction and, for the first time, implemented a one-qubit error-correcting quantum code. The fidelity achieved is well above what is needed to demonstrate preservation of entanglement in the presence of a non-trivial error-model. It is unlikely that the goal of showing improvement in error-correcting independent depolarizing errors is achievable with liquid state NMR. A device that achieves this challenging goal will be well on the way toward realizing robustly scalable quantum computation.

Acknowledgments. Supported by the Department of Energy (contract W-7405-ENG-36) and by the NSA. We thank the Stable Isotope Resource for providing equipment and support. We thank Ryszard Michalczyk, Brian MacDonald, Cliff Unkefer and David Cory for their help.


text:
# On particle production for high energy neutrino beams

## 1 Introduction

A renewed interest in the predictions of absolute flux, energy distribution and flavour composition of neutrino beams produced at accelerators is motivated by the long-baseline neutrino beams in construction at CERN (CNGS CNGS ) and FNAL (NuMI NUMI ) and derived from the decays in flight of mesons generated by 400 GeV/c and 120 GeV/c proton interactions on graphite, respectively.

The available experimental data on particle yields from high energy proton interactions on light nuclei targets, such as the ones used in the production of neutrino beams, are not extensive and one is often faced with the problem of extrapolating to different target materials and shapes and to different incident proton energies.

Motivated by the need to estimate the neutrino flux at the West Area Neutrino Facility at CERN WANF , a consistent set of measurements of particle yields on Beryllium targets of different lengths were performed at CERN by the NA20 Coll. NA20  and subsequently by the NA56/SPY Coll. SPY1 ; SPY2 , with 400 GeV/c and 450 GeV/c incident protons, respectively. Secondary particles ($`\pi ^+`$, $`\pi ^{}`$, $`K^+`$, $`K^{}`$, $`p`$, $`\overline{p}`$) were measured in a transverse momentum range from 0 to 600 MeV/c, with the NA20 measurements in the high momentum range (60 GeV/c $`<p<`$ 300 GeV/c) and with NA56/ SPY extending to the low momentum region (7 GeV/c $`<p<`$ 135 GeV/c), thus covering also the relevant kinematic range of secondary particles for the planned long-baseline (LBL) neutrino beams.

In this work, in order to make these measurements of general applicability, the measured single-particle yields on Beryllium have been converted to single-particle invariant cross sections, correcting for target efficiency and tertiary particle production. The resulting data has been empirically parameterized as a function of the transverse momentum ($`p_T`$) and the scaling variable $`x_R=E^{}/E_{max}^{}`$, defined as the ratio of the energy of the detected particle in the centre-of-momentum frame and the maximum energy kinematically available to the detected particle, with a formula based on general physical arguments. The choice of the variables used in the description of the invariant cross section is motivated by the approximate scaling behaviour of the invariant cross section, when expressed in terms of these variables. A satisfactory comparison of this empirical parameterization, based on data collected with 400 GeV/c and 450 GeV/c protons, with invariant cross section measurements with 100 GeV/c protons shows that such approach might be useful also for the prediction of the NuMI neutrino beam at FNAL.

Extension of the parameterization to target materials other than Beryllium is possible by the known dependence of the invariant cross section on the atomic number as shown in the measurements of Barton et al. Bar87 . This is particularly relevant for the planned long-baseline neutrino experiments at CERN and FNAL, both planning to use a graphite target.

Such parameterization of the invariant cross section for proton-nucleus interactions could be implemented in a full Monte Carlo simulation of the target geometry. On the other hand, we show that it is possible to compute analytically particle yields from a given target by taking into account the geometry of the target and the contribution of cascade processes. Based on the NA56/SPY measurements of particle yields from targets of different lengths, an empirical formula describing tertiary particle production is derived.

This approach, coupled to the tracking of the produced secondary hadrons through the neutrino beam-line elements (magnetic lenses, collimators and decay tunnel), has been implemented in a fast Monte Carlo for the explicit calculation of neutrino beam fluxes. Effects due to reinteractions in the beam-line material have been parameterized with good approximation. A comparison with the neutrino and anti-neutrino spectra measured by the CHARM II experiment at CERN is shown. Predictions for the planned CNGS long-baseline neutrino beam are given.

With respect to the standard approach using complex Monte Carlo simulations based on hadronic cascade codes, such as GEISHA GHEI  or FLUKA92 FL92 , as implemented in the GEANT package GEAN , CALOR95 CA95 , the stand-alone version of FLUKA FONE  or MARS MARS , a simple empirical parameterization of particle yields as described here has the advantage of a more transparent physical input to the calculations. Moreover, when coupled to the tracking of the secondary hadrons from the target, it becomes a fast neutrino Monte Carlo generator with the obvious advantage of enhanced statistics. Detailed studies of systematic effects due for example to horn misalignments in the neutrino beam-line are made possible even with limited computer time.

## 2 Parameterization of p-Be data

An empirical formula for particle production, based on the yields from thick Be targets measured by NA20, has long ago been proposed by Malensek Mal81 . This has already offered an alternative approach to neutrino beam computations NTeV . Extensions of this parameterization to include effects related to the target geometry have been proposed Bon96 . Malensek’s parameterization, however, fails to describe the particle production in the low momentum region covered by NA56/SPY measurements and cannot be stretched to consistently fit the whole set of available data on p-Be interactions. An alternative parameterization of NA20 data on particle yields, proposed by NA20 NA20 , gives an even less accurate description of the NA56/ SPY data.

At variance with these works, the particle yields measured by the NA20 and NA56/SPY collaborations have been converted to inclusive invariant cross sections for particle production, which have then been parameterized.

### 2.1 The data

In Ref. SPY2 , inclusive invariant cross sections for particle production in the forward direction in p-Be interactions have been reported by the NA56/SPY Collaboration. We have used the same method discussed in that paper to derive inclusive invariant cross section at all energies and angles from the published NA20 NA20  and NA56/SPY SPY1 ; SPY2  data on particle yields in p-Be interactions.

The inclusive invariant particle production cross section is related to the measured yields from targets of finite length by:

$$E\times \frac{d^3\sigma }{dp^3}=(100Y)\frac{E}{p^3}\frac{A}{N_0\rho \lambda _pF(L)}$$
(1)
where $`Y`$ is the yield per incident proton per sr ($`\mathrm{\Delta }p/p`$ %); $`A`$, $`\rho `$, $`L`$ the atomic mass number, density and length of the target; $`N_0`$ the Avogadro’s number; $`E`$ and $`p`$ the energy and momentum of the secondary particle, respectively. The target efficiency $`F(L)`$ can be estimated assuming that the produced secondary particles if absorbed in the target do not generate additional particles. Under this approximation and for forward production, one has

$$F(L)=_0^L\frac{dz}{\lambda _p}e^{z/\lambda _p}e^{(Lz)/\lambda _s}$$
(2)
representing the convolution of the probability that the primary proton interacts between $`z`$ and $`z+dz`$ in the target and that the produced secondary particle escapes from the target ($`\lambda _p`$ and $`\lambda _s`$ are the effective mean free paths of primary and secondary particles respectively). This model neglects the production of tertiaries. This effect can be accounted for, in a model independent way, by estimating the invariant cross sections for various target lengths and then extrapolating the results to zero target thickness SPY2 .

This procedure also allows to estimate at each momentum the contribution to the total yield from finite targets due to tertiary particles.

For particle production at angles different from zero, the method is easily generalised by taking into account the target geometry in eq. (2). At angles different from zero, however, the NA56/SPY collaboration has measured the production yields only at a fixed target length. In order to derive inclusive invariant cross sections from those data points, we have made the additional assumption that the correction for tertiary production be angle independent. This has been verified to be approximately true in the kinematic region covered by NA20 data, collected with a planar target setup identical to the one adopted in the NA56/SPY experiment.

In this way we have derived the single-particle inclusive invariant cross sections at all available energies and production angles from NA20 and NA56/SPY data. Table 1 summarizes our results. Errors include the experimental errors and the uncertainty on the extrapolation procedure.

### 2.2 The parameterization of inclusive invariant cross sections

One of the goals of this analysis is to adopt a simple functional form for inclusive particle production which will be appropriate for extrapolation to different centre of mass energies and/or secondary particle momenta.

Feynman has speculated in very general terms about the shape and energy dependence of inclusive processes Fey69 , suggesting an approximate scaling behaviour of the single-particle inclusive invariant cross section when expressed in terms of the transverse momentum ($`p_T`$) and the longitudinal variable $`x_F=2p_L^{}/\sqrt{s}`$, where $`p_L^{}`$ and $`\sqrt{s}`$ are the longitudinal momentum of the detected particle and the total energy in the centre-of-momentum frame. A factorization in $`x_F`$ and $`p_T`$ of the invariant cross section has also been advocated as an experimental fact.

An alternative scaling variable $`x_R=E^{}/E_{max}^{}`$, defined as the ratio of the energy of the detected particle in the centre-of-momentum frame and the maximum energy kinematically available to the detected particle, was suggested by Yen Yen74  and Taylor et al. Tay76  and shown to greatly extend the range of validity of scaling at sub-asymptotic energies. As discussed in section 4.2, in the comparison of NA56/SPY and NA20 data to $`pBe`$ data collected at lower energies, we have not observed an improved scale invariance when $`x_R`$ is used. Nonetheless, we have adopted this variable because of two practical advantages, which translates in a simpler analytical parameterization: $`x_R`$ is always positive by construction and can never be zero, except for infinite energy in the centre-of-momentum frame.

Early phenomenological analyses of $`pp`$ data, based on the constituent quark model, helped clarify that in hadronic processes the produced particles reflect the motion of the constituents Och77 . In particular, in the fragmentation region of the projectile, the longitudinal momentum distribution of produced hadrons reflects the momentum distribution of valence quarks inside the incident hadron. This suggests an $`x`$ dependence of hadron production at large $`x`$ 111At large $`x_F`$, $`x_F`$ and $`x_R`$ are equivalent. with a functional form similar to the one describing parton distributions. The functional shape of non-direct formation processes at small $`x`$, in which sea quarks are involved, however, is not easily described in this framework.

After some trials and considerations of the physical process, the following empirical parameterization of the inclusive invariant cross sections for positive sign secondary meson ($`\pi ^+`$, $`K^+`$) production in p-Be interactions has been found to give a suitable description of data:

$`(E\times {\displaystyle \frac{d^3\sigma }{dp^3}})`$ $`=`$ $`A(1x_R)^\alpha (1+Bx_R)x_R^\beta \times `$ (3)
$`(1+a^{}(x_R)p_T+b^{}(x_R)p_T^2)e^{a^{}(x_R)p_T}`$
where $`a^{}(x_R)=a/x_R^\gamma `$ and $`b^{}(x_R)=a^2/2x_R^\delta `$.

These formula assumes an approximate factorized scaling form in $`x`$ and $`p_T`$. The $`(1x)^\alpha `$ behaviour at large $`x`$ is theoretically motivated on the basis of quark counting rules Bla74 ; Gun79 . The $`x^\beta `$ behaviour empirically accounts for the non-direct hadron formation mechanism at small $`x`$.

The $`p_T`$ behaviour is modelled with the known exponential fall of soft interactions and a polynomial behaviour to interpolate the low $`p_T`$ part of the spectrum. The $`x`$ dependence of $`a^{}(x)`$ and $`b^{}(x)`$ is introduced to parameterize the violation of $`p_T`$ invariance observed in the data. Models based on the parton structure of the hadrons predict a $`p_T^n`$ dependence of the cross section at large $`p_T`$, where hard scattering processes take over. A reasonable parameterization of p-Be data with this functional form has not been found. The possibility that the proposed parametric form fail to describe particle production in the $`p_T`$ region not covered by present data is acknowledged.

The ratio r of positive to negative data ($`\pi ^+/\pi ^{}`$ or $`K^+/K^{}`$) has been instead parameterized with the empirical formulae:

$`r(\pi )=r_0(1+x_R)^{r_1}`$ (4)
$`r(K)=r_0(1x_R)^{r_1}`$ (5)
The shape of these ratios is supported by the phenomenological analysis of $`pp`$ data of Ref. Och77 , showing that $`r(\pi )1`$ for $`x0`$ and rises to about 5 for $`x1`$, closely following the $`u/d`$ ratio of valence quarks in the projectile proton, while $`r(K)`$ has a $`(1x)^3`$ behaviour for $`x1`$. NA56/SPY and NA20 data only cover the fragmentation region of the proton at large $`x`$ and the central region. At large $`x`$ a functional behaviour similar to the one exhibited by $`pp`$ data is expected.

In order to keep the number of free parameters limited, positive and negative mesons are assumed to have the same $`p_T`$ distributions. This has long been known to be only approximate in $`pp`$ data Allab .

Table 2 summarizes the results of our best-fit to the data. As indicated in the table, some of the parameters have been fixed in the fitting procedure, since they appeared to be redundant. In the $`K^\pm `$ fit, $`\delta =2\gamma `$ was chosen and $`B`$ was set to zero, since its fitted value was found to be consistent with zero within errors.

The comparison between the empirical parameterization and the experimental data is shown in figure 1. The accuracy of the parameterization of the $`\pi ^\pm `$ and $`K^\pm `$ data is displayed in figure 2, showing the relative discrepancy between our parameterization and the experimental data as a function of $`x_R`$. The proposed parameterization gives an accurate fit of $`K^\pm `$ data with a reduced $`\chi ^20.85`$, while the reduced $`\chi ^2`$ is somewhat larger ($`\chi ^2/N_{dof}=77.1/37`$) for $`\pi ^\pm `$ data. This partly reflects our difficulty to parameterize the $`x`$ dependence of the $`p_T`$ distribution, although about 1/3 of the $`\chi ^2`$ is contributed by the two data points (one measured by NA20 and the other by NA56/SPY) for positive pion production in the forward direction at $`x=0.3`$, which are about 30-40% off the best-fit prediction. A reduced $`\chi ^2`$ around 1 is obtained, if a relative error of 10% on each experimental point is assumed. We conclude that the proposed empirical formulae are adequate to describe NA20 and NA56/SPY data with a 10% accuracy.

Proton and anti-proton data, which are of less direct interest for neutrino beams, have been parameterized with the following empirical formulae:

$`(E\times {\displaystyle \frac{d^3\sigma }{dp^3}})_{pBepX}`$ $`=`$ $`A(1+Bx_R)(1x_R)^{bp_T^2}\times `$ (6)
$`(1+ap_T+{\displaystyle \frac{a^2}{2}}p_T^2)e^{ap_T}`$
$`(E\times {\displaystyle \frac{d^3\sigma }{dp^3}})_{pBe\overline{p}X}`$ $`=`$ $`A(1x_R)^\alpha x_R^\beta \times `$ (7)
$`(1+ap_T+{\displaystyle \frac{a^2}{2}}p_T^2)e^{ap_T}`$
For anti-protons a functional similar to the one given in (3) has been adopted, except that an exact factorization in $`x_R`$ and $`p_T`$ has been assumed, since this was sufficient to give a reasonable fit to data. For protons the “leading particle effect” had to be taken into account. A reasonable fit to data has been obtained by following the empirical observation that the longitudinal momentum distributions of the leading nucleon in $`pp`$ collisions is flat Coc71 , which translates into a linear rise of the inclusive invariant cross section as a function of $`x_R`$. The transverse momentum distribution is also affected by the leading particle effect, resulting in an enhanced leading particle production in the forward direction (see for example Tay76 ; Allab ). In the proton fit, this is empirically accounted for by the term $`(1x_R)^{bp_T^2}`$. The parameter $`a`$, that controls the shape of the $`p_T`$ distribution for non-leading particle production, has been assumed to be the same for protons and anti-protons.

Results of these fits are also given in Table 2 and the comparison between the empirical parameterization and the experimental data is shown in figure 3. Our parameterization gives a satisfactory description of proton and anti-proton inclusive production, with a reduced $`\chi ^2`$ about 1, in the range covered by NA56/SPY and NA20 data. At larger values of $`x`$, proton production by means of diffraction should occur. This is not described by our parameterization.

## 3 Production of neutral Kaons

The knowledge of neutral kaon production in p-Be interactions is important for the exact calculation of the $`\nu _e`$ background in neutrino beams coming from $`K_L^0`$ decays. As an example, in the case of the WANF beam serving the CHORUS and NOMAD experiments, the estimated contribution to the $`\nu _e`$ content of the beam due to $`K_L^0`$ is around 15%. In addition, the knowledge of the neutral kaon production is also of interest for the experiments on CP-violation using neutral kaon beams and it is particularly relevant in searches for rare neutral kaon decays.

A rather complete set of data on $`K_S^0`$ production at angles different from zero was collected by Skubic et al. Sku78 , using a 300 GeV/c proton beam hitting Be targets. Measurements of $`K_S^0`$ production in the forward direction from 200 GeV/c protons on beryllium were performed by Edwards et al. Edw78 . Both these experiments, however, covered a region of secondary momenta corresponding to $`x_F>0.25`$, which only partially overlaps to the momentum region important for present and future neutrino beams.

Additional information on neutral kaon production can be obtained starting from the available measurements of charged kaon production in p-Be interactions. Charged and neutral kaon production rates can be related in the frame of a simple quark parton model, in which valence quarks ($`q_v`$) and sea quarks ($`q_s`$) are considered.

In the simplest model, isospin symmetry may be assumed, giving $`u_v/d_v=2`$, $`u_s=\overline{u}_s=d_s=\overline{d}_s`$ and $`s_s=\overline{s}_s`$, from which one finds that the production of charged and neutral kaons should be related by

$$N(K_S^0)=N(K_L^0)=\frac{1}{4}(N_{K^+}+3N_K^{})$$
(8)
This model has been shown to provide an accurate description of the asymmetry of $`K^0`$/$`\overline{K^0}`$ production in the range $`0.18<x_F<0.36`$ in neutral kaons beam-lines Dob90 . We have verified that it also agrees within 15% with direct measurements of $`K_s`$ production up to around $`x_F=0.5`$. At larger values of $`x_F`$, i.e. in the fragmentation region of the proton, we find that a reasonable description of existing $`K_s`$ production data requires a “dynamic” model in which the $`u_v/d_v`$ ratio is let to be $`x`$-dependent and estimated from the $`\pi ^+/\pi ^{}`$ production ratio, following the arguments given by Ochs Och77 , who has demonstrated a remarkable empirical similarity between the $`\pi ^+/\pi ^{}`$ production ratio in $`pp`$ collisions and the $`u/d`$ quark ratio measured in deep-inelastic lepton-proton scattering. The exploitation of this model is beyond the aim of this paper222A similar “dynamic” prediction of neutral kaon production was already considered in ref. Bor96 . These models however seems less accurate than formula (8) in the description of the $`K^0`$/$`\overline{K^0}`$ production asymmetry at small/intermediate $`x_F`$. and neutral kaon production is always estimated from eq. (8) in the following.

## 4 Particle production for neutrino beams

### 4.1 Scaling to targets of different materials

Beryllium targets have been used in most neutrino beams derived from proton beams extracted from the Super Proton Synchrotron (SPS) at CERN, but the use of graphite as target material seems promising in view of operation with short spills in fast extracted proton beams. In particular, the use of a graphite target is foreseen in the design of both the CNGS and the NuMI beams CNGS ; NUMI . Prescriptions to rescale the inclusive invariant cross-sections to different target materials are given in the following.

Invariant cross sections $`E\frac{d^3\sigma ^{hA}}{dp^3}`$ for hadron-nucleus interactions ($`hAh^{}X`$) depend on the mass number $`A`$ of the target nucleus, via parameterizations of the type:

$$E\frac{d^3\sigma ^{hA_1}}{dp^3}=(\frac{A_1}{A_2})^\alpha E\frac{d^3\sigma ^{hA_2}}{dp^3}$$
(9)
where a value for $`\alpha =2/3`$ would correspond to the case where particle production off a nucleus is identical to the production off a single nucleon per inelastic collision. In accordance with the scaling hypothesis, $`\alpha `$ has been found to be weakly dependent on the incident beam momentum. It depends on the incident hadron type $`h`$ and it is a smooth function of $`p_T`$ and $`x_F`$ of the produced hadron. It has moreover been experimentally observed that, in first approximation, it is independent of the detected particle type, with perhaps the exception of anti-protons Bar87 .

A parameterization of $`\alpha `$ as a function of $`x_F`$ has been proposed by Barton et al. on the basis of their and previous data at $`p_T=0.3`$ GeV/c Bar87 . A $`p_T`$ dependence of $`\alpha `$ has been clearly observed in $`K_S^0`$ and $`\mathrm{\Lambda }`$ production by Skubic et al. Sku78 . A suitable representation of the whole set of data can be obtained with the parameterization:

$$\alpha (x_F)=(0.740.55x_F+0.26x_F^2)(0.98+0.21p_T^2)$$
(10)
where the $`x_F`$ dependence is taken from the fit of Barton et al. at $`p_T=0.3`$ GeV/c and the $`p_T`$ dependence is fitted to Skubic et al. data and normalized in such a way that it reduces to the parameterization of Barton et al. at $`p_T=0.3`$ GeV/c.

In conclusion, a conservative estimate of the uncertainty in the extrapolation from beryllium to carbon data sits (in the $`p_T`$ range of interest for neutrino beams: up to $`600`$ MeV/c) around 5%, on the top of a measurement error of 5–10%, depending on the secondary momentum, for the cross sections on beryllium. The estimate of this systematic uncertainty is based both on data collected by Barton et al. Bar87  and on the extensive compilation of J. Kuhn on nuclear dependence for $`pA\pi ^{}X`$ interactions Kuh78 .

### 4.2 Scaling to different centre of mass energies

The NuMI neutrino beam at FNAL is planned to be derived from a primary proton beam of 120 GeV/c momentum, resulting in a centre-of-mass energy about two times smaller than that available at NA56/SPY and NA20. Besides its general interest, a test of the scaling hypothesis of one-particle inclusive invariant cross-sections is thus relevant to assess the range of validity of the proposed parameterization. This has been studied by comparing our parameterisation to available $`pA`$ data collected at different centre-of-mass energies.

In figure 5, our prediction for one-particle inclusive invariant cross section are compared to data collected with 100 GeV/c protons on carbon target Bar87 , where the extrapolation from Be to C has been made using formulae (9) and (10). As discussed above $`x_R=E^{}/E_{max}^{}`$ has been used as scaling variable. The agreement between pion data shown in the figure and our empirical parameterization is excellent up to about $`x0.8`$. A good agreement is also found when kaon data are considered, although the precision of kaon data from ref. Bar87  is poorer.

Figure 5 shows a comparison of our parameterization with one-particle inclusive production data measured with 24 GeV/c protons hitting a Beryllium target Eic72 . In that work, Lorentz invariant particle densities $`\omega (p_{Lab},\theta )`$ were measured as a function of the particle momentum $`p_{Lab}`$ and production angle $`\theta `$. These data have been converted into invariant cross sections using the relation

$$E\times \frac{d^3\sigma }{dp^3}=\frac{1}{2}\sigma _{abs}\omega (p_{Lab},\theta )$$
(11)
where $`\sigma _{abs}`$ is the absorption cross section for $`pBe`$ Car79 . A reasonable agreement is observed in the shape of the distributions for $`\pi ^+`$ and $`K^+`$, although the estimated production of $`\pi ^+`$ is about 35 $`\pm `$ 15% lower than that measured. This is also true for negative pions, not displayed in the figure, while the agreement is somewhat worse for the other particles.

Given its general interest, we have tested the scaling hypothesis on the same data also using $`x_F`$ and $`x_{Lab}=p/p_{inc}`$, the latter defined as the ratio of the momentum of the detected particle in the laboratory reference frame to the momentum of the incoming proton, as longitudinal variables. A better agreement with $`\pi ^\pm `$ data is obtained, if $`x_{Lab}`$ is used, which however shows a worse agreement to kaon data333All the variables are off roughly by the same amount in the description of the $`K/\pi `$ production ratio.. The variable $`\mathrm{\Delta }y=y_{max}y`$, where $`y=\frac{1}{2}\mathrm{log}((E+p_L)/(Ep_L))`$ is the rapidity of the produced particle and $`y_{max}`$ is the maximum rapidity kinematically available to that particle, has also been tried. At variance with the previuos longitudinal variables, $`\mathrm{\Delta }y`$ is Lorentz invariant and scales the phase space at different centre-of-mass energies independently of the reference frame. Although we find a somewhat improved agreement to data in this case, scaling to data collected with 24 GeV/c protons is only approximate even with this variable.

Still, in the centre-of-mass energy range of interest for present high-energy neutrino beams, the agreement between our parameterization scaled according to $`x_R`$ and the data is satisfactory.

### 4.3 Yields from finite length targets

With reference to eq. (1), the differential particle production along the target can be parameterized as:

$$\frac{dY(E,p_T,z)}{dz}=\frac{N_0\rho \lambda _p}{100A}\frac{p^3}{E}\left(E\times \frac{d^3\sigma }{dp^3}\right)f(z)$$
(12)
where $`f(z)dz`$ is the probability that the outgoing particle be produced at a depth $`z`$ to $`z+dz`$ inside the target and the other quantities have been introduced in Section 2.1.

In general $`f(z)`$ will also depend on the production angle of the secondary particle Bon96 . In the naive reabsorption model introduced in Section 2.1, one has:

$$f(z,\theta )=\frac{1}{\lambda _p}e^{z/\lambda _p}e^{z^{}(z,\theta )/\lambda _s}$$
(13)
which can be obtained by differentiating the target efficiency $`F(L)`$ of equation (2), with $`z^{}(z,\theta )`$ representing the target length that the produced hadron has to cross to escape from the target ($`Lz`$ for forward production). This naive model is only adequate to give a first order description of particle production from relatively thin targets or at large $`x`$, but in general particle production by means of cascade processes cannot be neglected.

Since all types of hadrons can be produced when an energetic hadron of any flavour interacts in the target, the production of particles from a target of finite length is governed by a set of coupled transport equations that depend on the properties of the particles, their interactions and on the structure (geometry and material) of the target. Given their complexity, the solution of these equations is in general addressed numerically with hadron cascade/transport codes, calibrated on a large set of hadron production data. However, in the spirit of this work, the main features of the cascade process in needle-shaped targets, like those used for neutrino beams, can be described in a simple parametric form and tuned on NA56/SPY data on yields from targets of different length and geometry.

In hadronic interactions, about half of the available energy is typically dissipated in multiple particle production (inelasticity $`k0.5`$) with average transverse momentum $`p_T0.35`$ GeV/c . The remainder of the energy is carried by fast forward-going “leading” particles. Only these particles are responsible for the propagation of the cascade in thin and long targets, since they are almost collinear with the primary beam444The divergence of leading particles with respect to the primary beam is typically $`p_T/(1k)p_{inc}`$ few mrad for proton energies of order a few hundred GeV. and their energy is large enough to result in a sizeable yield of additional energetic particles.

In the collinear approximation and assuming that leading particles in $`pA`$ interactions are mostly protons, thus characterized by an effective mean free path $`\lambda _p`$ equal to that of primary protons, the naive reabsorption model of equation (13) can be improved by the expansion:

$$f(z)=\frac{1}{\lambda _p}e^{z/\lambda _p}[1+A_h(x)\frac{z}{\lambda _p}]e^{z^{}(z,\theta )/\lambda _s}$$
(14)
where the second term in brackets accounts for particle production by reinteractions of secondary particles (tertiary particle production), while higher rank contributions are neglected, since neutrino targets are typically of order $`2\lambda _p`$ and the less energetic is the reinteracting particle the lower is the yield of produced particles.

In equation (14), $`A_h(x)`$ weights the probability that reinteractions of secondary particles will result in a produced hadron $`h`$ of fractional longitudinal momentum $`x`$ and is given by:

$$A_h(x)=\underset{h^{}}{}_x^1𝑑x^{}G_{h^{}h}(x,x^{})\mathrm{\Phi }_h^{}(x^{})$$

where $`G_{h^{}h}(x,x^{})`$ is the $`p_T`$ integrated cross section for inclusive production of $`h`$ in $`h^{}`$ interactions and $`\mathrm{\Phi }_h^{}(x^{})`$ is the flux of $`h^{}`$ particles. Thus, $`A_h(x)`$ is a function of $`x`$ and $`h`$ only. Its dependence on these variables has been derived from NA56/SPY measurements of inclusive yields in the forward direction from finite length targets.

In the short target approximation, valid until the target length is $`L<<\lambda _p\lambda _s/(\lambda _s\lambda _p)`$, integration of equation (14) in the forward direction (i.e. $`z^{}=Lz`$) predicts the fraction $`t(L)`$ of tertiary particle production to be approximately a linear function of $`L`$ given by:

$$t(L)A_h(x)\frac{L}{2\lambda _p}.$$
(15)
This linear behaviour is in agreement with experimental data, as shown figure 6, where the fraction of tertiary production as a function of target length for $`\pi ^\pm `$ and $`K^\pm `$ as derived from NA56/SPY data is displayed, and was exploited in section 2.1 to derive inclusive invariant cross sections.

Compared to the naive absorption model, the excess of particle production present in the NA56/SPY data increases with decreasing momentum. Figure 7 shows the fraction of tertiary particles produced in the forward direction for a 100 mm target as a function of the fractional longitudinal momentum $`x_{Lab}=p_{Lab}/p_{inc}`$. The curves shown in the figure represent the best-fit of $`A_h(x)`$ to data according to the empirical parameterization:

$$A_h(x_{Lab})=A_h^{ter}(1x_{Lab})^{b_h^{ter}}$$
(16)
where a common value of $`A_h^{ter}`$ and $`b_h^{ter}`$ has been chosen for positive and negative particles. The best-fit results are reported in table 3.

In this simple model, no $`p_T`$ dependence of tertiary particle production is predicted. Indeed, assuming Feynman scaling and the (approximate) factorization in $`x`$ and $`p_T`$ of the invariant cross-section, the $`p_T`$ distribution of tertiary particles is expected to be the same as that of secondary particles. It follows that the divergence of leading particles with respect to the primary beam direction is in general negligible as compared to the typical production angle of lower momentum tertiary particles ( collinear approximation). These considerations are supported by NA56/SPY data that show no $`p_T`$ dependence of the transverse momentum distribution of tertiary particle production in a 300 mm Be needle-shaped target once mere geometric effects are considered (see figure 20 of Ref. SPY2 ).

Equation (14) is valid for NA56/SPY and NA20 targets, but fails to give the correct prediction for tertiary particle production when the target length is of order $`\lambda _{geom}Rp_l/p_T`$, where $`R`$ is the transverse dimension (radius) of the target and $`p_l=(1k)p_{inc}`$ is the typical momentum of the leading particle. In this limit, there is a sizeable probability that some of the leading particles also escape from the side of the target before interacting. This can be accounted for, modifying equation (14) into:

$$f(z)=\frac{e^{z/\lambda _p}}{\lambda _p}\left[1+A_h(x)_0^z\frac{dy}{\lambda _p}w_z(y)\right]e^{z^{}(z,\theta )/\lambda _s}$$
(17)
where $`w_z(y)`$ is an acceptance factor, which in general depends on the $`p_T`$ and energy distributions of the leading particles and on the target geometry555In equation (17), the collinear approximation has been again advocated to factorize the acceptance factor and the reinteraction probability of leading particles.. For cylindrical symmetry, one can write:

$$w_z(y)=\frac{1}{p_T}_0^{p_T^{max}}h(p_T)p_T𝑑p_T$$
(18)
where $`p_T^{max}p_lR/(zy)`$ in the collinear approximation and $`h(p_T)`$ is the $`p_T`$ distribution of the produced particles.

In the simulation of neutrino beams, discussed in the following, the computation of this acceptance factor has been performed under the approximation that leading particles carry on average half of the primary beam momentum and that their transverse momentum distribution is represented by a pure exponential of average transverse momentum $`p_T`$. For cylindrical symmetry this allows the analytical integration of equations (17) and (18). In case of continuous targets of radius $`R`$, the integration yields:

$$f(z)=\frac{e^{z/\lambda _p}}{\lambda _p}\left[1+A_h(x)\frac{z}{\lambda _p}(1e^{\lambda _{geom}/z})\right]e^{z^{}(z,\theta )/\lambda _s},$$
(19)
which reduces to equation (14), predicting a linear increase of the differential production of tertiary particles along the target, until $`z<<\lambda _{geom}`$ (short target approximation). At larger depth inside the target, the differential production of tertiary particles will tend to saturate to a value proportional to $`\lambda _{geom}/\lambda _p`$.

In case of a segmented target, made of several rods interleaved by air gaps, as that used in the WANF beam and considered for the forthcoming CNGS beam at CERN, the details of the target geometry have to be considered. The analytical integration of equations (17) and (18) is still possible under the same approximation as that adopted to derive equation (19) and it has been used in the evaluation of the efficiency of segmented targets.

In Figure 8 we present predictions for secondary meson production at target level, based on our parameterization of the inclusive invariant cross-sections corrected for target efficiency. Momentum spectra (top), angular distributions (middle) and exit point profiles (bottom) of the secondary mesons in the case of the CNGS target geometry and for an angular acceptance of 30 mrad (much larger than the horn aperture of the CNGS beam-line described in the next Section) are shown. For sake of comparison, the predictions of a FLUKA based full Monte Carlo program FONE ; CNGS  are also shown. The agreement is remarkable, provided that tertiary particle production is taken into account in the evaluation of the target efficiency, as discussed in this Section.

A slight disagreement in the $`p_T`$ distributions of all charged particles is visible mainly at large $`p_T`$, that translates in a slightly more forward peaked distribution of the exit points along the target. This disagreement appears in a kinematic region not covered by NA56/SPY and NA20 data. However, as remarked in the next Section, neutrino beam predictions are only slightly affected by this large $`p_T`$ region, since it falls outside the momentum and angular acceptance selected by the focusing optics. The systematic difference between the two models is thus negligible.

## 5 Simulation of neutrino beams

As it is well known, a “classic” wide band neutrino beam is produced from the decay of mesons, mostly $`\pi `$’s and $`K`$’s. Mesons are created by the interaction of a proton beam into a needle shaped target, they are sign-selected and focused in the forward direction by two large acceptance magnetic coaxial lenses, conventionally called at CERN horn and reflector, and finally they are let to decay into an evacuated tunnel pointing toward the detector position.

In case of positive charge selection, the beam content is mostly $`\nu _\mu `$ from the decay of $`\pi ^+`$ and $`K^+`$. Small contaminations of $`\overline{\nu }_\mu `$ (from the defocused $`\pi ^{}`$ and $`K^{}`$) and $`\nu _e`$ (from three-body decay of $`K`$’s and $`\mu `$’s) are present at the level of few percent.

As an example, a schematic layout of the future CERN to Gran Sasso neutrino beam, CNGS CNGS , is shown in Figure 9. Its main characteristics, relevant for beam simulation purposes, are listed in Table 4 together with those of the old West Area Neutrino Facility (WANF) at CERN in the configuration set up for the CHARM II experimentWANF .

The neutrino fluxes for such a kind of beam are relatively easy to predict once the secondary meson spectra are known, because the meson decay kinematics is well understood and the geometry of the decay tunnel is quite simple.

Uncertainty in the estimation of the neutrino fluxes could arise because secondary mesons are selected over a wide momentum range and over a wide angular acceptance ($`20`$ mrad).

Re-interactions of secondary mesons in the target and downstream material contribute to reduce the neutrino fluxes and increase the uncertainty in the calculations (mainly for the wrong sign and wrong flavour contaminations). They are generally minimized using a target made of a number thin rods of low Z material interleaved with empty spaces (to let the secondary mesons exit the target without traversing too much material). In addition the amount of material downstream of the target (i.e. horn and reflector conductor thickness) is kept to the minimum.

The parameterization of the secondary meson production from protons onto a thin target, proposed in this paper, is thus well suited to be used in neutrino beam simulations both because it extends its prediction over a wide range of longitudinal and transverse momenta and also because the small fraction of tertiary production from re-interactions in the target and downstream material can be accounted for with the approximations described in Section 4.3.

A comparison of the neutrino flux prediction based on the proposed parameterization with some measured spectra is thus an effective estimator of the quality of the secondary mesons parameterization. For this purpose, our parameterization has been coupled with a neutrino beam simulation program able to provide rapid and accurate predictions of neutrino spectra at any distance (i.e. short and long base line). The comparison has been performed both with already published data (CHARM II) and with predictions for the future CNGS long-baseline neutrino beam generated with GEANT and/or FLUKA based Monte Carlo programs.

### 5.1 The simulation program

The program is a stand-alone code developed as a tool that allows to vary and optimize all elements and the geometry (in 3-D) of the beam line providing the results in terms of neutrino spectra and distributions at large distance with high statistics and in short time. As a by-product, the program provides also spectra and distributions of secondary hadrons and muons along the beam-line and in the muon monitor pits after the hadron stop at the end of the decay tunnel.

The underlying idea is that in order to produce rapidly a neutrino spectrum at large distance over a small solid angle (typically $`d\mathrm{\Omega }10^{10}`$ for the future LBL beams of CERN and FNAL), one has to force all the mesons to decay emitting a neutrino, and force all neutrinos to cross the detector volume. A weight is then assigned to each neutrino, proportional to the probability that this process actually happened.

In practice this method is implemented by subdividing the simulation into four subsequent steps (as described in detail in the following). The weight assigned to each neutrino is the product of factors originated in each step times the solid angle subtended by the detector:

$$W_{tot}=\underset{i}{}W_id\mathrm{\Omega }_{det}.$$

#### 5.1.1 Mesons production along target

First the total number of secondary mesons ($`\pi ^+`$, $`\pi ^{}`$, $`K^+`$, $`K^{}`$ and $`K_L^0`$) to be generated is calculated by integration of the yield calculated in our model over the interesting range of longitudinal and transverse momenta. These numbers obviously depend on target material, density and length and on the number of protons on target.

Secondary mesons are generated along the target according to the distribution of proton interaction points. The latter depends on the proton beam size and divergence, on the target thickness, $`T_{targ}`$, and on the proton interaction length, $`\lambda _p`$, of the target material. Momentum and angular distributions of the mesons are sampled from the proposed parameterization. The weight associated with this step is

$$W_1=(1e^{T_{targ}/\lambda _p})e^{z/\lambda _p}$$

because protons are forced to interact in the target.

Meson trajectories in the target are calculated and their length, $`z^{}(z,\theta )`$ used to estimate the probability that the mesons exit the target without re-interacting: $`e^{z^{}/\lambda _s}`$, where $`\lambda _s`$ is the meson interaction length in the target.

Tertiary contribution is added following the parameterization described in Section 4.3:

$$W_2=e^{z^{}/\lambda _s}\left[1+A_h(x)_0^z\frac{dy}{\lambda _p}w_z(y)\right]$$

where $`A_h(x)`$, $`w_z(y)`$ and the approximation related to this approach have been discussed above. In case of the segmented target considered for the CNGS beam, this model gives the results of figure 8.

#### 5.1.2 Meson tracking in the neutrino beam-line

The trajectory of each meson in the beam-line is calculated, taking into account the tracking in the magnetic field of horn and reflector, until it hits the walls of the decay tunnel or the collimators. The amount of material crossed by the particle is also recorded.

The meson is forced to decay along its trajectory, $`traj`$, accordingly to its decay length, $`\lambda _{dec}`$. The weight associated to this process is

$$W_3=(1e^{traj/\lambda _{dec}})e^{Z_{int}/\lambda _s}$$

where $`Z_{int}`$ is the amount of material crossed up to the decay point and $`\lambda _s`$ is the interaction length in that material.

Contributions due to reinteractions of secondary particles in the material along the beam line are taken into account with a parameterization similar to that used for the target but in the “short target” approximation. The corresponding weight is:

$$W_4=1+A_h^{ter}(1x)^{b_h^{ter}}\frac{Z_{int}}{2\lambda _p}$$

as derived in formulae (15) and (16) of Section 4.3

#### 5.1.3 Neutrino production from mesons

For each meson a neutrino is produced; its flavour and its momentum distribution in the parent meson rest frame depend on the decay mode and branching ratio, $`B.R.`$. The neutrino direction in the laboratory frame is determined requiring that it crosses the detector volume. The angle, $`\theta _{s\nu }`$, between parent meson and neutrino directions allows to calculate the Lorenz boost of the neutrino from the meson rest frame to the laboratory frame. This in turns allows to obtain the neutrino momentum in the laboratory frame.

The weight associated to this process is proportional to the probability that the neutrino is emitted in the detector direction. This is obtained by simply boosting back the solid angle subtended by the detector in the meson rest frame, where the neutrino is emitted isotropically;

$$W_5=B.R.\times (m_s/(E_sP_s\mathrm{cos}\theta _{s\nu }))^2$$

where $`m_s`$, $`E_s`$ and $`P_s`$ are the mass, energy and momentum of the secondary meson.

#### 5.1.4 Neutrino production from muons

Muons are produced in the decay of secondary mesons taking into account the correct kinematics (branching ratio and polarisation). Muons are also tracked through the neutrino beam line and forced to decay to produce neutrinos in the detector direction. An additional weight, $`W_{3\mu }`$ (equivalent to $`W_3`$ for meson decays), is introduced for the neutrinos from muon decay. A weight $`W_{5\mu }`$ replaces $`W_5`$ accounting for the the muon decay kinematics (including polarisation).

### 5.2 Statistical accuracy

The statistical accuracy of this way of simulating neutrino beams does not depend much on the distance between the detector position and the neutrino source as it is the case for classical unweighted methods.

In the classical unweighted case, only a fraction of the pions (typically 5–10%) decay before interacting (either in the beam-line material or in the decay tunnel walls); in addition neutrinos are spread over a wide solid angle (about 1 mrad in the CNGS case) because of the decay kinematics. It follows that, to enhance statistics, at large distance the neutrino spectra need to be computed on a surface much wider than the actual detector area, relying on the fact that the spectra shapes varies slowly with the radius. In the CNGS case an accuracy of a few percent can be achieved with several millions protons on target if a detector area larger than $`10^4`$ m2 is considered.

In the parameterization case, since all mesons – within the focusing optics acceptance – are exploited to produce neutrinos in the detector, the statistical accuracy is independent from the detector distance and proportional to the inverse of the square root of the number of generated positive pions (for $`\nu _\mu `$ beams), namely about the number of generated proton interactions on target. An accuracy better than a percent is thus obtained with less than $`10^5`$ p.o.t., for any size of the detector surface.

### 5.3 Simulations of past and future neutrino beams

In order to give a more quantitative appreciation of the accuracy that one can obtain in the simulation of neutrino beams using the generator described in the previous Sections, we present the comparison with the published neutrino spectra measured with the CHARM II detector CHdet  exposed at the CERN-WANF beam.

In addition we present the comparison between the CNGS simulation based on our method and that based on the FLUKA stand-alone program for secondary meson production interfaced with GEANT for particle tracking.

#### 5.3.1 Comparison with CHARM II data

The WANF neutrino beam line at CERN, to which the CHARM II detector was exposed, is well described elsewhere WANF . Its main characteristics has been summarized in Table 4. The facility was run, during several years of operation, either selecting positive charged particles ($`\nu _\mu `$ beam) or negative ones ($`\overline{\nu }_\mu `$ beam). Neutrino/anti-neutrino interactions were collected in the CHARM II detector and fully reconstructed CHARM .

In figure Figure 10top, we show the comparison between the measured neutrino fluxes in case of positive mesons focusing, ($`\nu _\mu `$ beam with $`\overline{\nu }_\mu `$ contamination) and the simulation performed with our method and based on $`10^5`$ p.o.t.. On the left hand side a logarithmic scale is used to make evident the spectral behaviour at high energy as well as the wrong sign contamination; on the right hand side a linear scale is used for a better appreciation of the focusing/defocusing effect. In figure Figure 10bottom, the same plots are shown in the case of negative meson focusing ($`\overline{\nu }_\mu `$ beam with $`\nu _\mu `$ contaminations).

The overall agreement is at the percent level, with at most 10% disagreement on a bin per bin basis.

The high energy tails of the distributions are dominated by the production of high energy secondary mesons peaked in the forward direction, and are practically insensitive to the magnetic focusing. In fact high energy mesons, with small angular aperture, travel most likely through the neck of the horn where they are hardly deviated. The good agreement between simulation and data indicates that high $`x_F`$ production on target is well simulated and that re-interactions on the material along the beam-line is correctly taken into account.

In the focusing/defocusing energy range the agreement is an indication that low $`x_F`$ production is correctly generated at least up to $`10`$ mrad (the WANF optics acceptance). The fact that also the wrong sign contamination in the simulation behaves as the data, means that tertiary production in target and down-stream material (mainly the horn neck) is described to a sufficient level of approximation.

#### 5.3.2 The CNGS beam-line

As mentioned earlier, our neutrino beam generator was originally developed to allow a rapid optimization of LBL neutrino beams. To test the reliability of it, we have performed a detailed simulation of the CNGS LBL neutrino beam for an extensive comparison with the full beam simulation based on the FLUKA stand-alone package FONE  for secondary particle production and on GEANT for tracking along the beam line CNGS .

Apart from the tunnel geometry and the focusing optics, the main differences of the CNGS beam-line with respect to the WANF are the target material ($`carbon`$ instead of $`beryllium`$) and the proton beam energy ($`400`$ GeV instead of $`450`$ GeV). These differences have been taken into account in our simulation with the scaling laws proposed and described in Section 4.2 and 4.3.

As remarked in Section 4.3, the agreement at the target level between our calculation and the FLUKA based full Monte Carlo program is noticeable, provided that tertiary production is taken into account in the former case. The slight disagreement visible in the $`p_T`$ distributions of all charged particles mainly at large $`p_T`$ translates in a slightly more forward peaked distribution of the exit points along the target.

Remarkably enough, as far as the neutrino spectra at large distance are concerned, the discrepancies between the two models do not propagate with the same strength (see Figure 11). This is because in the momentum range and angular acceptance selected by the focusing optics both particle production models reproduce very well the available experimental data.

Before concluding, it is worth mentioning that early simulations of the CNGS beam line based on GEANT stand-alone disagreed with those presented here, being too optimistic by more than $`20\%`$. On the other hand FLUKA stand-alone gives results fully compatible with those presented here.

## 6 Conclusions

Empirical formulae for single-particle inclusive invariant cross sections in $`p`$-Be interactions have been derived, on the basis of single-particle inclusive production data collected by the NA20 NA20  and NA56/SPY SPY1 ; SPY2  experiments. These formulae reproduce the experimental data within a 10% accuracy.

The hypothesis of Feynman scaling has been verified to hold with our parametrisation giving a suitable representation of production data collected over a wide range of primary proton beam energies (from 24 GeV to 450 GeV).

Prescriptions to extrapolate this parameterization to finite targets and to targets of different materials have been given.

The results obtained have been used as an input for the simulation of neutrino beams. A comparison to data collected by the CHARM-II neutrino experiment at CERN has demonstrated the capability of this approach to predict the past.

These formulae can be of great practical importance for fast calculations of neutrino fluxes and for designing new neutrino beam-lines. Thus they may be used in fast simulations aiming at the optimisation of the long-baseline neutrino beams at CERN and FNAL. Predictions for the neutrino spectra of the CNGS beam from CERN to Gran Sasso have also been given.

## Acknowledgements

We gratefully acknowledge enlightening discussions with and fruitful suggestions of K. Elsener, L. Gatignon, P.G. Ratcliffe, S.I. Striganov. We are indebted to P. Gorbunov, for having provided us the files with neutrino spectra measured in CHARM II. We are grateful to S. Ragazzi, our spokesman in the NA56/SPY collaboration, and P. Picchi for their encouragement to this work and to all our colleagues in the NA56/SPY collaboration.