Prompt neutrinos and intrinsic charm at SHiP

We present a new evaluation of the far-forward neutrino plus antineutrino flux and number of events from charm hadron decays in a 400 GeV proton beam dump experiment like the Search for Hidden Particles (SHiP). Using next-to-leading order perturbative QCD and a model for intrinsic charm, we include intrinsic transverse momentum effects and other kinematic angular corrections. We compare this flux to a far-forward flux evaluated with next-to-leading order perturbative QCD, without intrinsic transverse momentum, that used the angular distribution of charm quarks rather than the neutrinos from their decays. The tau neutrino plus antineutrino number of events in the perturbative QCD evaluation is reduced by a factor of about three when intrinsic transverse momentum and the full decay kinematics are included. We show that intrinsic charm contributions can significantly enhance the number of events from neutrinos from charm hadron decays. Measurements of the number of events from tau neutrino plus antineutrino interactions and of the muon charge asymmetry as a function of energy can be used to constrain intrinsic charm models.


Introduction
The potential for the discovery of dark matter and hidden particles in theories beyond the standard model are the motivation for new beam dump experiments under review or design. The SHiP (Search for Hidden Particles) experiment is one such project, now in the comprehensive design phase [1,2]. The experiment plans to run with 2 × 10 20 protons on target over a 5-year run using the SPS 400 GeV proton beam at CERN. While a major goal is to detect or constrain physics beyond the standard model, the protons incident on the beam dump target will produce a significant number of neutrinos. Designed to minimize the production of pions and kaons to optimize the beam dump as a source of hidden particles, guaranteed outputs are neutrinos from charmed mesons that decay promptly, as well as some neutrinos from pion and kaon decays. The short lifetimes of charmed mesons mean that they decay before losing energy by hadronic scattering in the target.
Of particular interest are tau neutrinos, as the direct detection of tau neutrinos is thus far limited to 9 ν τ events at DONUT [3,4] and 10 tau neutrino candidate events at Opera [5,6]. In ref. [2], an evaluation of the tau neutrino flux with some simplifying approximations was made, predicting on the order of 1,000 tau neutrino and antineutrino interactions in the SHiP detector. Muon neutrino and antineutrinos are also produced from charm decays, at a level that is a factor of ∼ 10 larger that tau neutrinos and antineutrinos. The charged and neutral D mesons decay semileptonically to muon neutrinos and electron neutrinos. Cabibbo and kinematic suppression factors mean that the D ± mesons have only a tiny branching fraction to ν τ in two-body decays, and the D 0 ,D 0 , not at all. The D ± S is heavy enough, and not suppressed by the Cabibbo angle in its decay, so that D − S → τν τ and its charge conjugate have a branching fraction of (5.48 ± 0.23)% [7]. Subsequent tau decays, also prompt, produce a second tau neutrino for each tau antineutrino that comes directly from the D − S . The SHiP experimental configuration, with a beam dump and 51.5m distance between the front of the target and the detector, collects particles in the very forward direction, for a narrow range of angles θ relative to the beam direction. For the nominal SHiP detector, angles θ < 7.3 × 10 −3 − 1.9 × 10 −2 rad reflect the width and height of the detector 2m × 0.75m, respectively. In terms of pseudorapidity, this means η > 4.6 − 5.6. In ref. [2], the neutrino flux was evaluated in a collinear approximation using next-to-leading order (NLO) QCD for the pA → ccX cross section. In the collinear approximation used there, the charm quark momentum direction was used to determine η c , which was used as approximate cut for whether or not the neutrinos and antineutrinos were directed into the detector. Following the charm momentum angular (or pseudorapidity) cuts, the energies of the mesons, taus and eventually neutrinos were traced using analytic energy distributions where relative angles are already integrated.
In this paper, as a proxy for detector angular cuts, we restrict pseudorapidity values to η > 5. 3. When the cut is applied to the charm pseudorapidity (η c > 5.3), with the same parameter inputs as in ref. [2], we find a tau neutrino plus antineutrino event rate very close to the results in ref. [2] which is calculated with a more detailed detector geometry. The simplified restriction with pseudorapidity is 10% lower than the event rate using the detector geometry. A value of η = 5.3 corresponds to a distance of 0.5m from the beam axis at the detector. Here, we refine the evaluation of the flux of tau neutrinos and antineutrinos at SHiP to include the angular dependence of the charmed meson decay. We also incorporate intrinsic transverse momentum effects in the charmed hadron production. The primary impacts of the calculational improvements discussed here to allow for a cut on the neutrino pseudorapidity, equal to the neutrino rapidity, η ν = y ν with η ν > 5.3, are to reduce the flux of neutrinos into the detector. We discuss the impact of each correction below.
The prompt production of neutrinos from charm decays is also a topic of considerable interest in the context of the prompt atmospheric neutrino flux, for example, in refs. [8][9][10][11][12][13][14][15][16][17]. Because the spectrum of cosmic rays incident on air nuclei falls steeply with cosmic ray energy, φ(E CR ) ∼ E −3 CR , a high charm energy fraction is favored for the prompt neutrino flux, even if this kinematic regime is a relatively small portion of the phase space for the cross section at a fixed beam energy. Gluon fusion, gg → ccX is the dominant parton model process in perturbative production of charm [18][19][20][21][22]. Other approaches (e.g., [23][24][25]) including intrinsic charm [26][27][28][29][30][31][32][33][34] have been proposed as a mechanism to enhance the high energy atmospheric lepton flux [14][15][16][17]. In this approach, heavy and light quarks are coalesced to produce high momentum charmed hadrons. We explore here how signals of intrinsic charm may appear in the spectrum of the tau neutrino plus antineutrino flux at SHiP. Efforts to measure production of charmed mesons with the 400 GeV beam at CERN, e.g., the DsTau project [35], would constrain forward production of charmed hadrons directly.
Neutrino-antineutrino flux asymmetries from intrinsic charm are potentially observable in the muon neutrino and antineutrino rates. The muon neutrino and antineutrino fluxes from charm are evaluated in the NLO QCD approach and for intrinsic charm. The asymmetry in the NLO perturbative QCD evaluation is quite small, since gluon fusion dominates the cross section. On the other hand, for five or seven particle states (e.g., |uudcc , |uudccuū ) where valence quark momenta is brought to, for example, D 0 = uc, but not a D 0 =ūc, particle-antiparticle asymmetries can be large. We make an approximate evaluation of the muon neutrino-antineutrino asymmetries relying on our perturbative charm evaluation, intrinsic charm and the approximate flux of neutrinos at SHiP from pion and kaon decays based on fluxes from ref. [1]. We begin with the production of charmed quarks and hadrons in the parton model and with a model for intrinsic charm based on ref. [26][27][28] discussed in detail in ref. [29]. Our assumptions about the intrinsic transverse momentum are included in this section. In section 3, we outline the procedure used for modeling the charmed hadron decays. Section 4 shows our results for neutrinos from prompt decays of charm for proton beam energies of 400 GeV. We discuss how the energy dependence of the tau neutrino flux and asymmetry measurements of muon neutrino and antineutrino event rates could show evidence of intrinsic charm. Our conclusions appear in section 5. Details of the inputs to the intrinsic charm evaluation, and our approximations for the fluxes of muon neutrinos and antineutrinos from pions and kaons, are given in the appendix.

Production of charmed quarks and hadrons
Charmed quark production in the parton model in perturbation theory is evaluated with partons collinear with the incident hadrons. The factorization scale µ F characterizes this collinear approximation, with the interpretation that, roughly, parton transverse momenta below µ F are approximately collinear. At SHiP, because of the small angular size of the detector downstream of the target, even small transverse momentum effects are potentially important. We discuss below first, the standard next-to-leading order evaluation of perturbative charm production. Then, we show our implementation of "intrinsic" transverse momentum to model the low transverse momentum of the partons and parameters of quark fragmentation to hadrons. Finally, we show how intrinsic charm could contribute to charmed hadron production at SHiP.

Perturbative charm
The parton model evaluation of charm production involves the production of a charmed quark pair, followed by the hadronization of the quark and antiquark into a charmed meson or baryon. Charm production can be written in terms of the parton level hard scattering cross sectionσ ij that depends on the strong coupling constant α s (Q 2 ) evaluated at a renormalization scale Q 2 = µ 2 R . The differential cross section at the parton level is convoluted with the parton distributions for partons i and j via the parton distribution functions f H1 i (x 1 , µ 2 F ) and f H2 j (x 2 , µ 2 F ), with x 1,2 equal to the momentum fractions of the participating partons for hadronic collisions H 1 + H 2 → ccX. At leading order in QCD, gluon fusion dominates the production of cc pairs, as it does at next-to-leading order (NLO) [18][19][20][21][22].
The differential cross section for the inclusive c cross section for charm momentum p is for charm production, so σ(pA) Aσ(pN ) in this approximation.
Theoretical evaluations of perturbative production of charm have a long history, with evaluations through NLO in QCD in Refs. [18][19][20] that have been implemented in the HVQ computer program, as described in refs. [21,22]. Our perturbative evaluation uses the HVQ computer program as a basis for the NLO evaluation. For the results presented here, we take the charm quark mass to be m c = 1.27 GeV, and the central values of the factorization scale µ F = 2.1m c and renormalization scale µ R = 1.6m c , following the analysis presented in ref. [36]. The scale choices will permit us to make a direct comparison with the results in ref. [2]. We use the CT14 parton distribution functions (PDFs) [37]. The charm pair production cross section at next to leading order with these parameter choices for a proton beam energy of E p = 400 GeV incident on an isoscalar nucleon target is 19.7 µb.
Since perturbative charm is produced primarily through gluon fusion, the charmed hadron and anti-hadron production rates are equal. This will distinguish perturbative charm production from other production mechanisms like intrinsic charm, discussed below.

Intrinsic Transverse Momentum
For the large rapidities required by SHiP, intrinsic transverse momentum can deflect charm quarks from the detector direction, so it must be included. Schweitzer, Teckentrup and Metz [38] showed that a Gaussian model of intrinsic transverse momentum is compatible with semi-inclusive deep-inelastic scattering and with Drell-Yan. In the Gaussian approach, azimuthal symmetry around the beam axis permits us to write normalized to unity with the integral over d 2 k T . For the Gaussian distribution, The differential distribution for charmed quark production becomes In ref. [38], it was shown that for leading order Drell-Yan production, the Gaussian approximation provides a good model if k 2 T,DY = 1.4 GeV 2 for a 300 GeV proton beam scattering on platinum (FNAL-288). For pion scattering from tungsten (FNAL E615), the value is a little larger (1.7 GeV 2 ). To illustrate the impact of intrinsic transverse momentum, we take a standard value of k 2 T = 1.5 GeV 2 for the average intrinsic momentum squared. This is equivalent to k T 1.1 GeV. We compare with using a larger value, k T = 1.5 GeV instead.

Fragmentation
The Peterson fragmentation function is used to describe charm quark fragmentation into D mesons.
For hadrons H equal to charged and neutral D's and for Λ c , we can write  Table 1. Parameters for the charm quark fragmentation to H = D mesons and the H = Λc, from ref. [40] table I, and the fragmentation fractions of ref. [41] Kramer [40], with µ 0 = 1.5 GeV, and we rescale N H to normalize the fragmentation functions to the fragmentation fractions of ref. [41]. The parameters H and the fragmentation fractions are shown in table 1. We do not evolve the fragmentation functions.

Intrinsic charm
Beam dump experiments have the potential to constrain additional contributions to charmed hadron production, beyond perturbation theory, through neutrino measurements. The description of high momentum charmed hadrons have long been a topic of interest, e.g., as discussed in refs. [23][24][25][26][27][28][29][30][31][32][33][34]. Often cited are the data from the LEBC-MPS collaboration from 800 GeV protons incident on a bubble chamber [42], and the SELEX (E781) collaboration's measurements of high momentum Λ c 's in π − , Σ − and p beams incident on copper and carbon targets [43]. These data show an excess of events for high momentum fractions, beyond the perturbative prediction in this energy regime, however, there are large error bars. They also found asymmetries in the production of hadrons with components that could come from beam valence quarks compared with their antiparticles.
For definiteness to find charm contributions to neutrino fluxes at SHiP, we use an intrinsic charm model of Brodsky, Hoyer, Peterson and Sakai (BHPS) [26,27], described in detail in refs. [28,29,44]. There, intrinsic charm (IC) is modeled by considering protons with 5-quark and 7-quark Fock states in which two of the quarks are a cc pair. Heavy and light quarks can coalesce into high momentum charmed hadrons. While this is not the dominant charmed hadron production mechanism, it does dominate the large Feynman x F distribution.
Based on frame-independent probability distributions for intrinsic charm P ic , two types of contributions to the production of charmed hadrons are included: independent (uncorrelated) fragmentation (F), and coalescence distributions (C) specific to individual charmed hadrons. They depend on the n-particle probability functions [29], for n = 5 and n = 7. The uncorrelated fragmentation is assumed to have a delta function fragmentation function for c → H, so as in ref. [29], For coalescence, the probabillty function includes a delta function that combines the momentum fractions of the quarks that make up the final state hadron, In each of the above equations, x H is the Feynman x F variable for the given hadron H. In the large In what follows, we approximate x F x E in all of our evaluations of the neutrinos from intrinsic charm.
The relative contributions of the five quark (uudcc) and seven quark (uudccqq) terms for qq = uū, dd, ss are set according to ref. [29]. The normalizations are P 5 ic = 3.1 × 10 −3 , and P 7 icu = P icd = (m c /m u ) 2 P ic 7 , P 7 ics = (m c /m s ) 2 P 7 ic where P 7 ic = 4.4 × 10 −2 P 5 ic andm u =m d = 0.45 GeV,m s = 0.71 GeV andm c = 1.8 GeV. For completeness, we write the differential distributions for the probabilities for each of the charmed hadrons as a function of x F in the appendix.
The overall normalization of intrinsic charm remains to be fixed, the value of which is of great interest in the context of the prompt atmospheric flux [14][15][16]. Charm production in the atmosphere, followed by its prompt decay into leptonic or semileptonic decay modes, is the dominant contribution to the atmospheric neutrino flux at sufficiently high energies. The IceCube measurements of the diffuse astrophysical flux have ruled out an extremely large intrinsic charm component, as discussed in, e.g., ref. [15]. Nevertheless, intrinsic charm may play an important role in the high energy atmospheric neutrino flux because the incident cosmic ray hadrons have a very steep energy dependence. Laha and Brodsky in ref. [16] argue that IceCube will be able to constrain intrinsic charm in the future via their neutrino measurements in the 100 TeV to PeV energy range. Laha and Brodsky normalization of the intrinsic charm cross section in pN collisions for √ s = 39 GeV [16]: The intrinsic charm cross section in this approach is expected to scale with energy as the pA total cross section does, so only weakly with energy. We use this same normalization for our evaluation of the intrinsic charm contribution to the neutrino flux at SHiP (labeled IC25). The weak energy dependence of intrinsic charm means that the relative contributions of intrinsic to perturbative charm production is larger at beam energies of 400 GeV than at the LHC. To account for nuclear effects, we use the appropriate scaling of A α where α = 0.71 [45]. For molybdenum, A 0.71 /A = 0.27. We vary the overall IC normalization in our discussion of the potential signals of intrinsic charm in the neutrino signal, for example, considering a normalization of the pN differential cross section at x F = 0.32 of 7 µb (IC7), also considered in ref. [16]. Asymmetries between D − and D + , and between D 0 andD 0 , for the x F distributions are predicted in intrinsic charm models. Asymmetries are smaller for D ± S . The detailed probability distributions listed in the appendix quantify the differences in particle and antiparticle production as a function of x F .
In this paper, we will also explore the implications for asymmetries in the high energy ν µ and ν µ energy distributions. Neutrino and antineutrino production from pion and kaon decays are also asymmetric, but the beam dump will suppress pion and kaon production of muon neutrinos and antineutrinos at high energies. We make an approximate evaluation of the muon neutrino and an- tineutrino flux from pion and kaon decay for the SHiP target configuration based on fluxes presented in ref. [1]. We compare them with contributions from intrinsic and perturbative charm.

Charmed hadron production
In this section, we begin with our results for charmed hadron production via the perturbative approach for E p = 400 GeV. Fig. 1 shows the differential cross section as a function of the charmed hadron energy for each of the charmed hadrons considered here. To compare the distributions, the distributions have been normalized by 1/B c→H . The energy distribution, integrated over the full phase space, does not change significantly with the introduction of intrinsic transverse momentum, even for k T = 1.5 GeV. The energy distribution of the charmed hadron does change when an angular cut is applied. To illustrate, we show in fig. 2 the energy distribution of the D S when η D S > 5.3. The intrinsic k T along with the pseudorapidity cut eliminates approximately 50% of the low energy D S 's when no intrinsic k T is included with the cut. This translates to a lower flux of low energy tau neutrinos at the SHiP detector. This is not unexpected. Since k T 1 GeV and θ ∼ k T /E, θ 0.01 rad means the energies out to ∼ 100 GeV are impacted.
The charmed hadron energy distributions for both perturbative production and via intrinsic charm are shown in fig. 3 as a function of energy (upper) and lab rapidity (lower). The upper figure shows that intrinisic charm can dominate the energy distribution of charmed hadrons for hadron energies larger than E H ∼ 50 − 100 GeV, depending on the normalization of the intrinsic charm contribution. As noted above, we normalize the intrinsic charm according to eq. (2.11). Intrinsic charm clearly dominates at high rapidity, as anticipated.
The perturbative result for D + equals the result for D − , and similarly for the other hadrons, however, the intrinsic charm contributions are different for hadrons incorporating a charmed quark compared to charmed antiquark. The Λ + c are produced at the highest rapidities, as can be seen seen in fig. 4, where perturbative and intrinsic contributions to charmed hadron rapidities are shown on a linear scale.
Perturbative charm and intrinsic charm production D S 's may be measured in the near future. The DsTau experiment proposes to study D S production from a double kink signature, where the decay D S → τ ν τ and the τ → ν τ X will be observed in an emulsion detector. They aim for 2.3 × 10 8 proton interactions in a thin tungsten target to get an estimated 1000 double kink events when a 20% detector efficiency is included. Our evaluation of the differential distribution for D ± S for the tungsten target is shown in fig. 5. The lower dashed histogram comes from the perturbative QCD evaluation with intrinsic transverse momentum k 2 T = 1.5 GeV 2 , while the green dot-dashed histogram shows the intrinsic charm contribution to D ± S normalized by eq. (2.11). The sum of the two is shown with the solid blue histogram. The dotted histogram shows the sum with the intrinsic charm reduced by a factor of 7/25. For perturbative QCD production, with a branching fraction of B(D S → τ ν τ ) = 5.5% and the factor of 20% for detector efficiency, the number of double kink events is 917 with our evaluation, integrated over all energy, with 326 events coming from D ± S with E > 24 GeV. For the intrinsic charm contribution with the full energy range, the number of events from IC is 272, with 256 events having D ± S energies above 24 GeV. For E > 100 GeV, there are 21 perturbative events and 93 IC events of this type. The lower normalization factor of 7/25 for IC reduces the number of IC events by a factor of 0.28. Approximating has a power of n 10.6 at large x E for the perturbative result, and n 4.6 for the sum of perturbative and intrinsic charm as evaluated by this model. The DsTau experiment should be able to put strong constraints on the intrinsic charm model of ref. [29].
The direct and chain decay neutrino energy distributions from π ± and µ ± decays are discussed in detail in ref. [46][47][48]. The evaluation is directly transferable to D S → τ ν τ . We sketch the result here. In the D s rest frame, τ s are produced isotropically since the D S has zero spin. For this two-body decay, in the rest frame of the D S , the tau and neutrino energies are where r = (m τ /m D S ) 2 = 0.815, m Ds = 1.97 GeV and m τ = 1.78 GeV. We use primed variables to denote quantities in the D S rest frame. One can find the corresponding distributions in the lab frame through a combination of boosts and rotations. Ref. [49] has a discussion of the details of the boosts, for example. In the tau rest frame, the distribution of the massless (or nearly massless) leptons from the threebody decay of the tau is given by [46][47][48]  where x * ≡ 2E * l /m τ ∈ [0, 1], and θ * is the angle between the lepton and the spin of the tau. For ν τ or ν τ , f 0 (x * ) = 2x * 2 (3 − 2x * ) and f 1 (x * ) = 2x * 2 (1 − 2x * ). The starred variables denote the tau rest frame. Similar results apply to other tau decay modes.
The tau decay distribution in its rest frame, for angles of lepton l = ν τ (ν τ ) or D S momenta relative to an axis in the lab direction of the tau momentum, can be rewritten as where the opposite helicities of ν τ andν τ conspire to yield a single equation for both neutrino and antineutrino distributions, with [46,47] cos θ * where E D S and E τ are the energies of D S and τ in the lab frame. The results shown below use eq. (3.3) boosted to the lab frame for the tau leptonic decay. The remainder of the tau decays are approximated by two-body decays, to πν τ , ρν τ and a 1 ν τ . The decay distributions can be found in ref. [11]. As noted above, particle-antiparticle asymmetries in charm hadron production will not be evident for D ± S but will be for other D mesons and the Λ c . To estimate signals in SHiP, we consider the decays H → ν µ X. For the muon neutrino and antineutrino distributions shown here, we have used an approximate energy distribution equal to the muon energy distribution. As implemented in the HVQ program, the parameterization of the muon distribution follows a naive spectator model which depends on the mass ratio r = m 2 K /m 2 H , where m H is the charmed hadron mass.
4 Results for neutrinos from prompt decays of charm

Neutrino and antineutrino energy distributions
The NLO QCD result for the neutrino differential cross section is shown in fig. 6 for E p = 400 GeV. The distributions are the sum of neutrinos and antineutrinos. The neutrino or antineutrino rapidity is required to be larger than η ν > 5.3. The upper curve shows the charm contribution to the energy distribution of ν µ andν µ . The lower solid curve shows the tau neutrino plus antineutrino distribution, the sum of the dashed (direct decays) and dotted (chain decays from tau) curves. As noted above, the direct neutrino has a much smaller energy than the tau in the D S rest frame, so even after the tau decay, the direct neutrinos are less energetic than those from the D S → τ → ν τ chain decay. The crossover from direct decay dominated to chain decay dominated is at E ν ∼ 20 GeV. The prompt muon neutrino plus antineutrino differential cross section is approximately a factor of ten larger than that of prompt tau neutrinos plus antineutrinos. Since charm is produced with anticharm, the differential cross sections to produce ν andν are the same for a given . We approximate the electron neutrino and muon neutrino differential cross sections to be equal since the charged lepton masses m e and m µ are small compared to the charmed meson masses.
It is not surprising that with the normalization of intrinsic charm discussed in the previous section, intrinsic charm contributions to the neutrino differential cross section in the very forward region overwhelm the perturbative neutrino differential cross section, as shown in fig. 7 for the tau neutrino and antineutrino distributions for protons incident on a molybdenum target at SHiP. The lower dashed histogram shows the perturbative tau neutrino energy distribution (equal to the perturbative tau antineutrino energy distribution), while the solid lines show the sum of intrinsic and perturbative tau neutrino, and separately, tau antineutrino energy distributions. The IC is normalized with eq. (2.11), where with A = 96, the normalization factor is A 0.71 /A = 0.27. The default intrinsic transverse momentum k 2 T = 1.5 GeV 2 is used, with η ν > 5.3. Similarly, the dashed line shown in fig. 8 shows the ν µ energy distribution dσ/dE per nucleon in the target, equal to theν µ energy distribution evaluated using NLO perturbative QCD of cc with k 2 T = 1.5 GeV 2 for η ν > 5.3. The solid lines show the sum of perturbative charm and intrinsic charm contributions to the muon neutrino (lower solid histogram), and muon antineutrino (upper solid histogram), energy distributions for the same k 2 T and rapidity restriction. Since charmed hadrons D − , D 0 and Λ c are favored over their charged conjugates with this model for intrinsic charm momentum distributions, the production rate for anti-neutrinos is larger than for neutrinos. Above E ∼ 20 GeV, differences in the asymmetry between particle and antiparticle production with and without intrinsic charm will be apparent if the intrinsic charm contribution is normalized by eq. (2.11). We discuss the muon neutrino-antineutrino asymmetry below, where we show that the enhanced flux of muon antineutrinos compensates in part for the lower antineutrino cross section in the detector.
Intrinsic charm dominates muon neutrino and antineutrino energy distributions for most of the energy range. The intrinsic charm cross section normalization can be constrained at SHiP from high energy muon neutrino and antineutrino event rate measurements. Indeed, SHiP is better positioned to constrain the normalization of the intrinsic charm contribution than current LHC experiments. At LHC center-of-mass energies, the intrinsic charm cross section is much smaller fraction of the perturbative charm cross section.

Neutrino and antineutrino cross sections
To convert the differential cross section to event rates with a charged lepton, the neutrino charged current cross section is required. We use PDFs for a lead neutrino target, modeled using the deep inelastic scattering (DIS) cross section with the nCTEQ parton distribution functions for lead [37] according to the charged lepton mass m corrected differential cross section for ν (k)+N (p) → (k )+X process with Q 2 = −(k − k ) 2 = −q 2 [50][51][52][53]. We reproduce the expression here, in terms of the structure functions F i = F i (x, Q 2 ), x = Q 2 /2p · q, y = p · q/p · k and the neutrino energy in the fixed nucleon target frame E ν , We use the low momentum transfer extrapolation of the DIS structure functions outlined in refs. [53,54]. Shown in fig. 9 are the neutrino and antineutrino cross sections per nucleon per neutrino energy as a function of incident neutrino energy for the production of muons and taus. We require that the hadronic final state energy W be larger than W min = m p + m π for DIS. At low energy for ν µ interactions, the impact of W > W min can be seen in fig. 9 for muon neutrinos and antineutrinos. For tau neutrinos, the suppression of σ/E comes from both W > W min and tau mass kinematic corrections. The tau mass corrections persist out to relatively high energies. The muon neutrino charged current cross section is about a factor of 1.5 larger than the tau neutrino cross section at E ν = 50 GeV. The ratio reduces to ∼ 1.10 for E ν = 400 GeV. Tables 2 and 3   . The deep-inelastic neutrino and antineutrino lead scattering cross section, per nucleon, for muon and tau (anti-)neutrinos as a function of neutrino energy. We require that the hadronic final state energy be larger than Wmin = mp + mπ.

Event rates
The event rates depend on the neutrino flux and the neutrino cross section. For the flux of neutrinos from charm decays, we write As in ref. [2], we take N p = 2 × 10 20 protons, σ(pA → νX) = Aσ(pN → νX) and σ(pA)/A = σ pN = 10.7 mb. The event rate is Event rate = φ ν L λ ν (4.5) where λ ν = (N A σ(νA)/A) −1 and L = 66.7 g/cm 2 for lead in the detector. Fig. 10 shows the impact of intrinsic transverse momentum and the full three-dimensional kinematics for the sum of tau neutrino plus antineutrino events for η ν > 5.3, as a function of neutrino energy, evaluated using NLO perturbative QCD. The upper blue histogram shows the number of events for tau neutrinos plus antineutrinos, evaluated using the same approximations as in ref. [2]: that the differential cross section for charmed quark production, evaluated with CT10 PDFs with no intrinsic transverse momentum, in the approximation that the charm quark momentum is in the same direction relative to the incident proton beam as the neutrino from its decay (the "collinear approximation"). Here we show for η ν > 5.3 as an approximate representation of the detector solid angle. The total number of events shown in fig. 10 is within 10% of the number of events evaluated in ref. [2] where a more detailed detector geometry is used. Using CT14 PDFs reduces the event rate by ∼ 10% compared to the rate with CT10 PDFs. With intrinsic transverse momentum k 2 T = 1.5 GeV 2 , the event rate is 65% of the rate without intrinsic transverse momentum. The three dimensional  decay kinematics further reduces the overall event rate. The number of events goes from 972 to 307 events for N p = 2 × 10 20 with these parameter choices for NLO perturbative production of charm with fragmentation to D S mesons that decay to τ ν τ .
The choice of the scale factor has an impact on the number of events. As in ref. [2], we use a range of scales guided by ref. [36], namely (µ R , µ F ) = (1.48, 1.25)m c and (µ R , µ F ) = (1.71, 4.65)m c to bracket the predicted charm production cross sections. Fig. 11 shows the event rates with the bands bracketed by the scale choices. The range in the number of events is 222 events to 404 events (sum of tau neutrino plus antineutrino). The fluxes vary by less than ∼ ±40% with these scale choices below E = 100 GeV, less than ∼ ±30% below 50 GeV. Intrinsic charm, at the level with normalization in eq. (2.11), makes a large contribution to the overall event rate of tau neutrinos (middle, orange histogram) and tau antineutrinos (lower, green histogram) shown in fig. 12. The upper (blue) histogram is the sum of tau neutrino plus antineutrino events. The total number is 1688 events with these parameter choices. If the normalization of the intrinsic charm contribution in pN collisions is lowered such that dσ/dx F = 7 µb, the event rate from intrinsic charm reduces by a factor of 0.28. The double peak structure comes from the separate direct tau neutrinos from the D S decay (lower energy peak) and the chain decay from D S → τ → ν τ with the higher energy peak. Including smaller η ν values (including larger angles around the beam direction) merges the two peaks. The asymmetry between the numbers of tau neutrino and antineutrino events comes from the differences in the cross sections. The tau neutrino and antineutrino fluxes are nearly identical for intrinsic charm, although this is not the case for muon neutrinos and antineutrinos.
We now turn to the number of muon neutrino and antineutrino events. We begin with NLO perturbative QCD predictions, including intrinsic transverse momentum and three-dimensional kinematics, for muon neutrinos and antineutrinos from charm production and decay. Fig. 13 shows the number of muon neutrino events from charm as a function of energy and antineutrino events for the 400 GeV proton beam with the standard parameter choices. The two dashed histograms show the intrinsic charm number of events. While for the perturbative evaluation, the difference in neutrino and antineutrino number of events is due to the different neutrino and antineutrino cross sections, the number of events for intrinsic charm reflects enhanced forward production of hadrons with a charmed quark or antiquark that combines with valence quarks. Charm leptonic decays go to neutrinos, and anticharm decays to an antineutrino. The meson components dominate over Λ ± c , so the flux of antineutrinos is higher than neutrinos. The asymmetry in particle and antiparticle production with intrinsic charm is an interesting feature to use as a diagnostic for intrinsic charm.
While muon neutrinos and antineutrinos come from charm decays, they also come from pion and kaon decays, even in the beam dump environment. A full simulation of the muon neutrino plus antineutrino flux at SHiP is beyond the scope of this work, but we can use the approximate flux evaluated in ref. [1] for SHiP. The details of our normalization of their results (presented in ref. [1] with arbitrary normalization) for the muon, electron and tau neutrino plus antineutrino fluxes are discussed in the appendix. We review them briefly here. Our strategy is to use the tau neutrino plus antineutrino flux to set the overall normalization. Assuming simple power laws and taking the electron neutrino plus antineutrino flux from charm, we can extract the contribution from kaons to ν e +ν e . The kaon contribution to ν µ +ν µ is approximated, and the pion contribution to the ν µ +ν µ flux is extracted. Charged particle production ratios [55] are used to approximate the separate neutrino and The flux of ν µ +ν µ from pions falls quickly with energy relative to those from kaons because of the relatively long pion lifetime. While the figure shows an energy range to very small E ν , the simple power law certainly fails at low energy. Below 40 GeV neutrino energy, our approximate flux is less reliable. The number of muon neutrino and antineutrino events from pQCD production of charm and our modeled pion and kaon neutrinos is shown in fig. 15 with the green: the upper solid line for ν µ charged current events and lower dashed line green line forν µ events (no intrinsic charm). Low energy contributions come from pions, however, for ν µ , the kaon contribution dominates for E larger than ∼ 40 GeV. For antineutrinos, the crossover between pion and kaon contributions is also at ∼ 40 GeV, at which point the prompt contribution from charm is nearly equal to the kaon contribution to theν µ event rate. A shaded band with a ±40% error estimate on the antineutrino flux is shown. The shaded band for the same error on the neutrino flux does not appear here because of the log plot.
The nearly flat orange curves in fig. 15 show the ν µ (solid) andν µ (dashed) events where they come from only intrinsic charm production. The intrinsic charm event rate for ν µ 's is approximately equal to that for muon neutrinos from pions, kaons and perturbative charm, and forν µ from intrinsic charm. One can define the charge asymmetry as where dN/dE = φ · L/λ and λ = (N A σ CC ) −1 . Fig. 16 shows the charge asymmetry defined in eq. (4.6). The highest histogram shows the asymmetry evaluated with only pion, kaon and perturbative charm contributions to the ν µ +ν µ flux. The lowest curve shows the case where the intrinsic charm is normalized according to eq. (2.11) for protons incident on molybdenum added to the perturbative charm, pion and kaon contributions. Intermediate curves show the results when the intrinsic charm is decrease by a factor of 7/25 following the discussion of Laha and Brodsky [16], and decreased by a factor of 1/10. The error bands arise from the estimate of ±40% error on the perturbative ν µ +ν µ flux. Therefore, it is our estimate that the intrinsic charm event rate is large enough at high energy that even with 10% of the nominal normalization, the muon charge asymmetry will reflect its presence at that level.

Discussion
Tau neutrino and antineutrino detection at SHiP would show not just evidence of tau neutrino charged current production of taus, but also allow for the study of tau neutrino interactions and the details of tau neutrino production. By approximating the detector acceptance for the SHiP detector a distance of 51.5m downstream from the front of the target with an angular cut equivalent to η ν > 5.3, we can reproduce our earlier results in ref. [2] to within 10%. Our new NLO perturbative QCD evaluation of charmed mesons, in particular for D S that includes intrinsic transverse momentum effects in the hard scattering and the angular effects in the decays of the D S → τ ν τ and τ → ν τ X, reduces the tau event rate by a factor of ∼ 3 compared to the event rates in ref. [2], predicting more than 300 events with our central renormalization and factorization scale choices. As in ref. [2], scale uncertainties introduce a ∼ ±30% uncertainty in the event rate.
Forward production of charm that is not perturbative can enhance neutrino production in beam dump experiments. We have used the BHPS model for intrinsic charm [26][27][28][29] as a representative model to show the impact that forward charm production can have on the tau neutrino and antineutrino event rates. With an intrinsic charm normalization suggested in the context of the prompt atmospheric neutrino flux, which also emphasizes forward production, we find that intrinsic charm production and decay to tau neutrinos and antineutrinos would significantly enhance the event rate. Indeed, intrinsic charm is relatively more important at low energies than for high energy interactions because the cross section is assumed to scale like the total hadronic cross section, so weakly dependent on energy, unlike the perturbative charm cross section. With the intrinsic charm normalization suggested in ref. [16], the event rate for tau production from intrinsic charm alone is almost 1,700 tau neutrino plus antineutrino charged current events.
One of the important features of intrinsic charm models is the asymmetry in particle production because of the valence momenta in the proton beam. The valence quarks are not relevant to D ± S production, so particle-antiparticle asymmetries won't appear in the tau neutrino and tau antineutrino event rates, but they will appear in muon neutrino and antineutrino charged current event rates. To estimate the µ − − µ + asymmetry as a function of energy, we used the Monte Carlo evaluation of the neutrino fluxes from ref. [1], then used particle production asymmetries [55] to estimate the separate pion and kaon contributions to the muon neutrino and antineutrino fluxes. Signals of intrinsic charm would be offsets in the particle asymmetry predicted from just pion and kaon contributions to the 0 20 40 60 80 100 120 140 We have used the BHPS model for intrinsic charm to compare with what has been used to evaluate prompt atmospheric neutrino fluxes in ref. [16], however, this is not the only model with forward charm production or for intrinsic charm. Within the BHPS model, there are choices such as a simplified fragmentation in the uncorrelated fragmentation contributions (F) and equal weights for the uncorrelated and coalescence contributions. These approximations can be changed and their impacts on predictions can be assessed should experimental data point to intrinsic charm.
Presented here are neutrino fluxes with a single proton interaction, responsible for the bulk of the high energy neutrinos. Secondary and tertiary interactions will produce more, lower energy, neutrinos to increase the event rate. A full Monte Carlo of pion, kaon and charmed hadron production with multiple interactions within the beam dump target is beyond the scope of this work. Ideally, a full simulation can be performed after direct measurements of charmed hadron production are made, for example, by the DsTau experiment. While perturbative QCD has a strong theoretical foundation, the value of the average intrinsic transverse momentum, implemented here with a Gaussian function, is phenomenological. If our default value of k 2 T = 1.5 GeV 2 is increased to k T = 1.5 GeV, the tau neutrino plus antineutrino event rate decreases by a further factor of ∼ 0.82. A smaller value of k 2 T increases the event rate.
In lieu of the next-to-leading order QCD evaluation of charm pair production with additional intrinsic transverse momentum, the Pythia Monte Carlo can be used to model charm production. A disadvantage of Pythia is that it is tuned to central production rather than forward production. Comparisons of Pythia modeling of charm production with the E791 results in π − N collisions with a  500 GeV beam energy [56] showed that the data do not match well with the default tune, as noted in ref. [57]. Dijkstra and Ruf, in ref. [57] use a retuned Pythia and consider multiple interactions within the target to produce charm. They find that the total amount of charm produced is 2.3 times the amount produced in the primary pN interaction, although, with a softer spectrum, so the event rate will not increase by as large a factor.
Neutrino events in a beam dump experiment like SHiP will provide important tests of lepton flavor universality in tau neutrino and antineutrino interactions. Both tau and muon flavors of neutrino and antineutrino event rates, as a function of neutrino energy, will enable a disentangling of perturbative and non-perturbative elements of charm production. The small angular region around the proton beam direction required for neutrinos to intercept the SHiP detector probes very forward charm production, complementary to the higher energy measurements of charm in collider experiments. Better measurements of charm production, either directly or through their neutrino decays, will help narrow uncertainties for neutrino fluxes from charm produced by cosmic ray interactions in the atmosphere, a background to the diffuse flux of astrophysical neutrinos. π, K, pQCD π, K, pQCD + IC25 π, K, pQCD + IC7 π, K, pQCD + 1 10 IC25 Figure 16. Neutrino-antineutrino asymmetries in the event rates, as a function of energy, for muon neutrinos and antineutrinos, as defined in eq. (4.6). Error bars reflect a ±40% factor multiplying the perturbative charm contribution to the asymmetry, as an estimate of the uncertainty due to the scale dependence.

Acknowledgments
This work is supported in part by a US Department of Energy grant DE-SC-0010113. MHR acknowledges discussions with F. Tramontano and collaboration with Y. S. Jeong for the development of the collinear approximated fluxes and event rates based on the HVQ code.

A Intrinsic charm probability distributions
For completeness, we include the differential distributions of the x F distributions for intrinsic charm production in the BHPS model described in detail by Gutierrez and Vogt in ref. [29]. The x F probability distributions depend on an uncorrelated fragmentation function (F) and coalescence distributions (C) that are specific for each charmed hadron produced. We consider incident protons, where the valence quarks together with QQ = cc fluctuations contribute to the 5 quark configuration, and an additional qq = uū, dd, ss pair contributes to 7 quark configurations. As outlined in the appendix of ref. [29] for D ± , D ± S and Λ ± c , and extended to D 0 ,D 0 , the probability distributions are If the sum of charged and neutral D's (but not D ± S and Λ ± c ) for pN interactions gives dσ/dx F (x F = 0.32) = 25 µb, the cross sections are those shown in table 4.

B Neutrino flux at SHiP from pions and kaons
The muon neutrino-antineutrino asymmetry has contributions from pions and kaons that decay to neutrinos since pA production ratios of positively charged mesons to negatively charged pions and kaons is not unity. A full evaluation of the neutrino flux from pions and kaons is beyond the scope of this work, but we estimate their contributions to the electron and muon neutrino fluxes and muon asymmetry using the flux estimation of ref. [1]. In fig. 5.25 of ref. [1], the neutrino fluxes into a detector 56.5 m from the target, with an area of 1.3 m 2 , are shown from an evaluation that includes the full kinematics from D S decays and a GEANT4 evaluation of pion and kaon production by 400 GeV protons incident on a molybdenum target. Their figure is normalized to 100 neutrinos, and it is approximated by the solid curves shown in fig. 14, except for the low energy fluxes. We use our evaluation of the high energy flux of ν τ +ν τ to normalize the flux of neutrinos from pions and kaons in ref. [1], as described below. The tau neutrino flux is approximately φ c ντ +ντ = 4.0 × 10 −2 N ν exp −E/29.4 GeV . (B.1) The tau neutrino plus antineutrino flux determines the overall normalization constant to be approximately N ν = 8.7 × 10 13 neutrinos/GeV, by matching our perturbative φ ντ +ντ with eq. (B.1) at for energies between E 50 GeV and η ν > 5.3. The comparison of our flux evaluation of tau neutrinos plus antineutrinos from perturbative charm production and the normalized curve from eq. (B.1) is shown in fig. 17. Given the flux of ν τ +ν τ , from charm, we can subtract the charm contributions from the total ν e +ν e and ν µ +ν µ fluxes. We take the ratio of the electron neutrino flux at E ν = 200 GeV to the tau neutrino flux at the same energy equal to a rescaling constant to relate the charm contributions to the fluxes φ c νe+νe = φ c νµ+νµ to φ ντ +ντ , since the tau neutrino flux only comes from charm. The high energy contributions of kaons (and pions to ν µ +ν µ ) is negligible due to meson reinteractions. In fig. 14, extrapolated to 200 GeV, the ratio is N c = 8.9. The flux φ c νe+νe = φ c νµ+νµ is shown with the (blue) dot-dashed line.
The difference between the total neutrino fluxes and the charm contributions gives the contributions from pions and kaons. Since pions don't contribution to ν e +ν e , the difference in curves is due only to kaons. After the charm subtraction, then, 3) The kaon contribution to the electron neutrino flux is shown with the dashed blue curve in the figure.
To determine the kaon contribution to the muon neutrino flux, the ratio of charged to neutral kaons is needed. Kaon production is taken in proportions N (K S ) = N (K L ) = (N (K + ) + 3N (K − ))/4 according to an approximation of valence u v = 2d v and sea u s =ū = d s =d = s =s in a simple isospin symmetric model of production [55]. The ratios of positive to negative particle production, in terms of the ratio x R = E * /E * max , the ratio of the energy of the kaon (and for reference, the pion) in the center of momentum frame to the maximum energy, approximately equal to the ratio of the energies in the lab frame, are [55,58] r π (x R ) = 1.05 · (1 + x R ) 2.65 (B.4) The charge ratio for pions saturates around ∼ 5, while the charge ratio for kaons increases with kaon energy [58], which translates to more neutrinos than antineutrinos from pions and kaons. The charge ratios and branching fractions to ν e and ν µ allow for an approximate energy dependent scaling of the kaon contribution to the electron neutrino flux already determined, to the kaon contribution to the muon neutrino flux. The kaon contribution to the muon neutrino flux is shown with the upper (red) dashed line in fig. 14. For x R in eq. (B.5), we have made the approximation that E ν ∼ (1/3 − 1/2)E K to find the energy dependent ratio between electron and muon neutrino fluxes for the dominant decays: K 0 e3 , K 0 µ3 , K + e3 , K + µ3 , and K + → µ + ν µ . The muon neutrino plus antineutrino flux from pion production, with its subsequent decay to π + → µ + ν µ , with a positive to negative charge ratio from eq. (B.4), is shown in fig. 14  The relative normalizations between the kaon contributions to the muon neutrino and electron neutrino fluxes depend on r K via N K (E ν ) = 0.27N K L + 0.64 N K + (2x ν ) + N K − (2x ν ) 0.41N K L + 0.051 N K + (3x ν ) + N K − (3x ν ) (B.10) where x ν = E ν /E p . For kaons, we take the energy fraction carried by the muon neutrino to be on average a factor of 1/2 of the kaon energy, while for the electron neutrino, we use a factor of 1/3, but the results are not very sensitive to this choice. The numerical factors in eq. (B.10) come from the dominant branching fractions to muons and electrons in kaon decays. To determine the asymmetry in muon neutrino and muon antineutrino event rates, we use the overall fluxes from pions and kaons in eqs. (B.8-B.9) and the charge ratios in eqs. (B.4-B.5).