Forward Neutrinos from Charm at Large Hadron Collider

The currently operating FASER experiment and the planned Forward Physics Facility (FPF) will detect a large number of neutrinos produced in proton-proton collisions at the LHC. In this work, we estimate neutrino fluxes at these detectors from charm meson decays, which will be particularly important for the $\nu_e$ and $\nu_\tau$ channels. We make prediction using both the next-to-leading order collinear factorization and the $k_T$-factorization approaches to model the production of charm quarks as well as different schemes to model their hadronization into charm hadrons. In particular, we emphasize that a sophisticated modeling of hadronization involving beam remnants is needed for predictions at FASER and FPF due to the sensitivity to the charm hadron production at low transverse momenta and very forward rapidity. As example, we use the string fragmentation approach implemented in \texttt{Pythia~8}. While both standard fragmentation functions and \texttt{Pythia~8} are able to describe LHCb data, we find that \texttt{Pythia~8} predicts significantly higher rate of high energy neutrinos, highlighting the importance of using the correct hadronization model when making predictions.

The currently operating FASER experiment and the planned Forward Physics Facility (FPF) will detect a large number of neutrinos produced in proton-proton collisions at the LHC.In this work, we estimate neutrino fluxes at these detectors from charm meson decays, which will be particularly important for the ν e and ν τ channels.We make prediction using both the next-to-leading order collinear factorization and the k T -factorization approaches to model the production of charm quarks as well as different schemes to model their hadronization into charm hadrons.In particular, we emphasize that a sophisticated modeling of hadronization involving beam remnants is needed for predictions at FASER and FPF due to the sensitivity to the charm hadron production at low transverse momenta and very forward rapidity.As example, we use the string fragmentation approach implemented in Pythia 8.While both standard fragmentation functions and Pythia 8 are able to describe LHCb data, we find that Pythia 8 predicts significantly higher rate of high energy neutrinos, highlighting the importance of using the correct hadronization model when making predictions.

I. INTRODUCTION
The forward production of charm quarks in high-energy proton-proton collisions at the LHC provides an excellent probe of the strong interactions.In the far forward region, corresponding to pseudorapidity η ≳ 7, which is beyond the coverage of the main LHC detectors, this process is sensitive to parton distribution functions at small momentum fractions x ∼ 10 −7 and at a scale Q ∼ m c .In this region of very small-x and small-Q 2 , that is not accessible to the direct measurement at HERA [1], deviations from the collinear factorization approach may be expected and novel small-x dynamics can occur, see e.g.Ref. [2].In particular, the non-linear contributions, which lead to saturation effects of the gluon density are expected to play an important role [3].In addition, the fragmentation functions needed to predict D-meson production are not well known in this regime, as they are usually constrained in e + e − collisions, while the hadronic environment of the proton-proton collisions introduces some new dynamical features which may lead to the factorization breaking.Future measurements of forward charm production will therefore provide a unique opportunity to study and test different aspects of QCD in this novel kinematic regime.
The study of forward charm production is also important in the context of large neutrino telescopes, such as IceCube.Here, charmed hadrons can be produced in cosmic ray collisions and their decay constitutes the source of prompt atmospheric neutrinos and hence a background to the extra-galactic neutrino signal [4][5][6].There are currently large uncertainties on the associated production rate and flux, underpinned by the lack of input data both at colliders as well as at neutrino telescopes.Indeed, IceCube has not seen evidence of prompt atmospheric neutrinos and only sets an upper limit on the corresponding flux [7].New input on forward charm production from the LHC will also improve the predictions for prompt atmospheric neutrino production.
The distribution of charmed hadrons at the LHC have been measured in the central region by ATLAS [8], CMS [9, 10] and ALICE [11,12] and in the forward region by LHCb [13][14][15].Together, these measurements cover the pseudorapidities |η| < 4.5, while at higher values charm production remains as yet unconstrained.This situation will soon change due to a new set of far forward experiments which will be able to detect and study neutrinos produced at the LHC.Many of these LHC neutrinos originate from the decay of charmed hadrons, and hence a measurement of the neutrino spectrum allows us to indirectly constrain forward charm production.The first two experiments, FASERν [16,17] covering η > 8.9 and SND@LHC [18,19] covering 7.2 < η < 8.7, have started their operation with the beginning of LHC Run 3 in summer 2022.Together, they will detect about ten-thousand neutrino interactions.Larger detectors with the ability to detect about a million neutrino interactions have been proposed in the context of the Forward Physics Facility (FPF) [20][21][22], which would operate during the HL-LHC era.
In anticipation of first data from the LHC neutrino experiments, it is important to have a dependable modeling of forward charm production and reliable predictions for the resulting neutrino fluxes and their uncertainties.Such estimates of the neutrino flux are needed as input for a variety of planned measurements, for example that of the neutrino interaction cross section.In addition, a comparison of theoretical predictions with the neutrino flux measurements will then allow to constrain QCD parameters, such as the mass of the charm quark, factorization and renormalization scales, parton distributions at small-x, and the fragmentation of the charm into D-mesons.In this paper we address these questions and study charm production at the LHC using two different QCD approaches: the perturbative collinear approach at next-to-leading order (NLO) and the k Tfactorization approach.In particular we constrain our models using available data from LHCb and make predictions for the LHC neutrino experiments.
Our main focus is on investigating the sensitivity of our calculations to the different modeling of the fragmentation of charm quarks into hadrons.This is especially important since the charmed hadrons are produced at very forward rapidity and at low transverse momenta.In this region additional effects may occur due to the interactions with beam remnants, and thus the standard fragmentation approach which is suitable for high transverse momenta may not be an applicable description in this kinematics.We base our analysis of different fragmentation schemes on the two different QCD models of charm pair production mentioned above to ascertain where the major source(s) of uncertainties and model dependence lie.In particular, our detailed analysis using the Pythia Monte Carlo event generator to model the fragmentation with color reconnection shows significant differences in the forward rapidity region with respect to the calculations using different fragmentation functions from the literature.This demonstrates that the forward particle production and the resulting high energy neutrino flux is particularly sensitive to the physics of fragmentation.
The paper is organized as follows.In Sec.II we review the experimental setup of the LHC neutrino experiments.We then discuss the modeling of charm production via the perturbative NLO calculation and the k T -factorization approach in Sec.III.Different approaches to modeling of the fragmentation, including standard fragmentation function approach and Pythia, are discussed in Sec.IV.The main results are presented in Sec.V. Finally in Sec.VI we present our summary and conclusions.

II. FORWARD NEUTRINO EXPERIMENTS AT THE LHC
The production of flavored hadrons has been extensively studied at all four main LHC experiments.This data provides a crucial input, for example for the modeling of high energy cosmic ray collisions and atmospheric neutrino fluxes.However, the most energetic hadrons are typically produced in the far forward direction, which lies outside of the coverage of the main LHC detectors.These particles are particularly relevant for modeling of cosmic ray collisions, since they carry a large fraction of the air showers energy and are also the source of the most energetic atmospheric neutrinos.While there are some measurements on far forward hadron production using additional LHC detectors, for example on photons and neutrons from LHCf [23,24], no data exists so far on strange and charm hadrons.Such input would, however, be desirable to address the cosmic ray muon puzzle [25,26] as well as to improve predictions for prompt atmospheric neutrino flux at neutrino telescopes [4][5][6].This situation is changing with the start of the LHC neutrino experiments, which will provide novel constraints on the far forward production of flavored hadrons.
Already in the 1980s it was noticed that the LHC would produce a large number of neutrinos through the decay of hadrons [27].Indeed, at each collision point the LHC generates an intense and strongly collimated beam of high-energy neutrinos along the beam collision axis.About 480 m downstream from the ATLAS interaction point this neutrino beam passes through the TI12 and TI18 tunnels, which housed the injector during the LEP era but remained empty during the LHC era.These locations provide unique opportunities to access the neutrino beam and study its properties.The first measurement illustrating the potential was performed by the FASER collaboration, which reported the first neutrino interaction candidates using a small pilot detector in 2021 [28].Following this proof of feasibility, two dedicated detectors have been installed in these locations.Located in TI12 is the FASER experiment [29][30][31].While it is mainly designed to search for light long-lived particles predicted by models of new physics [32][33][34][35][36], it also contains a dedicated emulsion neutrino detector called FASERν [16,17].This detector is centered on the beam collision axis and covers the pseudorapidity range η > 8.9.Located in TI18 on the opposite site of ATLAS is SND@LHC [18,19], which also contains an emulsion target as well as additional electronic components.Unlike FASER, it is positioned slightly off-axis and covers 7.2 < η < 8.7.Both detectors have the ability to distinguish neutrinos of different flavors and measure their energies.With the start of Run 3 of the LHC in summer 2022, both experiments have now started their operation and recently reported the observation of the first collider neutrinos [37,38].During LHC Run3, which is expected to last until 2025, the two experiments are expected to detect about ten thousand neutrino interactions with TeV scale energies.
Upgraded detectors to continue the LHC neutrino program are envisioned for the HL-LHC era.These would be located in the FPF [20][21][22], which is a dedicated cavern to be constructed 620 m downstream of ATLAS with the space to house a suite of experiments.In particular, three dedicated neutrino detectors have been proposed in this context: the emulsion based neutrino detector FASERν2, the electronic neutrino detector AdvSND, and the liquid noble gas based neutrino detector FLArE.Due to a tenfold increase in both target mass and luminosity these detectors have the potential to see more than a million neutrino interactions and study their properties in greater detail.
A first estimate of the neutrino flux has been provided in Ref. [39], taking into account both the prompt flux component from charm decays occurring the interaction point as well as a displaced component from the decay of long-lived light hadrons occurring further downstream from the interaction point.It uses a variety of different Monte Carlo event generators to model the production of hadrons at the LHC and employs a dedicated fast simulation to model the propagation and decay of long-lived hadrons when passing through the LHC beam pipe and magnetic fields.The results show that a majority of muon neutrinos and electron neutrinos at low energy originate from the displaced decay of light hadrons, while high energy electron neutrinos and tau neutrinos are mainly produced in the prompt decay of charmed hadrons.It was also noted that (i) there are large differences between the Monte Carlo generator's predictions for the prompt neutrino flux component of about an order of magnitude at high energies, and (ii) most of the generators have not yet been tuned or validated for charm production.More reliable predictions for forward charm production are needed.
Unlike for light mesons, the forward production of charm quarks can be described using perturbative QCD.This provides a different approach to obtain predictions for the LHC neutrino flux, which also offers a deeper connection to the underlying theory of QCD.Several recent studies have presented perturbative calculations for forward charm production at the LHC and derived corresponding predictions for the associated neutrino fluxes.In Ref. [40], the authors employed the collinear factorization approach at NLO, supplemented by additional k T -smearing and fragmentation functions to account for hadronization effects.Subsequent work by the same authors explored the associated PDF uncertainties [41] and the connection to astroparticle physics [42], also see Refs.[43][44][45].In Ref. [46,47], the authors used the k T -factorization approach, both in the full and hybrid realization, with fragmentation functions and a recombination model for hadronization.They also investigated the impact of an additional intrinsic charm component on the forward neutrino flux.
In the present analysis, we consider both of these perturbative QCD approaches and we particularly focus on the modeling of fragmentation.We present our predictions for the charm production at LHCb and forward neutrino fluxes at FASER.In the following sections, we provide detailed descriptions of the forward charm production modeling employing both QCD approaches.

III. CHARM QUARK PRODUCTION IN HADRONIC COLLISIONS
The standard routine for calculating the charmed hadron production cross section σ H is to fold the hadronic charm quark cross section σ c with a fragmentation function 1. Left: gluon-gluon fusion process for charm production in hadron-hadron collisions in the collinear factorization approach.f 1 , f 2 are the integrated gluon distribution functions which depend on the longitudinal momentum fractions x 1 , x 2 and the hard scale of the partonic sub-process.Right: the same process, illustrated for the case of forward production in the k T -factorization.The gluon x 1 is treated on-shell, and the gluon x 2 is off-shell with transverse momentum k T .σ is the partonic cross section which is on-shell (left panel) and takes into account off-shellness of one gluon (right panel).
In this section, we first focus on the perturbative calculation of charm quark production in hadronic collisions.In particular, we will describe and utilize two QCD approaches to calculate charm quark production: the NLO collinear factorization formalism and the k T -factorization formalism.A detailed discussion of fragmentation into hadrons will then be provided in Sec.IV.
The production of charm in hadronic collisions is dominated by the gluon-gluon scattering.In this process, gluons from two colliding hadrons fuse and produce a charm quark-antiquark pair which subsequently fragments into the hadrons.The generic diagram for gluon-gluon fusion process in hadronic collision is illustrated in the left panel Fig. 1, where the cross section can be factorized into two gluon distribution functions, f 1 and f 2 , and the perturbatively calculable partonic cross section σ.This is the framework at the root of collinear factorization approach.
In the forward region, this process probes the kinematics where the two incoming partons have very different longitudinal momenta.The longitudinal momentum of the forward charm quark at high energy is approximately equal to x F ≃ E c /E p where E p is the energy of the incident proton and E c is the energy of the charm quark.Since we are interested in TeV-energy neutrinos from TeV-energy charm decay, the corresponding forward charm production kinematics probes values of x F of order 0.1 or higher.This in turn means that the longitudinal momentum fraction of one of the gluons is large, x 1 ∼ x F , and the other one is very small.To be precise the longitudinal momentum fraction x 2 ≃ m 2 cc /(x F s), where m cc is the invariant mass of the produced charm-quark pair and √ s is the center of mass energy of the hadronic collision.This means that for high energies the forward production is particularly sensitive to the gluon density at very low values of x 2 ≳ 10 −7 , which is not constrained very well in this region.Thus the forward production offers unique possibilities for tests of novel QCD dynamics in the region of small-x.

A. Collinear Factorization at NLO
The double differential NLO cross-section for charm pair production is given by the expression where m c is the charm mass, √ ŝ = √ x 1 x 2 s is the partonic CM energy, µ F and µ R are the factorization and renormalization scales respectively, and f i,j represent the quark and gluon parton distribution functions (PDFs) as appropriate.As noted previously, we compute the cross-section to the next-to-leading order in perturbation theory.The double differential cross-sections for charm quark production are calculated using the FONLL code [48,49], which provides an interface to LHAPDF [50,51], thus allowing one to use a variety of up-to-date PDFs.We choose to use the central CT14nlo PDF set [52] from the LHAPDF database as a representative set for our analysis.While there are more recent PDF sets available in the literature, including those that have been fit to 13 TeV LHCb data and consequently have reduced uncertainties in their predictions at low-x [53,54], we find that uncertainties in the cross-section from scale variation dwarf those from using different PDFs.Instead, our choice of the central CT14nlo PDF allows us to maintain compatibility with results obtained in Ref. [4], while also using m c = 1.3 GeV, consistent with the PDG best-fit.
To obtain best-fits to current charm data, we choose to vary the factorization scale µ F and renormalization scale µ R while keeping the charm mass fixed.Assuming the scales vary proportionally to the charm transverse mass, , it has been the norm to vary these parameters independently within a range from (0.5-2.0) ∝ m T .However, when restricting ourselves to this narrow range, we find that at high energies √ s ⩾ 7 TeV fits to data become progressively worse with increasing rapidities.Furthermore, determining uncertainties around the best-fit scales also requires one to extend the search beyond this range.Therefore, we allow these parameters to vary independently over a broader range ∈ [0.5, 8.0] unencumbered by theoretical preferences, allowing the best-fit parameters to be instead determined by fitting to data.We also determine the parameters defining a 1σ uncertainty band around the best-fit cross-section.
We compute cross-sections for a range of parameters (µ R , µ F ) and obtain, for each choice of fragmentation scheme, the meson cross-section that may be fit to data from LHCb.The end result of this fitting exercise is that we obtain different sets of best-fit (µ R , µ F ) for different fragmentation scheme.We defer the details of our fitting procedure to Appendix A.

B. k T -factorization
In the forward regime, one should apply a framework which incorporates resummation of the large logarithms α s ln 1/x.This is accomplished through the k T -factorization formalism [55][56][57].The k T -factorization formalism involves off-shell matrix-elements for partonic scattering and unintegrated gluon distribution 1 functions F(x, k T ) which depend on the transverse momentum vector k T of the off-shell gluons.The unintegrated gluon distribution functions encode more detailed information about the dynamics of the partons, and can be especially important in providing information about the details of the kinematics of the event.The k T -factorization approach in hadroproduction of heavy quarks has been considered in Refs.[55][56][57] where the off-shell matrix element for heavy quark production have been derived.The expression for the cross section in the k T -factorization formalism has the following form, see e.g.[55] where the off-shell partonic cross section σ off contains contributions from gluon-gluon scattering, dominant for the high energy limit.For the specific case of forward charm production considered here, due to the fact the kinematics is very asymmetric and one gluon has large longitudinal momentum fraction x 1 it is appropriate to use an approach in which this gluon is treated on-shell and satisfies the DGLAP evolution.Therefore the formula Eq. ( 3) in this limit becomes where σ on−off can be obtained from σ off by setting one gluon on-shell, see Ref. [55].This is illustrated in the right panel of Fig. 1.The gluon with large longitudinal momentum fraction x 1 is indicated together with schematically drawn collinear cascade originating from one proton.On the other hand, the gluon with very small x 2 has transverse momentum k T and it is produced as a result of a very long cascade of emissions from the other proton.These emissions are not collinear, hence their transverse momenta are not ordered.Therefore such cascade leads to the diffusion in the transverse momentum distribution.This approach was used in Ref. [2] with the large x 1 gluon in the DGLAP collinear regime, which is on-shell and the small x 2 gluon off-shell, with appropriate approximation of the matrix element.
In this work we are interested in the differential distributions in rapidity, which can be obtained by generalizing collinear formula Eq. ( 2) to include the transverse momentum dependence.Since we are using expressions from [55], which are formally lowest order, the differential cross section can be taken as √ s exp(−y 4 ) as well as the transverse masses m 2 3,4T = p 2 3,4T + m 2 c of the quark and antiquark (see also [46]).The unintegrated gluon distribution functions within the high-energy formalism need to be computed from the appropriate evolution equations which incorporate the small-x dynamics.The unintegrated parton densities within the high energy formalism are usually computed from the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation which resums the powers of α s ln 1/x [61,62].It has been computed at leading logarithmic (LL) and next-to-leading logarithmic order (NLL) in QCD.For the phenomenological applications it needs to be supplemented by the additional corrections which take into account higher orders in the form of kinematical constraints and the constraints from matching to the DGLAP evolution [63].In addition, in the limit of high energies, or very small-x, other corrections are expected to occur, which are related to the parton saturation phenomenon [3].In this regime, the gluon densities are so large that recombination effects need to be taken into account which are expected to slow down the growth of the gluon densities.These corrections lead to the appearance of the non-linear terms in the small-x evolution equations.The non-linear evolution leads to the taming of the gluon distribution in the region of very small-x and moderate to small values of scales k T .To be specific, these evolution equations generate the x-dependent saturation scale Q 2 s (x).Whenever the relevant scale of the process, say the k T of the gluon, is smaller than Q 2 s (x) non-linear terms are very important, while for k 2 T > Q 2 s (x) they can be neglected and the non-linear evolution equations give results which coincide with thus obtained form the linear evolution.
The effective theory for high density at small-x is the Color Glass Condensate [64][65][66][67][68][69], with the corresponding JIMWLK evolution equations.In the multicolor limit the hierarchy of JIMWLK equations reduces to the Balitsky-Kovchegov equation [70,71], the latter being the BFKL equation supplemented with the nonlinear term in the gluon density.
The small-x unintegrated gluon density for the present paper was taken from Ref. [72] as well as from Ref. [73].The gluon in Ref. [72] which was based on the unified BFKL+DGLAP evolution supplemented with small-x resummation [74].Two sets of gluon distributions were used: based on linear evolution as well as non-linear evolution cast in the momentum space [75,76].The latter one includes the non-linear term in density which is responsible for the saturation effects.Both sets of distributions were fitted to the data on F 2 structure function at HERA.The non-linear term is important for low-x and low values of transverse momenta and leads to taming of the gluon distribution and therefore the resulting observable cross section.We also used the gluon extracted from more recent fit in Ref. [73] to HERA data, which was based on the full resummation [63,77] including the BFKL at NLO.

IV. CHARM FRAGMENTATION
In the previous section we have discussed the perturbative aspects of charm production.We now turn to question of fragmentation of charm quarks into charm hadrons, which is a non-perturbative process and requires a separate treatment.Here we first review the standard fragmentation function formalism as well as its short-comings.We then present an alternative approach based on the modeling of hadronization in Monte Carlo generators.

A. Fragmentation Functions
Many studies of charm production at the LHC make use of the factorization theorem to separate the charm production and fragmentation process.In the literature, the latter is then modeled via fragmentation functions that have been extracted from lepton collider data, assuming that they are also applicable at hadron colliders.As we will explain later, this may not be appropriate at hadron colliders, especially in forward and low transverse momentum region that is most relevant for FASER.In this approach, one uses the fact that charm quarks in electron-positron annihilation are produced with a known momentum, for example with p c = m Z /2 at LEP.One can then measure the flavor and momentum of charmed hadrons p H to constrain the fragmentation process.This is typically parameterized in terms of fragmentation fractions f H , describing the probability of a charm quark to form a specific charm hadron H, and a fragmentation function D H (z), describing the distribution of fractional energy inherited by the hadrons z = p H /p c .In a later comparison of fragmentation approaches, we use the fragmentation fractions f D + = 0.244, f D 0 = 0.606, f D + s = 0.081 and f Λ + c = 0.061 as obtained in Ref. [78] and the Peterson fragmentation function [79].It has the form where we choose ϵ = 0.035 following Ref.[80].Note that the same fragmentation function is used for all charmed hadrons.Simply for illustration, we will also consider the unphysical case with no fragmentation beyond fragmentation fractions.This means that quark and hadron momenta are identical, implying D H (z) = δ(z − 1).
Although the above-mentioned fragmentation functions approach has been successfully applied to measurements of charm production in the central and high-p T region of the LHC, it faces additional challenges in the forward and low-p T regime.There are a variety of hadron collision measurements that contradict the predictions obtained using fragmentation function; see Sec. 6.2.2 of Ref. [22] for a pedagogical overview.In the following, we summarize three important observations that are particularly relevant for the modeling of forward charm production: • The first observation concerns the production asymmetry of charmed mesons and their antiparticles.While the fragmentation function approach predicts equal production rates of charmed hadrons and their anti-particles, an excess of D − compared to D + has been observed at high x F in π − -nucleus fixed target collisions recorded by WA82 [81], E769 [82] and E791 [83].Such production asymmetries in the forward direction are typically explained by charm hadronization involving the beam remnants [84].In the case of π − -nucleus collisions, the c can hadronize with the valence d from the pion and form an energetic D − meson.In contrast the formation of a D + requires a c and d.Since the d cannot be a valence quark, but either a sea-quark or produced otherwise, the D + mesons are expected to be less energetic.This effectively induces a production asymmetry at high x F .
• The second observation regards the energy spectra.Using the same data from pion fixed target experiments, it has been found that the momentum spectrum charm of hadrons are about as hard as or even harder than the charm quark spectra obtain from perturbation theory [81][82][83].This contradicts the fragmentation functions approach, which predict the hadrons to be softer than the charm quarks.In contrast, the above-mentioned mechanism of hadronization with other light quarks in the event, especially valence quarks from the beam remnant, would naturally allow the hadrons to be more energetic than the charm quarks and explain this observation.
• The third observation relates to the baryon to meson production ratios.Recently, ALICE has measured the ratio between the Λ c baryon and D 0 meson production rates in the central region and found that this ratio increases from about 10% at high transverse momentum to about 50% at low transverse momentum [85][86][87].A similar enhancement was also seen by CMS [88].This disagrees with the expectation from fragmentation functions applied in the lab frame and extracted from LEP, which predict a roughly constant Λ c to D 0 ratio of around 10%.
The observations above illustrate that fragmentation functions extracted from lepton colliders are not sufficient to describe charm production at hadron colliders.

B. Hadronization using MC Generators
One way to address the abovementioned problems is to use more sophisticated models of fragmentation which are typically implemented in Monte Carlo generators.Here, we will take advantage of these efforts and use Pythia 8 [89,90] to model hadronization.Pythia uses the Lund string model [91,92] in which colored objects are connected by a color string containing the field lines of the strong force.This model can intuitively explain two of the above observations: a charm quark connected to a beam remnant valence quark will be pulled forward, and hence gain energy, or even hadronize together with the valence quark, leading to a production asymmetry.By default, Pythia uses the Monash tune [93].While broadly used to describe phenomena at the LHC, we note that it is not able to properly describe the baryon enhancement observed at ALICE.This problem is addressed by a newer QCD-inspired color reconnection scheme introduced in Ref. [94].
It allows for different string topologies, such as junctions of three strings, which leads to a higher baryon production rates in high-multiplicity regions.It has been also recently suggested [95], using modeling with Pythia, that this QCD-inspired color reconnection mechanism might be essential for the proper description of the J/ψ production at the LHC.Throughout this work, we use the mode 2 configuration introduced in Ref. [94].
One practical complication is that the tools we use to model the perturbative production of charm quarks do not generate events that can be used as input to Pythia, but only provide the charm quark distribution d 2 σ c /(dp T,c dy c ).We bypass this problem by using a re-weighting approach which is inspired by Refs.[96,97].To understand the underlying idea, let us recall that, conceptually, we can write the charm hadron distribution d 2 σ H /(dp T,H dy H ) as a convolution of the charm quark distribution d 2 σ c /(dp T,c dy c ) and a (unitary) transfer function f (⃗ p c , ⃗ p H ) describing the hadronization process: In general, the transfer function would depend on both the quark and hadron momenta as well as the collider setup.In the fragmentation function approach, we assumed that and that it is independent of the collider setup.In Monte Carlo generators, the hadronization procedure is more complex and f cannot be parameterized by a simple function.However, the transfer function is encoded in a generated event output: the charm production process of Pythia provides a sample of events, where each event is characterized by the parton momentum ⃗ p c , the hadron momentum ⃗ p H , a hadron ID and an event weight w.The events in the sample implicitly follow a distribution d 2 σ P 8 c /(dp T,c dy c ) for the charm quarks and d 2 σ P 8 H /(dp T,H dy H ) for the charm hadrons related via a Eq. ( 6) through a transfer function f .
To apply the same hadronization to a different model of charm production, we use the reweighting procedure and adjust the weights By construction, the events will then follow a d 2 σ c /(dp T,c dy c ) at quark level.The hadrons follow the desired distribution which we can extract from the event sample.Let us summarize our approach.The usual fragmentation function approach assumes that the charm hadronization process is described by a transfer function of the specific form f = D H (z), which is universal for all colliders, applicable to all predictions of charm quark production, and independent of the charm quark kinematics and hadronic environment.We saw, however, that this assumption is invalid at hadron colliders.For example, hadronization with beam remnants, that is not captured in the fragmentation functions, leads to a harder forward charm hadron energy spectra and a charge asymmetry.This has been observed at past beam dump experiments and is expected to be important for forward charm hadron production at the LHC.
We therefore propose an alternative approach to model charm hadronization using Pythia, which only assumes that the underlying transfer function f is the same for different predictions of charm quark production, and that Pythia provides a reasonably good prediction of hadronization especially in the forward direction.We note that the accuracy of Pythia's description of forward charm hadronization, especially with beam remnants, has not yet been experimentally tested the hadronization process due to a lack of experimental data.However, Pythia's good description of charm hadrons at beam dumps as well as light hadrons in the forward direction of the LHC [23] provides some confidence in its overall description of hadronization.

V. RESULTS
In this section, we present and discuss the results of our different charm production models.We start this by systematically varying the modeling.For each considered setup, we shall show comparisons of our predictions to the double differential cross section of D 0 meson measured at 13 TeV by LHCb as well as the expected neutrino event rates at FASERν.To determine the neutrino flux, we follow the same approach as Ref. [39].Initially, the charm hadrons are decayed in their rest frame according to the decay branching fractions and energy distributions obtained with Pythia.Subsequently, the resulting neutrinos are boosted into the laboratory frame and recorded if they pass through the detector's cross-sectional area.To obtain the anticipated number of neutrino interactions in the target volume, we convolute the neutrino flux with the interaction cross-sections obtained by GENIE [98].Here, we consider FASERν to consist of a 25 cm × 25 cm × 1 m tungsten target [16].
In the following, we will present results for collinear factorization in Sec.V A and k T -factorization in Sec.V B. We will compare both approaches and show additional distributions in Sec.V C.

A. Collinear Factorization at NLO
We first consider the calculation using the NLO collinear factorization.As described in Sec.III A, we obtain multiple best-fit cross-sections corresponding to different fragmentation schemes.We find that a variation of the scale parameters (µ F , µ R ) mainly influences the normalization of the crosssection predictions, while the shape of the p T distribution remains largely unchanged.In contrast, the latter is more significantly affected by the choice of the fragmentation scheme.
We show a comparison of these results to the LHCb data in the left panel of Fig. 2 for three different modeling approaches for fragmentation.The green dotted line shows the best fit prediction obtained using a constant fragmentation factor.The best-fit cross-section in this case is obtained for (µ F , µ R ) = (2.1, 1.6) m T consistent with results from [4].However, we find that the p T shapes of the corresponding double differential cross-sections are inconsistent with LHCb data, consistently overestimating at high p T .With change of scales primarily affecting cross-section normalizations, and not the shape, there is no way to improve the fit within the realm of our analysis when using constant factors for fragmentation.Thus, this demonstrates the importance of including more realistic fragmentation schemes.The blue dashed lines show the best fit results using the Peterson fragmentation function, obtained for (µ F , µ R ) = (3.75,1.75) m T .These agree reasonably well with LHCb data for all rapidity regions, while still overestimating the data at low p T somewhat.Finally, best-fit results obtained using Pythia for fragmentation are shown as red solid lines.These correspond to (µ F , µ R ) = (2.25,1.5) m T .We observe that this setup produces similar results to those using the Peterson fragmentation function in the regime accessible to LHCb, with slight differences mainly at low p T < 2 GeV.
We proceed to evaluate the electron neutrino flux from charm hadron decay at FASERν from these simulations.The results are shown in the right panel of Fig. 2. With Peterson's fragmentation, the obtained flux has lower rates and peaks at lower energies compared to the scenario without any In the left panel, we compare these predictions with measurements of the double differential neutral D-meson production rate obtained by LHCb at 13 TeV.We present results for three different rapidity regions, where the results at higher rapidity were scaled.In the right panel, we show the resulting number of electron neutrinos from charm hadrons decay that interact with the FASERν detector as a function of the neutrino energy.For context, we also display in black the event rate resulting from neutrinos from light hadron decays as obtained in Ref. [39].See the main text for a detailed discussion.fragmentation.This outcome is expected since in the fragmentation function approach, the charm hadron is always less energetic than the charm quark.In contrast, using Pythia for fragmentation increases the neutrino flux and shifts it to higher energies compared to the scenario without any fragmentation.As discussed in Sec.IV, this outcome is consistent with observations at beam dump experiments, where hadronization with beam remnants plays a role.We emphasize that despite both fragmentation choices providing similarly good descriptions of the LHCb data, they lead to a significant difference in neutrino event rates at FASERν, differing by about one order of magnitude.This highlights the importance of properly modeling fragmentation for forward charm and, consequently, neutrino flux predictions for FASERν and other LHC neutrino experiments.For comparison, we also show the event rate from light hadron decays in black, as obtained in Ref. [39], using various generators.The solid line represents the central prediction, while the shaded band shows the range of predictions from different generators.This line is meant to provide optical guidance and to illustrate regions where light and charm hadron decay contributions dominate the electron neutrino flux.While our prediction already agrees reasonably with the LHCb data, we observe an underestimation of events at intermediate p T ∼ 8 GeV and a mild overestimation at p T ∼ 1 GeV when compared to experimental measurements.As pointed out in Ref. [40], including an additional k T smearing, which aims to capture both an intrinsic transverse momentum of the initial state partons as well as some soft gluon emission effects, can help improve the agreement with data.The authors achieve this by introducing a Gaussian smearing with width ⟨k T ⟩ of the transverse momentum of the charm, while keeping its rapidity constant.However, we note that this approach does not conserve energy and can lead to charm quarks that are more energetic than the proton beam.Indeed, this leads to an unphysical order of magnitude increase of the neutrino event rate at high energies.To address this issue, we modify the smearing such that the z component of charm quark momentum is kept constant and the rapidity is allowed to change.
By iterating over a range of values of ⟨k T ⟩ (see Appendix B), we find that the best agreement to data is for a combination of ⟨k T ⟩ = 1.5 GeV and (µ F , µ R ) = (1.75, 1.25) m T .This value of ⟨k T ⟩ is consistent with Ref. [42], which uses ⟨k T ⟩ = 1.2 GeV, as well as with the default transverse momentum for hard interactions used within Pythia, which is 1.8 GeV.The corresponding neutrino fluxes are not highly sensitive to the choice of ⟨k T ⟩.
In order to illustrate why different fragmentation approaches give similar results for the LHCb data, but very different neutrino flux in the forward region, we show p T distribution for large y > 6 and the rapidity distributions for low p T in Fig. 3.We note that for lowest value of p T and large rapidity, calculations with various fragmentation schemes differ significantly.
Up to now, we've only shown our central prediction, which uses scale choices (µ F , µ R ) = (1.75, 1.25) m T that were obtained by fitting the data with ⟨k T ⟩ = 1.5 GeV.The same fit also allows to define scale uncertainties in a data driven way (see Appendix A for details).To illustrate this, we present in Fig. 4 our results for two additional scale choices, which provide an error band that encompasses the LHCb data.Looking at the right panel, the corresponding neutrino fluxes show only mild sensitivity to the choice of scales.

B. k T -Factorization
We have observed that introducing an additional k T smearing improves the agreement of the collinear factorization prediction with data.This smearing effectively simulates intrinsic transverse momentum and soft-gluon emissions in the initial state.These effects are naturally included in the k T -factorization approach due to the presence of the unintegrated gluon distribution function and the off-shell matrix element which depend on transverse momentum k T .As discussed in Sec.III B, we are using a hybrid approach which utilizes an unintegrated PDF for the low-x gluon and an integrated PDF for the high-x gluon.This is because ultimately we are interested in the very forward region where one x is very small and the other very large.
As the basic setup we choose the unintegrated gluon distribution from the Kutak-Sapeta (KS) calculation [72] using the nonlinear evolution, and for the large-x we use the CT14nlo gluon.Since the KS gluon has been fitted to the HERA data using the leading order strong coupling constant, we use the same setup for the one power of strong coupling in the formula for the cross section.The second power of the coupling is taken at NLO consistent with the CT14nlo PDF used for large-x gluon.As before we are are modeling the hadronization using Pythia with the QCD-inspired color reconnection scheme.The results are shown in Fig. 5 by the blue curve.We observe that the calculation has the right shape in p T but it significantly underestimates the experimental data.This was also observed in calculation of [46].This is not totally unexpected since the off-shell partonic cross section used in k T -factorization is effectively computed at the LO [55].Therefore when compared with NLO collinear calculation it does not have virtual terms as well as final state gluon emissions from the quarks.It also has an off-shell gluon only on the small-x side.Given that the NLO calculation in the collinear approach resulted in K-factor of the order of 2.5 with respect to the LO result, see e.g.[99], it is expected that the k T -factorization will likely have large K-factor as well.
In order to get the normalization to agree with LHCb data, and therefore make our extrapolations from LHCb to FASERν more reliable, we introduce a normalization factor which we refer to as k-factor2 in this calculation, determined by a fit to the data (with additional weights that ensure each rapidity bin contributes identically to the χ 2 measure).We find a best fit of k = 2.32; the resulting double-differential cross-section is illustrated by red line in Fig. 5.This is in excellent agreement with the LHCb data over the full p T and rapidity range.This is encouraging since it means that the x dependence of the unintegrated gluon, correctly reproduces the rapidity dependence, and also the p T dependence is correctly captured.We shall also see, that the k-factor does not change between the 7 and 13 TeV.We also determine an uncertainty of the fit, as illustrated by the orange and magenta curves in the same figure (using a rescaled χ 2 for this following the PDG procedure described in Refs.[80,100]).These variations form a nice envelope around the data with a width of about a factor 2 at low values of p T .The right panel in Fig. 5 shows the electron neutrino flux obtained in this approach.A similar size band is also obtained at FASERν, see right panel.
Next, we study the dependence of the results on the choice of the low-x unintegrated gluon distribution.In Fig. 6 we show our results for three choices of unintegrated PDFs: two choices for the KS gluon with linear evolution and with non-linear effects that describe saturation effects, and third choice of gluon from [73] obtained from the linear evolution including the resummation using the Ciafaloni-Colferai-Salam-Stasto (CCSS) approach [63].We find that the prediction which includes saturation effects is in excellent agreement with the LHCb data over the full p T range.In contrast, the linear cases overshoot the data at low p T .However, given that the results include the k-factor effectively added by fitting as explained before, it is not possible to conclude at this moment about the importance of the saturation effects in the LHCb data.Looking at the right panel, including saturation effects results in a reduction of the flux by a factor of approximately three compared to the linear case.This is due to the fact that the nonlinear effects are largest at very low p T .
We have also tested the sensitivity of the k T factorization calculations to the choices of the large x gluon distribution, the running coupling order and the scale choice.The results of these studies are collected in Appendix C. We have found rather small differences between the calculations for these various choices.

C. Comparison of Approaches
Based on the previous discussion, we identify central predictions for both factorization approaches.In particular, we consider the following configuration • collinear factorization at NLO with CT14nlo for the gluon parton distribution, renormalization scale µ R = 1.75 m T , factorization scale µ F = 1.25 m T , a k T smearing with ⟨k T ⟩ = 1.5 GeV, and fragmentation modeled with Pythia with the QCD-inspired color reconnection scheme • k T -factorization using KS unintegrated distribution for the low-x gluon including saturation effects, the CT14nlo parton distribution for the high-x gluon, a k-factor of 2.32 and fragmentation modeled with Pythia with the QCD-inspired color reconnection scheme In Fig. 7, we compare the corresponding distributions from both approaches.We note that k Tfactorization with saturation gives slightly better description of the p T shape of the LHCb data than the NLO case.However, we again remind the reader, that this has to be taken with caution since this calculation includes the fitted k-factor which is not needed for the NLO collinear approach.We find that both approaches give good description of D 0 + D0 data but when compared with LHCb data for D + + D − , and for D s + D − s , the low p T region is overestimated.We show distributions as a function of rapidity for different p T regions, and we find that NLO and k T -factorization with saturation give similar values for central rapidity, but they differ at large rapidity, by about a factor of 2, especially for 0 < p T < 0.5 region.For large values of p T , this difference is reduced.
In Fig. 8, we also show comparison of both approaches with the LHCb data at 7 TeV, and the description is very good in both cases.It should be stressed that the used scales for the NLO calculation and the k-factor for k T -factorization at 7 TeV are the same as extracted from 13 TeV data.As mentioned previously, this is encouraging since it means that the energy dependence of the data, which is driven mainly by the x evolution of the unintegrated gluon density is captured correctly.The latter one has been taken from the resummed approaches [63,74,77] which aim to reproduce both small-x and collinear dynamics.
The neutrino flux obtained using both QCD approaches is presented in Fig. 9.The upper row shows the number of interacting neutrinos in FASERν operating during LHC Run 3 with an integrated luminosity of 150 fb −1 while the bottom row shows the neutrino events rate at FLARE at the FPF during the HL-LHC with a luminosity of 3 ab −1 .The three columns correspond to the three neutrino flavors.The shape of the neutrino flux remains similar for all neutrino flavors in both approaches, with the NLO contribution slightly lower than that of the k T -factorization.However, the two approaches are very close and fall within the range of uncertainty, which is approximately a factor of two.The black lines represent the contribution to the neutrino flux from decays of light  hadrons.Notably, we find that the dominant contribution to neutrinos occurs above 500 GeV for ν e and above 1 TeV for ν µ .Detecting ν τ would serve as a direct test of charm production, as there is no contribution from pions and kaons decays.
Based on our calculation, we predict that FASERν during LHC Run 3 is expected to observe approximately 4000 ν e , 4000 ν µ , and 120 ν τ charge current interactions originating from decays of charm hadrons.The FPF, proposed to house larger neutrino detectors during the HL-LHC era, aims to record a significantly larger sample of neutrino interaction events [21,22].Specifically, we consider the FLARE detector housed within the FPF, for which we assume a 1 m×1 m×7 m liquid argon target [22].We can see that FLARE will detect approximately 1.4 × 10 5 ν e , 1.4 × 10 5 ν µ , and 6000 ν τ from charm hadron decays.This substantial increase in statistics will enable FPF experiments to conduct more detailed tests on forward charm production and provide the necessary data to distinguish between different predictions.

VI. CONCLUSION
Forward charm production at hadron colliders has long been recognized as an sensitive tool for probing the strong interaction.However, until recently, it has remained beyond the reach of the existing LHC experiments.This situation is now changing with the start of operation of the FASERν and SND@LHC experiments, which are strategically positioned in the far-forward direction of the LHC and specifically designed to detect collider neutrinos.Many of these neutrinos originate from the decay of charm hadrons, presenting a unique opportunity to investigate forward charm production.Together, FASERν and SND@LHC are projected to observe approximately ten thousand neutrinos during the LHC's Run 3, spanning from 2022 to 2025.Looking forward, a continuation of this collider neutrino program is envisioned for the HL-LHC era from 2029 to 2042: by utilizing larger detectors situated in the FPF it will be possible to detect millions of collider neutrinos.
In this work, we have predicted neutrino fluxes from charmed mesons in these forward neutrino experiments.To this end, we have modeled charm hadron production from pp collisions at 13 TeV using different QCD and hadronization models, fitting our hadron cross-sections to charmed meson data from LHCb to ascertain the values of parameters involved in our models.This also allows us to determine which QCD and hadronization models are well tailored to describing physics at the forward rapidities that will be probed at FASERν.When evaluating hadron cross-sections against current collider data, we have placed particular emphasis on the hadronization models used to convert charm cross-sections to hadronic ones.We have discussed how current fragmentation function based models in the literature are not especially well motivated to describing far forward physics, because, among other things, they omit the potential for involving beam remnants when hadronizing.With the end goal of accurately forecasting neutrino fluxes at FASERν, we have, instead, devised a scheme that employs the string fragmentation model implemented in Pythia 8, resulting in a more realistic representation of hadronization.This Pythia-based scheme naturally overcomes most of the theoretical shortcomings of fragmentation function based models.We also demonstrate that the use of this hadronization scheme leads to a significantly enhanced flux of forward neutrinos compared to those obtained using established fragmentation functions, which results from allowing the hadronization with beam remnants.This underscores the importance of utilizing an accurate fragmentation modeling.However, we also note that the topic of forward charm hadronization warrants further theoretical investigation.
To obtain the charm cross-sections that underpin our analysis, we have investigated two distinct QCD models and made noteworthy improvements to each insofar as they apply to forward physics: a) collinear factorization, where we use factorization and renormalization scales as free parameters to be determined by fitting to LHCb data, as is typically done in the literature, but in addition apply a k T smearing on the charm transverse momentum in an energy conserving way; and b) k T factorization, which is more suitable for the description of the forward particle production at high energy since it resums contributions due to the small x effects in the parton density, and where we include a k-factor to account for a mismatch in the normalization against LHCb data.When using the former, we find that -no matter the variation of scales -the agreement of the shape of final hadron differential cross-sections vis-à-vis LHCb data is noticeably improved by the allowing a Gaussian smearing of the charm transverse momentum with some mean k T .In contrast with Ref. [40], where this k T smearing effect has been first discussed, our analysis explicitly conserves energy when applying the transformation by keeping the charm z-momentum constant and allowing its rapidity to vary.We find a best-fit to 13 TeV LHCb data is obtained for {µ F , µ R } = {1.75,1.25}m T alongwith ⟨k T ⟩ = 1.5 GeV.When using the k T factorization scheme, our central prediction incorporates the KS unintegrated distribution with a non-linear evolution for the low-x gluon, and the CT14nlo distribution for large-x gluons.A salient feature of our analysis is that, in order to describe the LHCb data, we introduce a constant k-factor of 2.32, determined by fitting the overall normalization to data.The need for the inclusion of the normalization k-factor in k T -factorization approach likely stems from the fact that the off-shell partonic cross section for production of heavy quarks is only available at lowest order.To theoretically ascertain the value of the k-factor, higher orders of the off-shell partonic cross section will need to be computed, and possibly resummed.We further examined the impact of systematically varying the underlying QCD parameters, such as scale selection and parton distribution function choices at low and high x on these predictions.
We note that both QCD approaches provide a good description of the LHCb data when paired with the PYTHIA hadronization scheme.They exhibit similar energy dependence in the neutrino flux with slightly different overall normalizations, which can be attributed to the specific QCD parameter choices.More theoretical work with respect to the underlying uncertainties relevant to each model is needed to improve the precision of these calculations (in particular the small x approach), as well as further experimental input in order to distinguish between NLO collinear and k T -factorization, and especially to see an onset of parton saturation.
Using our best-fit QCD models, we have shown predictions for neutrino fluxes of all three flavors for the ongoing FASER experiment as well as the proposed FLARE detector at the FPF and compared them against those from the decay of lighter mesons.We find that, depending on the choice of QCD scheme, the electron neutrino flux from charmed mesons dominates over those from pions and kaons starting at neutrino energies between 400 and 500 GeV.Furthermore, muon neutrinos from charmed meson decay become comparable to those from pion and kaon decays at energies above 1 TeV for both QCD approaches.Tau neutrinos, produced exclusively from heavy meson decays, provide a background-free channel to investigate heavy meson QCD.Our models predict between 4000 and 6000 charged current tau neutrino interactions at FLARE with energies around 1 TeV, depending on whether one uses the collinear NLO scheme or the k T -factorization scheme respectively.
The first observation of collider neutrinos at FASER [37] heralds the opening of a new frontier towards significantly improved understanding of forward QCD.Further measurements at both FASER and SND@LHC will provide a unique opportunity to gather valuable information about small-x QCD, validity of k T -factorization and NLO collinear approach, and validity of different hadronic fragmentation scenarios at forward rapidities.In the future, the planned experiments at the FPF, with significantly improved statistics, will become the ideal place to unravel these most important facets of QCD.In addition, we expect that measurements of the forward neutrino production at the LHC will provide valuable inputs for estimating the prompt neutrino flux, reducing its theoretical uncertainties and thus providing a better understanding of the main background for the detection of ultra-high energy neutrinos be it from extragalactic astrophysical sources or from beyond standard model physics including dark matter decays and annihilation.cross-sections observed at LHCb [14].These include cross-sections for D 0 , D ± , and D s at rapidities between 2 ⩽ y ⩽ 4.5 binned by 0.5, i.e. five bins in y for each meson.We focus on this subset of collider data to determine the scales, µ R and µ F , that best describes it.
To compare our theoretical cross-sections against charmed meson cross-sections from LHCb, we first compute the double differential cross-section d 2 σ cc /dydp T for bare cc pair production using specific values for the fragmentation and renormalization scales.We then assume a specific fragmentation scheme, without any additional free parameters, to hadronize the charm quarks into hadrons.The resulting differential cross-section distribution at this stage may now be compared against corresponding LHCb data for a measure of its goodness of fit, which we achieve by means of a simple χ 2 analysis.Since accurately forecasting high rapidity cross-sections is critical towards obtaining predictions for forward experiments like FASER, it becomes important to ensure that the goodness of fit analysis is not skewed by the availability of significantly more data at LHCb's lower rapidities 2 ⩽ y ⩽ 3 rather at, say, y ⩾ 4. We therefore use a χ 2 measure that is normalized to the number of p T bins with cross-section measurements for each rapidity bin in the LHCb data, ensuring that each bin carries equal weight towards the measure.
Repeating this procedure for multiple (µ R , µ F ) values, we generate a range of cross-sections and, for a given fragmentation scheme, we ascertain the best-fit value of these parameters as the one that minimizes the χ 2 /d.o.f.Likewise, we also obtain the parameters corresponding to a 1σ region of variation around the best-fit cross-section.The gives us a set of best-fit (µ R , µ F ) parameters for each choice of fragmentation scheme.using different fixed values for ⟨k T ⟩ = 0, 0.5, 1.0, 1.5, 2.0 . .., we find that the best agreement with data is obtained for ⟨k T ⟩ = 1.5 GeV.As shown in the right panel of Fig. 10, the corresponding neutrino fluxes are not highly sensitive to the choice of ⟨k T ⟩.

FIG. 2 .
FIG.2.Modeling of Fragmentation: Predictions obtained using collinear factorization at NLO using the CT14nlo parton distribution functions.We show three different modeling approaches for fragmentation using only fragmentation fractions (green dotted), using the Peterson fragmentation function (blue dashed), and using Pythia with the QCD-inspired color reconnection scheme (red solid).For each approach, the scales were obtained using a fit to LHCb open charm data resulting in (µ F , µ R ) = (2.1, 1.6) m T (no fragmentation function), (µ F , µ R ) = (3.75,1.75) m T (Peterson fragmentation function) and (µ F , µ R ) = (2.25,1.5) m T (Pythia).In the left panel, we compare these predictions with measurements of the double differential neutral D-meson production rate obtained by LHCb at 13 TeV.We present results for three different rapidity regions, where the results at higher rapidity were scaled.In the right panel, we show the resulting number of electron neutrinos from charm hadrons decay that interact with the FASERν detector as a function of the neutrino energy.For context, we also display in black the event rate resulting from neutrinos from light hadron decays as obtained in Ref.[39].See the main text for a detailed discussion.

FIG. 3 .
FIG.3.Predictions at High Rapidity: Predictions obtained using the same QCD parameters and the same fragmentation functions as in Fig.2, but for rapidity y > 6 (left panel) and the rapidity distributions for small p T (right panel).

FASERv 13 75 FIG. 4 .
FIG. 4. Scale Variation in Collinear Factorization: Predictions using collinear factorization at NLO with different choices of scales µ F and µ R .All prediction use the CT14nlo parton distribution function, k T smearing with ⟨k T ⟩ = 1.5 GeV and Pythia with the QCD-inspired color reconnection scheme to model fragmentation.See the main text for a detailed discussion.

FIG. 5 .
FIG.5.Normalization in k T -Factorization: Predictions using k T -factorization before (blue) and after (red) applying an overall k-factor.These predictions use we use the KS (non-linear) unintegrated distribution for the low-x gluon, CT14nlo for the high-x gluon and use Pythia with the QCD-inspired color reconnection scheme to model fragmentation.See the main text for a detailed discussion.

FIG. 6 .
FIG.6.Low-x Gluon Distribution in k T -Factorization: Predictions using k T -factorization using the KS (blue dot-dashed) and CCSS (blue dashed) unintegrated distribution with a purely linear evolution as well as the KS unintegrated distribution including non-linear effects that describe saturation effects (red solid) for the low-x gluon.All predictions use a constant k-factor of 2.32, the CT14nlo parton distribution function for the high-x gluon, use Pythia with the QCD-inspired color reconnection scheme to model fragmentation.See the main text for a detailed discussion.

ddd
FIG. 7. Comparison of Charm Hadron Distribution at 13 TeV: Predictions using collinear factorization at NLO and k T -factorization.The shaded band around the NLO predictions corresponds to the scale variations shown in Fig. 4 while the shaded band around the k T -factorization prediction corresponds to a varation of the k-factor as shown in Fig. 5.In the top row, we show the p T distributions for all three charmed mesons in comparison to LHCb data.The bottom row show the rapidity distribution for D ± mesons in three transverse momentum regions.

d
FIG. 8. Comparison of Charm Hadron Distribution at 7 TeV: Transverse momentum distributions for all three charmed mesons in comparison to 7 TeV LHCb data using the same collinear factorization at NLO and k T -factorization setups as in Fig. 7.

FIG. 9 .
FIG. 9. Comparison of Forward Neutrino Distributions: Forward neutrino flux predictions using the same collinear factorization at NLO and k T -factorization setups as in Fig. 7.The energy spectra of neutrinos interacting in FASERν at Run 3 of the LHC are shown in the top row for all three neutrino flavors.Similar distributions for proposed FLARE detector at FPF during the HL-LHC era are shown in the bottom row.

FASERv 13
FIG. 10.Smearing of k T in Collinear Factorization: Predictions using collinear factorization at NLO including k T smearing for different values of ⟨k T ⟩.All predictions use fixed scales (µ F , µ R ) = (1.75, 1.25) m T , the CT14nlo parton distribution function and Pythia with the QCD-inspired color reconnection scheme to model fragmentation.See the main text for a detailed discussion.

FIG. 11 .FASERv 13 FIG. 13 .
FIG. 11.High-x Gluon Distribution in k T -Factorization: Predictions using k T -factorization based on different parton distribution of the high-x gluon.All predictions use a constant k-factor of 2.32, the KS (non-linear) unintegrated distribution for the low-x gluon and Pythia with the QCD-inspired color reconnection scheme to model fragmentation.See the main text for a detailed discussion.