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 charm hadrons that decay promptly, as well as some neutrinos from pion and kaon decays. The short lifetimes of charm hadrons 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 neutrinos and antineutrinos are also produced JHEP02(2019)077 from charm hadron decays. The charged and neutral D mesons decay semileptonically to muon neutrinos and electron neutrinos. Cabibbo and kinematic suppression factors mean that the D ± mesons decay to ν τ in two-body decays with a branching fraction < 1.2 × 10 −3 [7], 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 . Overall, perturbative production of charm results in a flux of muon neutrinos and antineutrinos at a level that is a factor of ∼ 10 larger than the flux of tau neutrinos and antineutrinos.
The SHiP experimental configuration, with a molybdenum target 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-toleading order (NLO) perturbative 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 to approximate whether or not the neutrinos and antineutrinos from charm decays were directed into the detector. After applying a charm momentum angular (or pseudorapidity) cut, the energies of the mesons, taus and eventually neutrinos were traced using analytic energy distributions where relative angles are already integrated.
In this paper, our primary focus is a refined evaluation of the flux of neutrinos and antineutrinos at SHiP that includes the angular dependence of charm hadron decays to neutrinos and intrinsic transverse momentum effects in charm hadron production. We evaluate the forward flux of tau neutrinos plus antineutrinos for SHiP using NLO perturbative QCD production of charm to compare to the evaluation in ref. [2]. In this paper, we also evaluate the flux of neutrinos using a model of intrinsic charm. Models of enhanced charm production (see, e.g., refs. [8][9][10]) including intrinsic charm [11][12][13][14][15][16][17][18][19] have been explored in the context of the high energy atmospheric neutrino flux [20][21][22][23][24]. Using the model of intrinsic charm described in refs. [11][12][13][14] and normalized according to ref. [22] for their evaluation of the atmospheric neutrino flux, we find a significant increase in the number of events from tau neutrinos and antineutrinos at SHiP.
Intrinsic charm introduces particle-antiparticle asymmetries at SHiP that we evaluate here for muon neutrino and antineutrino events. Gluon fusion, gg → ccX, is the dominant parton model process in perturbative production of charm [25][26][27][28][29], so particle-antiparticle asymmetries are not significant. On the other hand, for five or seven particle states (e.g., |uudcc , |uudccuū ) where valence quark momenta are brought to, for example, D 0 = uc, but not a D 0 =ūc, particle-antiparticle asymmetries can be large. They are significant for D + /D − and D 0 /D 0 , but considerably less significant for D + s /D − s . Thus, asymmetries in the muon neutrino and antineutrino fluxes, but not the tau neutrino and antineutrino fluxes, are predicted with intrinsic charm. We make an approximate evaluation of the muon neutrino-antineutrino asymmetries relying on our perturbative charm evaluation, intrinsic charm and an approximate flux of muon neutrinos and antineutrinos at SHiP from pion and kaon decays based on fluxes evaluated in ref. [1].

JHEP02(2019)077
We begin in section 2 with the production of charm quarks and hadrons in the parton model and with a model for intrinsic charm based on refs. [11][12][13] discussed in detail in ref. [14]. Our assumptions about hadronization and the incorporation of intrinsic transverse momentum are included in section 2, where we compare our NLO perturbative QCD transverse momentum distribution of D + + D − + D 0 +D 0 in pp collisions with E p = 400 GeV to LEBC-EHS measurements [30]. We show energy and rapidity distributions of charm hadrons to illustrate some of the features of charm production that translate through charm hadron decays to the neutrino energy and pseudorapidity distributions, at least approximately. Efforts to measure production of charm hadrons with the 400 GeV beam at CERN, e.g., the DsTau project [31], would help refine the perturbative calculation of the neutrino flux and constrain forward production of charm hadrons from intrinsic charm directly.
In section 3, we show our results for neutrino energy distributions from prompt decays of charm hadrons produced by a beam of 400 GeV protons directed to the SHiP detector. Given N p = 2 × 10 20 protons and a SHiP detector represented by a column depth of 66.7 g/cm 2 of lead, we evaluate the number of events. 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 the number of tau neutrino plus antineutrino events is 10% lower than the number of events determined using the more detailed detector geometry of ref. [2]. Given this good agreement, all of the neutrino and antineutrino results presented here approximate the detector geometry by requiring η ν > 5.3. The primary impact of the calculational improvements discussed here is to reduce the number of neutrinos incident on the SHiP detector. We discuss the impact of each correction and their associated uncertainties in section 3.
We also discuss how an enhanced number of tau neutrino and antineutrino events and their energy dependence could show evidence of intrinsic charm. As noted above, intrinsic charm models predict particle-antiparticle asymmetries that translate to different muon neutrino and antineutrino fluxes. Section 3 also includes an evaluation of the asymmetry in the number of muon neutrino and antineutrino events that shows how intrinsic charm changes the perturbative prediction for charm. With an approximate number of muon neutrinos and antineutrinos from pion and kaon decays, we show the muon charge asymmetry at SHiP as a function of energy. The discussion in section 4 summarizes our results for the number of events for tau neutrinos and antineutrinos, and the muon charge asymmetry, at SHiP. We point to uncertainties in our evaluations and to where refinements to the calculation can be made.
Relevant formulas for intrinsic charm are collected in appendix A. Appendix B gives some details of the tau neutrino energy distributions that come from the decays of the D s directly to ν τ and the chain decay D s → τ → ν τ . Our treatment of the muon neutrino and antineutrino fluxes from pion and kaon decays needed for the muon charge asymmetry, based on the results of ref. [1], appears in appendix C.

Production of charm quarks and hadrons
Charm quark production in perturbation theory is evaluated using the collinear parton model. The factorization scale µ F characterizes this collinear approximation, with the JHEP02(2019)077 interpretation that, roughly, partons with 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 charm hadron production at SHiP.

Perturbative charm
The parton model evaluation of charm production involves the production of a charm quark pair, followed by the hadronization of the quark and antiquark into a charm 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 , 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) [25][26][27][28][29].
The differential cross section for the inclusive c cross section for charm momentum p is for parton cross sectionσ. Implicit in eq. (2.1) are integrals over thec momentum for the leading-order and NLO 2 → 2 processes, and over thec and quark or gluon momentum in the 2 → 3 processes. Explicit expressions for the one-particle inclusive heavy quark differential cross section, after these integrals are performed, appear in ref. [26]. At SHiP, a proton beam, H 1 = p, is incident on a molybdenum target, so H 2 =Mo. The molybdenum atomic mass is A = 96 and charge Z = 42. The cross section for charm production is dominated by gg interactions which are independent of the number of neutrons and protons. There are potential valence quark effects in the forward production of charm, however. We approximate the parton distribution function (PDF) for partons in the nucleus by for charm production, where we assume isospin symmetry to relate the up and down quark distributions in the neutron to those in the proton. Theoretical evaluations of perturbative production of charm have a long history, with evaluations through NLO in QCD in refs. [25][26][27] that have been implemented in the HVQ computer program with the HVQMNRPHO main program and supporting fortran subroutines. Our perturbative evaluation uses the HVQ computer program as a basis for the NLO evaluation. We use the single inclusive charm quark evaluation provided there,

JHEP02(2019)077
with the addition of intrinsic transverse momentum and an implementation of the D s decay D s → τ ν τ and subsequent tau decays, as described below. For the main 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 ref. [32] where m c = 1.27 ± 0.09 GeV is used. The scale choices with m c = 1.27 GeV permit us to make a direct comparison with the results in ref. [2]. With m c = 1.27 GeV, we can also take advantage of the detailed data-constrained analysis of charm production over a wide energy range performed in ref. [32] to set the range of scales to be µ F = 1.25 − 4.65 m c and µ R = 1.5 − 1.7 m c .
We use the nCTEQ15 PDFs [33], available for both free and bound nucleons. The nCTEQ15 distribution of grids does not include Mo PDFs, but the PDFs for Kr (A = 84) and Ag (A = 108) give essentially the same cc cross section for E p = 400 GeV, where we use eq. (2.2) for Mo.
Before we discuss the nuclear PDF corrections, we note that the charm pair production cross section per nucleon at next to leading order with µ F = 2.1 m c , µ R = 1.6 m c and m c = 1.27 GeV for a proton beam energy of E p = 400 GeV incident on a Mo target is 16.3 µb using free nucleon nCTEQ15 PDFs, whether one includes the specific Z and A dependence for Mo or uses an isoscalar nucleon target. The nCTEQ15 free nucleon PDFs give a charm pair production cross section that is 14% lower than the cross section using the CT14 PDFs [34] for the same charm quark mass and scale choices. The measured cross section for pp scattering with E p = 400 GeV is 18.1 ± 1.7 µb as determined by Lourenco and Wohri in ref. [35] from data reported in ref. [30]. For our range of renormalization and factorization scales, the pp cross section with a 400 GeV beam energy is σ cc = 13.3 − 19.1 µb, with the lower cross section within 3σ and a central value of 16.2 µb that is about 12% lower than the central value of the measured pp charm pair cross section at this energy.
In general, it is not the case that the upper and lower theoretical bounds on the cross section are obtained by taking both the lower renormalization and factorization scales, or both the upper scale choices. However, for E p = 400 GeV, it is the case. At this energy, the lower renormalization scale gives a larger cross section because of a larger α s . Larger factorization scales result in PDF evolution that pushes the PDFs to smaller x values. For E p = 400 GeV, the larger factorization scale reduces the PDFs in the x region of interest. At higher incident beam energies, the evolved enhanced small x PDF distributions are more important than they are for E p = 400 GeV.
Distributions are often evaluated using the transverse mass, m T = p 2 T + m 2 c . In our single inclusive distributions, we can evaluate the cross section with this scale choice. With m c = 1.27 GeV, using µ F = 2.1m T and µ R = 1.6m T , the pp cross section is σ cc = 13.3 µb, on the lower side of the cross section range. In the results presented here, we use a default scale dependence on the charm quark mass rather than the transverse mass. Below, we use data from ref. [30] to set the parameter in a Gaussian intrinsic transverse momentum distribution. A different parameter for intrinsic transverse momentum would be obtained if the transverse mass scale is used instead.
At high energies where very small x PDFs are important, nuclear shadowing can have a large impact on forward production, even for fairly low A targets like nitrogen, as shown JHEP02(2019)077 in the context of the atmospheric neutrino flux from charm [36]. At the energy considered here, there are two largely compensating effects: the shadowing at smaller x values and the EMC effect at large x are offset by anti-shadowing for x in the range of ∼ 0.01 − 0.3 (see, e.g., refs. [34,37,38] and references therein). The net effect is that charm pair cross section is enhanced with nuclear corrections. The cross section is 18.5 µb for bound nucleons in Mo using the central PDF of nCTEQ15 for Ag, so an increase of about 13% due to nuclear corrections. By varying the nuclear PDFs over the 33 sets of nCTEQ15, the cross section on Mo ranges from 17.8 µb-19.4 µb, about ±5%. For reference, the gg initial state accounts for 94% of the cross section on bound nucleons for the central PDF choice in molybdenum for E p = 400 GeV.
The charm hadron and anti-hadron production rates are essentially equal for perturbative production of charm. This will distinguish perturbative charm production from intrinsic charm production, discussed below.

Intrinsic transverse momentum and fragmentation
For the large rapidities required by SHiP, intrinsic transverse momentum can deflect charm quarks from the detector direction, so it should be included. The fragmentation of charm quarks to charm hadrons also influences the transverse momentum distributions, as was noted, for example, already in ref. [39] and more recently for LHC energies in ref. [40]. The intrinsic transverse momentum depends on √ s [32]. The LEBC-EHS experiment [30] with 400 GeV protons incident on a hydrogen target is our benchmark for setting the intrinsic transverse momentum distribution value given the fragmentation functions used here. Our default values for fragmentation parameters, the same as in ref. [2], allow us to compare the intrinsic transverse momentum improved evaluation with the collinear approximated evaluation [2].
For hadrons H equal to charged and neutral D's and for Λ c , we can write according to the Peterson fragmentation form [41], where we take z as the ratio of the momentum of the meson to the momentum of the charm quark, p H = z p c , in the lab frame. We use H from Kniehl and Kramer [42], with µ 0 = 1.5 GeV, and we rescale N H to normalize the fragmentation functions to the fragmentation fractions of ref. [43]. The parameters H and the fragmentation fractions shown in table 1 are used as default values for the results shown here. We do not evolve the fragmentation functions. The Peterson fragmentation functions are used for all of the distributions of charm hadrons and neutrinos evaluated using NLO perturbative QCD presented here, unless specifically noted. Schweitzer, Teckentrup and Metz [44] showed that a Gaussian model of intrinsic transverse momentum is compatible with semi-inclusive deep-inelastic scattering and with Drell-Yan production. In the Gaussian approach, azimuthal symmetry around the beam axis permits us to write

JHEP02(2019)077
normalized to unity with the integral over d 2 k T . For the Gaussian distribution, (2.5) The differential distribution for charm quark production becomes In ref. [44], it was shown that for leading order Drell-Yan production, the Gaussian approximation provides a good model if for a 300 GeV proton beam scattering on platinum (FNAL-288) [45]. For pion scattering from tungsten (FNAL E615), the value is a little larger (1.7 GeV 2 ) [46]. Using the Peterson fragmentation function for c → D ± , D 0 andD 0 with D = 0.10, our NLO evaluations with k 2 T = 1.1 − 1.7 GeV 2 provide good fits with the D meson p 2 T distribution measured by the LEBC-EHS group [30]. The best fit with the data has k 2 T = 1.4 GeV 2 . The best fit distribution (solid red histogram) and the range of predictions for k 2 T = 1.1 − 1.7 GeV 2 (dashed black histograms) with the LEBC-EHS data are shown in figure 1.
We find a similar range of k 2 T values when we compare our predictions with the charm meson data from a 250 GeV proton beam incident on copper [47] when we make a 10% correction to the overall normalization of the p 2 T distribution. To illustrate the separate fragmentation and intrinsic transverse momentum effects, we also show in figure 1 the charm distribution without fragmentation or intrinsic transverse momentum (solid blue histogram), and with intrinsic transverse momentum but without fragmentation (dashed orange histogram). The green dashed histogram shows the result of including fragmentation, but not intrinsic transverse momentum.
In the evaluations of charm hadron production and the resulting neutrino energy and rapidity distributions, we use k 2 T = 1.4 GeV 2 for the average intrinsic transverse momentum squared. This is equivalent to k T 1.05 GeV. We discuss the impact of the range of acceptable values of k 2 T on our predicted number of tau neutrino plus antineutrino events at SHiP in section 4.
The choice of D = 0.10 for D 0 and D + production is larger than a more typical choice of D = 0.06. With our default value, we can make direct comparisons with the muon neutrino results evaluated in the collinear approximation in ref. [2]. As noted above, the choice of k 2 T depends on the fragmentation function. With D = 0.06, k 2 T = 0.9 − 1.4 GeV 2 gives a similar band of predictions to the choice of D = 0.10, k 2 T = 1.1 − 1.7 GeV 2 . Intrinsic transverse momentum choices outside the respective ranges of values are problematic for the larger values of p 2 T shown in figure 1, despite the large error bars that span a factor of ∼ 2 in the highest p 2 T bin. The large errors in the highest p 2 T bins are not reflected proportionately in the error estimates for the neutrino events because the detector angular size is small. At low p 2 T , the experimental error bars are of order ±15%, while the spread in predictions with this fragmentation, for the range of k 2 T , is less than ±10%.

Intrinsic charm
Beam dump experiments have the potential to constrain additional contributions to charm hadron production, beyond perturbation theory, through neutrino measurements. The description of high momentum charm hadrons have long been a topic of interest, e.g., as discussed in refs. [8][9][10][11][12][13][14][15][16][17][18][19]. Often cited are the data from the LEBC-MPS collaboration from 800 GeV protons incident on a bubble chamber [48], and the SELEX (E781) collaboration's measurements of high momentum Λ c 's in π − , Σ − and p beams incident on copper and carbon targets [49]. 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.

JHEP02(2019)077
For definiteness to find charm contributions to neutrino fluxes at SHiP, we use an intrinsic charm model of Brodsky, Hoyer, Peterson and Sakai (BHPS) [11,12], described in detail in refs. [13,14,50]. 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 or heavy quarks can independently fragment into high momentum charm hadrons. While this is not the dominant charm hadron production mechanism, it may dominate the large Feynman x F distribution. Our goal here is not to provide a definitive model of IC. We use this model to illustrate the potential for SHiP to explore the contributions of IC, without introducing specific fragmentation functions, binding energy and mass effects beyond what is outlined below.
Based on frame-independent probability distributions for intrinsic charm P ic , two types of contributions to the production of charm hadrons are included: independent (uncorrelated) fragmentation (F), and coalescence distributions (C) specific to individual charm hadrons. They depend on the n-particle probability functions [14], 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. [14], For coalescence, the probability function includes a delta function that combines the momentum fractions of the quarks that make up the final state hadron, In eqs. (2.9) and (2.10), x H is the Feynman x F variable for the given hadron H. In the large where σ in pN is the inelastic pN cross section and µ 2 /4m 2 c sets the scale of the intrinsic charm cross section. In what follows, we approximate x F x E in all of our evaluations of the neutrinos from intrinsic charm.
Alternatives to the delta functions in eqs. (2.9) and (2.10) can be used to relate the quark momentum to the hadron momentum, for example, the Peterson fragmentation function can be used for the uncorrelated fragmentation. We use the delta functions for intrinsic charm to be consistent in how we include the uncorrelated and coalescence contributions. There are many more uncertainties in the theoretical treatment of intrinsic charm than for perturbative production of charm. A more detailed treatment of fragmentation, binding JHEP02(2019)077 energy, mass effects and other approximations in this BHPS model of intrinsic charm would be required should evidence of intrinsic charm be observed experimentally.
The relative contributions of the five quark (uudcc) and seven quark (uudccqq) terms for qq = uū, dd, ss are set according to ref. [14]. The normalizations are P 5 ic = 3.1 × 10 −3 , and ic and m 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 charm 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 [20][21][22]. 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 [20-24, 36, 51-56]. The IceCube measurements of the diffuse astrophysical flux have ruled out an extremely large intrinsic charm component, as discussed in, e.g., ref. [21]. 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. [22] 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's normalization of the intrinsic charm cross section is based on the highest x F data in pp collisions for E b = 800 GeV ( √ s = 39 GeV) measured by the LEBC-MPS collaboration [48]. Using this highest bin at x F = 0.32, Laha and Brodsky choose [22]: a value above the central value in the bin, but within the error bar [48]. 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, also considered in ref. [22]. The IC7 normalization is such that the differential cross section is at the lower end of the LEBC-MPS error bar in the x F = 0.32 bin for E b = 800 GeV [48]. The intrinsic charm cross section in this approach is expected to scale with energy as the pp total cross section does, so it depends only weakly on energy. We use these same normalizations (IC25 and IC7) for our evaluation of the intrinsic charm contribution to the neutrino flux at SHiP where E p = 400 GeV. Intrinsic charm added with the IC25 normalization is still within the error bars for the lowest and highest p 2 T shown in figure 1, but the fit to the data with k 2 T = 1.4 GeV 2 is not as good as with the perturbative contribution alone. A larger intrinsic charm normalization would over-produce D mesons

JHEP02(2019)077
in the lowest p 2 T bin. The sum of the IC7 plus perturbative production of D mesons as a function of p 2 T has a similar fit to the data shown in figure 1. 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 or for the prompt atmospheric flux at IceCube. The lower beam energy at SHiP measurements is an advantage for constraining the intrinsic charm component.
To account for nuclear effects, we use the appropriate scaling of A α where α = 0.71 [57]. For molybdenum, A 0.71 /A = 0.27. We note that a change in the value of α can be absorbed in a change to the overall normalization of the intrinsic charm contribution. For example, α = 0.67 gives an IC contribution that is ∼ 80% of the IC contribution with α = 0.71. Experimental evidence of intrinsic charm would require a more detailed study of its A dependence.
Asymmetries between D − and D + , and between D 0 andD 0 , for the x F x E 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 antineutrino flux from pion and kaon decay for the SHiP target configuration based on fluxes presented in ref. [1], with normalizations outlined in appendix C. We compare them with contributions from intrinsic and perturbative charm.

Charm hadron production
In this section, we begin with our results for charm hadron production via the perturbative approach for E p = 400 GeV. Since we are interested in the neutrino and antineutrino energy distributions, most of the figures in this section display the charm hadron energy distributions to show the impact of elements of the calculation that will eventually translate to the lepton energy distributions. They also serve as a reference for measurements of charm hadron production, e.g., the DsTau project [31]. Figure 2 (upper) shows the differential cross section as a function of the charm hadron energy for each of the charm hadrons considered here. To compare the distributions, the distributions have been normalized by 1/B c→H . Figure 2 (lower) shows the ratios of the normalized distributions relative to the D − s distribution as a function of charm hadron energy. The D 0 and D − histograms overlap.
The charm hadron energy distributions integrated over the full phase space do not change significantly with the introduction of intrinsic transverse momentum, however, the energy distributions of the charm hadrons do change when an angular cut is applied. To illustrate, we show in figure 3 the energy distribution of the D s when η Ds > 5.3. For η Ds > 5.3, the cross section for D s production including intrinsic transverse momentum is 7% of the total D s cross section. The intrinsic k T has a significant impact. Using the pseudorapidity cut η > 5. 3  approximately a factor of two higher than with intrinsic k T . The inclusion of intrinsic k T 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, the requirement θ 0.01 rad means the energies out to ∼ 100 GeV are impacted. distribution of charm hadrons for hadron energies larger than E ∼ 50−100 GeV, depending on the normalization of the intrinsic charm contribution. Figure 4 shows distributions from intrinsic charm with the IC25 normalization of eq. (2.13). The lower figure shows that intrinsic charm clearly dominates at high rapidity, as anticipated. When a pseudorapidity cut of η > 5.3 is applied to the D s distribution from intrinsic charm, the cross section is 40% of the cross section for all of phase space.
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 charm quark compared to charm antiquark. The Λ + c are produced at the highest rapidities, as can be seen in figure 5, where perturbative and intrinsic contributions to charm hadron rapidities are shown on a linear scale.
Perturbative charm and intrinsic charm production of 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 figure 6, where we have used the tungsten PDFs from nCTEQ15 [33] in the evaluation. The lower blue dotted histogram comes from the perturbative QCD evaluation with intrinsic transverse momentum k 2 T = 1.4 GeV 2 , while the orange dotted histogram shows the intrinsic charm contribution to D ± s normalized by eq. (2.13) (IC25). The sum of the two is shown with the solid green histogram. The red dashed histogram shows the sum with the intrinsic charm reduced by a factor of 7/25, labeled by IC7. For perturbative QCD production, with a branching fraction of B(D s → τ ν τ ) = 5.5% and the factor of 20% for detector efficiency, the power is n 10.0 at large x E for the perturbative result and n 5.0 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. [14].

JHEP02(2019)077 3 Neutrinos and antineutrinos at SHiP
In this section we show energy distributions of neutrino and antineutrinos from charm hadron decays. In section 3.1, the NLO perturbative QCD contributions are evaluated first, followed by the inclusion of the much more uncertain component from intrinsic charm. The neutrino energy distributions from charm hadron decays are converted to the number of events per GeV at SHiP for muon and tau neutrinos and antineutrinos. We represent the SHiP detector as a column depth L = 66.7 g/cm 2 of lead and convert the neutrino energy distributions to events using the neutrino charged current cross section. These results are shown in section 3.2. In section 3.3, we describe an approximation for the number of ν µ andν µ from pion and kaon decays (with more details in appendix C). We evaluate the muon charge asymmetry from muon neutrino and antineutrino charged current interactions, where the muon neutrino and antineutrino fluxes come from pion, kaon and charm hadron decays, where the charm contribution is computed both with and without intrinsic charm.

Neutrino and antineutrino energy distributions from prompt decays of charm
The neutrinos from prompt decays of charm hadrons at the SHiP detector, including the full three-dimensional kinematics, are shown in this section. The D s decay D s → τ ν τ , with the τ → ν τ X that depends on the tau polarization, is an interesting feature. We refer to the tau neutrino produced directly in the D s → τ ν τ as the "direct neutrino" and the neutrino that comes from the tau decay as the "chain decay neutrino." For reference, we include in appendix B some details of our numerical implementation of the production of direct and chain decay tau neutrinos from D s decays. The NLO QCD results for the neutrino differential cross section as a function of neutrino energy are shown in figure 7 for E p = 400 GeV. The distributions are the sum of neutrinos and antineutrinos.
The neutrino or antineutrino rapidity is required to be η ν > 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. The direct neutrino has a much smaller energy than the tau in the D s rest frame. In the lab frame, even though the neutrino from the tau decay D s → τ → ν τ has a fraction of the tau energy, itself a fraction of the D s energy, the chain decay neutrino has more energy on average than the direct neutrino. The cross-over 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 charm 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  figure 9 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.4 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 charm hadrons D − andD 0 are favored over their charged conjugates with this model for intrinsic charm momentum distributions, the production rate for muon anti-neutrinos is larger than for muon neutrinos. Above E ∼ 20 GeV, the asymmetry between particle and antiparticle perturbative production is modified by intrinsic charm. It will be apparent if the intrinsic charm contribution is normalized by IC25 (eq. (2.13).) 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 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 a much smaller fraction of the perturbative charm cross section.  JHEP02(2019)077

SHiP neutrino events from prompt decays of charm
The conversion of the differential neutrino energy distribution to events with a charged lepton requires the neutrino charged current cross section. The SHiP detector will have thin emulsion films interleaved with lead bricks [1]. For the flux of neutrinos from charm decays, we write As in ref. [2], we take N p = 2 × 10 20 protons and hadronic cross section per nucleon σ pA /A = σ pN = 10.7 mb. The number of events per GeV is where λ ν = (N A σ(νA)/A) −1 and L = 66.7 g/cm 2 for lead in the detector.
To evaluate the neutrino cross sections with lead relevant for SHiP, we use the nCTEQ15 parton distribution functions for lead [34]. For tau neutrino interactions, we include the m τ /E ν corrections to the differential cross section [58][59][60][61]. We use the low momentum transfer extrapolation of the deep-inelastic scattering (DIS) structure functions outlined in refs. [61,62]. Shown in figure 10 are the neutrino and antineutrino chargedcurrent 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 figure 10 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. For reference, table 2 lists the charged current cross sections for scattering with a lead target by muon neutrinos and antineutrinos, and tau neutrinos and antineutrinos for energies up to E ν = 400 GeV. This is beyond the energy range of neutrinos at SHiP, but it is potentially of interest for E p = 800 GeV protons incident on a beam dump. The muon neutrino cross section is still 10% larger than the tau neutrino cross section for E ν = 400 GeV, even though E ν m τ . The uncertainties in the cross sections for neutrino and antineutrinos incident on lead can be approximated by considering the PDF uncertainties as quantified by the 32 variants to the nCTEQ15 PDFs for lead targets. The cross section uncertainties for ν τ andν τ DIS interactions in lead using these PDF variants are smaller than ±3% for the energies considered here, except for the lowest energies, below ∼ 10 GeV. Figure 11 shows the impact of intrinsic transverse momentum and the full threedimensional 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 per GeV for tau neutrinos plus antineutrinos, evaluated using most of the same approximations as in ref. [2]: that the differential cross section for charm quark production, evaluated with CT14 PDFs with no intrinsic transverse momentum but with fragmentation, in the approximation that the charm quark momentum JHEP02(2019)077 is in the same direction relative to the incident proton beam as the neutrino from its decay (the "collinear approximation" that η ν η c ). Here we show for η ν > 5.3 as an approximate representation of the detector solid angle. Using CT14 PDFs reduces the number of events by ∼ 10% compared to the number of events with CT10 PDFs used in ref. [2].
The three lower histograms in figure 11 are evaluated using the best fit nCTEQ15 PDFs for the free proton incident on the fixed molybdenum target. With only the addition of intrinsic transverse momentum k 2 T = 1.4 GeV 2 , the number of events is 62% of the number of events without intrinsic transverse momentum. The three dimensional decay kinematics further reduces the overall number of events. The number of events goes from 939 to 272 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 τ ν τ .
Our main results for tau neutrinos and antineutrinos at SHiP, evaluated with NLO perturbative QCD, are shown in figure 12 The number of events vary by less than ∼ ±40% below E = 100 GeV, less than ∼ ±30% below 50 GeV. Figure 12 also shows a broader uncertainty band for the separate tau neutrino and antineutrino number of events per GeV (upper plot) and for the sum (lower plot). The full uncertainty band includes, in addition to the scale dependence and choice of k 2 T , a number of other uncertainties added in quadrature. The 33 sets of nuclear PDFs provided in the nCTEQ15 distribution translate to an uncertainty in the number of tau neutrinos plus antineutrinos that increase with energy, from a few percent error at low energies, to ∼ 10% for E ν = 50 GeV and ∼ 20% for E ν = 100 GeV. Using the nCTEQ15 set for krypton rather than silver decreases the distribution as a function of energy with a smaller, but similar correction. The nCTEQ15 results are about 14% lower than using CT14 PDFs for the free proton, so we add (in quadrature) an additional 14% uncertainty to the upper band. The uncertainty band also includes decreases in the distribution that come from using m c = 1.4 GeV and µ R,F ∼ m T . Each of these two changes to the inputs of the (blue) histogram is the sum of tau neutrino plus antineutrino events. The total number is 1985 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 at x F = 0.32 (IC7), the number of events 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. A similar double peak structure is seen in the NLO perturbative QCD evaluation for η ν > 6.0, but the peaks are overlapping for η ν > 5.3, shown in figure 12. Again, 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. The number of events from tau neutrinos and antineutrinos with η ν > 5.3 and interactions in the SHiP detector, modeled as a lead detector with column depth L = 66.7 g/cm 2 , from a beam of 2 × 10 20 protons with E p = 400 GeV incident on molybdenum are summarized in table 3. We do not evaluate a detailed uncertainty for intrinsic charm here because the intrinsic charm contribution to forward charm production is only mildly constrained. As remarked in section 2, the A dependence uncertainty can be absorbed into the overall normalization. More detailed analyses with fragmentation functions, binding energy and mass effects would be merited in the event that an intrinsic charm effect is measured.
We now turn to the number of muon neutrino and antineutrino events. We begin with NLO perturbative QCD predictions, including intrinsic transverse momentum and threedimensional kinematics, for muon neutrinos and antineutrinos from charm production and decay. Figure 14    ν µ (IC25) ν µ (IC25) ν µ (π, K, pQCD) ν µ (π, K, pQCD) charm as a function of energy for the 400 GeV proton beam with the standard parameter choices. The two dashed histograms show the intrinsic charm number of events per GeV. While for the perturbative evaluation, the difference in neutrino and antineutrino number of events per GeV is due to the different neutrino and antineutrino cross sections, the number of events per GeV for intrinsic charm reflects enhanced forward production of hadrons with a charm 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. For the perturbative evaluation of charm, the number of events from muon neutrino and antineutrino interactions at SHiP with η ν > 5.3 are a factor of 20 larger than the number of events from tau neutrino plus antineutrino interactions, a total of 5738 events for the sum of solid histograms (labeled pQCD) in figure 14. The intrinsic charm IC25 number of events for muon neutrinos and antineutrinos, a total of 21055 for the sum of dashed histograms in figure 14, is a factor of ∼ 10 larger than the number of events for tau neutrinos and antineutrinos from the IC25 evaluation. Where there is no difference in the tau neutrino and antineutrino fluxes, the fraction ofν τ events of the total number of ν τ +ν τ events is about 30%, as is the fraction ofν µ events to the total number of ν µ +ν µ events evaluated from perturbative QCD. The fraction ofν µ to the total of ν µ +ν µ from IC25 alone is 44%.

Muon charge asymmetry at SHiP
The asymmetry in particle and antiparticle production with intrinsic charm is an interesting feature to use as a diagnostic for intrinsic charm. As noted above, the asymmetry in charm hadron particle-antiparticle production will translate to the numbers of ν µ andν µ per GeV, but not for ν τ andν τ . 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 appendix C. 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 [63] are used to approximate the separate neutrino and antineutrino fluxes. Figure 15 shows these fluxes with the separate pion, kaon and charm distributions. 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 so we have shaded the region of E ν < 40 GeV in figure 15.
The number of muon neutrino and antineutrino events per GeV from pQCD production of charm combined with our modeled pion and kaon neutrinos is shown in figure 16 with the green histograms: the upper solid histogram for ν µ charged current events per GeV and lower dashed green histogram forν µ events per GeV (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ν µ number of events per GeV. Since we have not done a full error analysis of the D → ν µ modeling, we have expanded the error band for the perturbative contributions to ν µ +ν µ . The uncertainty in the total number of events for the perturbative QCD evaluation of the ν τ +ν τ events in table 3 is 0.5-1.4 times the central value. The shaded bands reflect a more conservative error estimate for the perturbative muon neutrino and antineutrino fluxes, taking a factor of 0.3-2 times the central perturbative result. The green error band reflecting charm production of ν µ decreases with energy because at higher energies, the kaon contribution dominates. The green error band reflecting charm production ofν µ increases with energy as the charm contribution increases in importance with increasing energy.
The nearly flat orange histograms in figure 16 show the ν µ (solid) andν µ (dashed) number of events per GeV where they come from only intrinsic charm production. The intrinsic charm events per GeV for ν µ 's is approximately equal to that of muon neutrinos from the combined contributions from pions, kaons and perturbative charm, and forν µ from JHEP02(2019)077 intrinsic charm, all much larger than the π, K and pQCD contributions to the number of ν µ events per GeV. One can define the charge asymmetry as where dN/dE = φ evts (E) for the respective interactions. Figure 17 shows the charge asymmetry defined in eq. (3.3). 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.13) for protons incident on molybdenum (IC25) added to the perturbative charm, pion and kaon contributions. Intermediate curves show the results for IC7 (eq. (2.14)) and when the IC25 normalization is decreased by a factor of 1/10. The error band shown with the dashed histogram for π, K and perturbative charm production of muon neutrinos and antineutrinos shows an approximate energy independent ±40% error on the perturbative ν µ +ν µ flux that comes from scale dependence. The shaded bands show a conservative error estimate of the perturbative error with a range of 0.3-2.0 times the central perturbative ν µ +ν µ result. Even with the larger error band on the perturbative contributions, the change in the muon charge asymmetry when intrinsic charm contributes at the level of the IC25 or IC7 normalization is evident. A full study of perturbative and nonperturbative contributions to the ν µ andν µ flux is beyond the scope of this work. However, with our approximations for the pion and kaon contributions to the number of ν µ +ν µ , it is our estimate that the JHEP02(2019)077 number of intrinsic charm events normalized according to IC7 or larger will be reflected in a muon charge asymmetry different from the results from pion, kaon and perturbative charm contributions to that asymmetry.

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 charm 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 number of tau events by a factor of ∼ 3 compared to the number of events in ref. [2], predicting 272 events with our central renormalization and factorization scale choices and k 2 T = 1.4 GeV 2 . Scale and intrinsic transverse momentum uncertainties introduce a ∼ ±30 − 40% uncertainty in the number of events. A larger charm quark mass or scale choice that depends on m T rather than m c results in a smaller number of events. Combining mass, scale factor and other uncertainties leads to a lower value of 134 events and an upper value of 378 events, as summarized in table 3. The wide range of perturbative QCD predictions for the number of tau neutrino events for η ν > 5.3, here estimated to be a factor of 2.8, emphasizes the importance of additional experimental measurements of charm productions, for example, with the proposed DsTau experiment [31].
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 [11][12][13][14] as a representative model to show the impact that forward charm production can have on the number of tau neutrino and antineutrino events. 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 number of events. Indeed, intrinsic charm is relatively more important at low energies, e.g., at SHiP, than for high energy cosmic ray interactions in the atmosphere because the cross section is assumed to scale like the total hadronic cross section, so it is weakly dependent on energy, unlike the perturbative charm cross section. With the higher intrinsic charm normalization suggested in ref. [22] (IC25), the number of events for tau production from intrinsic charm alone is almost 2,000 tau neutrino plus antineutrino charged current events. The energy distribution of the number of tau neutrino and antineutrino events, different for perturbative production of charm and for intrinsic charm, will be essential to untangle the their relative contributions given the uncertainties in the total number of events.
Beyond the enhancement in the number of 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 different from the NLO perturbative QCD predictions JHEP02(2019)077 won't appear in the tau neutrino and tau antineutrino events, but they will appear in muon neutrino and antineutrino charged current events. 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 [63] 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 muon neutrino and antineutrino fluxes.
We have used the BHPS model for intrinsic charm to compare with what has been used to evaluate prompt atmospheric neutrino fluxes in ref. [22], however, this is not the only model with forward charm production or for intrinsic charm. Within the BHPS model, we have chosen simplified fragmentation in the uncorrelated fragmentation and coalescence contributions and equal weights for the two contributions. These approximations can be changed and refinements to the treatment of fragmentation can be made. 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 number of events. A full Monte Carlo of pion, kaon and charm hadron production with multiple interactions within the beam dump target is beyond the scope of this work. Ideally, a full simulation would be performed and compared to direct measurements of charm hadron production.
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 [64] showed that the data do not match well with the default tune, as noted in ref. [65]. Dijkstra and Ruf, in ref. [65] 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 number of events 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 number of events, 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. The high pseudorapidity and relatively low beam energy of SHiP will permit stronger constraints on models of intrinsic charm, with their normalizations only weakly dependent on energy, than at high energy colliders. 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.

JHEP02(2019)077
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.

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. [14]. The x F probability distributions depend on an uncorrelated fragmentation function (F) and coalescence distributions (C) that are specific for each charm 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. [14] for D ± , D ± s and Λ ± c , and extended to D 0 ,D 0 , the probability distributions are

B Charm hadron decays to neutrinos
Tau neutrinos and antineutrinos come from D s production and their decay D + s → τ + ν τ (direct neutrino) and tau decay τ + →ν τ X (chain decay neutrino), and charge conjugate processes. The branching fraction is B(D s → τ ν τ ) = (5.55 ± 0.23)% [7]. We have implemented the D ± s and τ ± decay distributions into the HVQ computer code. We sketch some of the inputs here. The direct and chain decay neutrino energy distributions from π ± and µ ± decays are discussed in detail in refs. [66][67][68]. The evaluation is directly transferable to D s → τ ν τ .
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 Ds ) 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 boosts which we evaluate numerically. Ref. [69] has a discussion of the details of the boosts, for example. The ν τ 's in the two-body D s decays are the direct neutrinos.
To get the energy and momentum of the chain decay neutrinos, we begin with the tau in its rest frame. In the tau rest frame, the distribution of the massless (or nearly massless) leptons from the three-body decay of the tau is given by [66][67][68] 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.

JHEP02(2019)077
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 [66,67] cos θ * Ds = where E Ds and E τ are the energies of D s and τ in the lab frame. The results shown below use eq. (B.3) rotated and 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. [36]. 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 charm hadron mass and m K is the effective mass of the hadronic final state in the decay.

C 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 figure 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 figure 15, 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 with eq. (C.1) for energies between E 50 − 100 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. (C.1) is shown in figure 18.
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 figure 15, 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, 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 = u = d s =d = s =s in a simple isospin symmetric model of production [63]. 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 JHEP02(2019)077 to the maximum energy, approximately equal to the ratio of the energies in the lab frame, are [63,70] r π (x R ) = 1.05 · (1 + x R ) 2. 65 (C.4) The charge ratio for pions saturates around ∼ 5, while the charge ratio for kaons increases with kaon energy [70], 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 figure 15. For x R in eq. (C.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. (C.4), is shown in figure 15 with the steeply falling (red) dot-dashed line.
The curves shown in figure 15, with the overall normalization, N ν = 7.1 × 10 13 neutrinos/GeV, can be parameterized by (see also eq. (C.1)): φ c νµ+νµ = φ c νe+νe = 8.9 φ ντ +ντ (C.6) φ π νµ+νµ = 50.0 N ν exp −E/9.63 GeV . (C.9) 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 ν ) (C.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. (C.10) come from the dominant branching fractions to muons and electrons in kaon decays. To determine the asymmetry in muon neutrino and muon antineutrino events, we use the overall fluxes from pions and kaons in eqs. (C.8)-(C.9) and the charge ratios in eqs. (C.4)-(C.5).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.