Energy dependence of ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} meson production at forward rapidity in pp collisions at the LHC

The production of ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} mesons has been studied in pp collisions at LHC energies with the ALICE detector via the dimuon decay channel in the rapidity region 2.5<y<4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.5< y < 4$$\end{document}. Measurements of the differential cross section d2σ/dydpT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{d}^2\sigma /\mathrm{d}y \mathrm{d}p_{\mathrm {T}}$$\end{document} are presented as a function of the transverse momentum (pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm {T}}$$\end{document}) at the center-of-mass energies s=5.02\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=5.02$$\end{document}, 8 and 13 TeV and compared with the ALICE results at midrapidity. The differential cross sections at s=5.02\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=5.02$$\end{document} and 13 TeV are also studied in several rapidity intervals as a function of pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm {T}}$$\end{document}, and as a function of rapidity in three pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm {T}}$$\end{document} intervals. A hardening of the pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm {T}}$$\end{document}-differential cross section with the collision energy is observed, while, for a given energy, pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm {T}}$$\end{document} spectra soften with increasing rapidity and, conversely, rapidity distributions get slightly narrower at increasing pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm {T}}$$\end{document}. The new results, complementing the published measurements at s=2.76\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sqrt{s}=2.76$$\end{document} and 7 TeV, allow one to establish the energy dependence of ϕ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi $$\end{document} meson production and to compare the measured cross sections with phenomenological models. None of the considered models manages to describe the evolution of the cross section with pT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{\mathrm {T}}$$\end{document} and rapidity at all the energies.


Introduction
Measurements of production cross sections and kinematic distributions of strange hadrons represent an effective tool to investigate strangeness production in high energy hadronic collisions, testing predictions from phenomenological models inspired by quantum chromodynamics (QCD). In this context, the hardness of the specific partonic processes roughly separates two different regimes. At low transverse momentum ( p T 2 GeV/c), non-perturbative processes dominate, described by phenomenological models where different approaches may be considered, such as the rope hadronization and color reconnection mechanisms implemented in PYTHIA 8 [1]. In this regime, strangeness production in hadronic collisions with no strange valence quark component in the initial state, such as proton-proton (pp), e-mail: alice.publications@cern.ch proton-nucleus and nucleus-nucleus collisions, depends on the s-quark content of the sea-parton wave function in nucleons. At high p T ( p T 5 GeV/c), strangeness production can typically be described in terms of hard partonic scattering processes via flavor creation and excitation, and gluon splitting, for which predictions can be obtained from perturbative calculations [2]. The intermediate p T region (2 < p T < 5 GeV/c) is characterized by a smooth transition between the production mechanisms dominating at low and high p T , whose implementation in theoretical models is, however, not well defined and therefore data in this p T range are of particular importance.
In addition, strangeness production is also addressed within the phenomenological statistical model approach [3]. In small hadronic systems, in particular, the effect of canonical suppression can play a significant role in determining the relative abundances of strange and lighter flavor hadrons, with the φ meson playing a special role due to its hidden strangeness composition [4,5]. A precise theoretical characterization of this mechanism in pp collisions, however, has yet to be established for the φ meson, and predictions are currently not available.
Measurements of φ meson production in small hadronic systems like pp collisions also provide the mandatory reference for the measurements in nucleus-nucleus collisions, where a precise pp baseline is needed to single out and characterize hot medium effects affecting particle production.
In this paper, results on the transverse momentum, rapidity and energy dependence of the φ meson production cross section at forward rapidity in pp collisions at the LHC energies are presented. Results are based on the data samples collected by ALICE (A Large Ion Collider Experiment) at various energies during the LHC Run 1 and Run 2 and are compared with the predictions from the QCD-inspired models PYTHIA 8 [1], PHOJET [6,7], and EPOS [8][9][10].
The production of φ meson is studied by reconstructing the 2-body decay φ → μ + μ − with the ALICE muon spectrometer. Results are reported in the forward rapidity (y) interval 2.5 < y < 4 and for p T values ranging from 0.75 to 10 GeV/c, probing both soft and hard regimes of φ meson production.

Experimental apparatus
A full description of the ALICE detector can be found in [11,12]. The results presented in this paper have been obtained using muon pairs detected with the forward muon spectrometer, which covers the pseudorapidity region −4 < η < −2.5. Throughout this paper, the sign of η is determined by the choice of the ALICE reference system, while the kinematics of the reconstructed φ meson is referring to the positive rapidity hemisphere. The other detectors relevant for the current analysis are the Silicon Pixel Detector (SPD) of the Inner Tracking System (ITS), the V0 detector and the T0 detector.
The muon spectrometer is composed of a hadron absorber, followed by a set of tracking stations, a dipole magnet, an iron wall acting as muon filter, and a set of trigger stations. The hadron absorber, made of carbon, concrete and steel, is placed 0.9 m away from the interaction point. Its total material budget corresponds to 10 hadronic interaction lengths. The dipole magnet provides an integrated magnetic field of 3 T m in the horizontal direction, perpendicular to the beam axis. The muon tracking is provided by five tracking stations, each one composed of two cathode pad chambers. The first two stations are located upstream of the dipole magnet, the third one in the middle of its gap and the last two downstream of it. A 1.2 m thick iron wall, corresponding to 7.2 hadronic interaction lengths, is placed between the tracking and trigger detectors and absorbs the residual secondary hadrons emerging from the hadron absorber. The combined material budget of the hadron absorber and the iron wall stops muons with total momentum lower than ∼ 4 GeV/c. The muon trigger system consists of two stations, each one composed of two planes of resistive plate chambers (RPC), installed downstream of the muon filter.
The SPD consists of two silicon pixel layers, covering the pseudorapidity regions |η| < 2.0 and |η| < 1.4 for the inner and outer layer, respectively. It is used for the determination of the primary interaction vertex position. The V0 detector is composed of two scintillator hodoscopes covering the pseudorapidity regions −3.7 < η < −1.7 and 2.8 < η < 5.1. The T0 detector is composed of two arrays of quartz Cherenkov counters, covering the pseudorapidity ranges −3.3 < η < −3 and 4.6 < η < 4.9. The coincidence of a signal in both sides of either the T0 (8 TeV) or the V0 detectors (5.02 and 13 TeV) is used to define the minimum bias (MB) trigger and serves as input for the luminosity determination. It also allows for the offline rejection of beam-halo and beam-gas interactions.

Data analysis
The analysis presented in this paper is based on the data samples collected by ALICE at √ s = 5.02, 8 and 13 TeV. They complement the results already published at √ s = 2.76 [13] and 7 TeV [14].

Signal extraction
The data considered for the signal extraction were collected with a dimuon trigger, defined as the coincidence of a MB trigger and at least one pair of track segments reconstructed in the muon trigger system. The muon trigger system is configured to select muon tracks with a transverse momentum above a lowp μ T threshold, resulting in the conditions p μ T 1 GeV/c for the data sample at √ s = 8 TeV, and p μ T 0.5 GeV/c for the data samples at √ s = 5.02 and 13 TeV. 1 The number of dimuon trigger events thus selected are ∼ 2×10 7 at √ s = 5.02 TeV, ∼ 8.4×10 5 at √ s = 8 TeV and ∼ 2.8 × 10 8 at √ s = 13 TeV. Track reconstruction in the muon spectrometer is based on a Kalman filter algorithm [16,17]. Tracks reconstructed in the tracking chambers are requested to match a track segment reconstructed in the trigger chambers. In order to remove the tracks close to the acceptance borders, the muon pseudorapidity is required to be within the interval −4 < η μ < −2.5. Dimuons are formed by combining a pair of selected muon tracks, and their rapidity is explicitly imposed to be within the range 2.5 < y < 4.
The opposite-sign dimuon invariant mass spectrum contains a contribution of both uncorrelated and correlated pairs. The former mainly comes from decays of pions and kaons, which constitute the combinatorial background. This background is evaluated with an event mixing technique, in which a muon coming from an event is paired with a muon from a different event, so that the resulting muon pairs are uncorrelated by construction. This technique is described in detail in [14].
The processes contributing to the correlated dimuon mass spectrum in the low mass region, after combinatorial background subtraction, are the 2-body and Dalitz decays of the light resonances (η, ρ, ω, η , and φ), usually referred to as the hadronic cocktail, superimposed onto a continuum mainly originating from semi-muonic decays of charm and beauty hadrons. In extracting the φ → μμ signal, however, no attempt is made to describe the underlying correlated contin- 1 Because of the design of the muon trigger system, the selection on the muon transverse momentum does not correspond to a sharp value. The reported values are the ones for which the trigger efficiency is ∼ 50% [15]. The different trigger thresholds are due to the fact that the data taking at √ s = 8 TeV was optimized for high invariant mass studies. The higher p T threshold was chosen to reduce the bandwidth in the data acquisition. uum in terms of open charm and open beauty processes [18]. Instead, a fit is applied to the mass spectrum after subtraction of the hadronic cocktail with an empirical function chosen among three options, all providing data description of similar quality: a polynomial of an appropriate degree, a single or double exponential 2 and a Gaussian function whose width varies as a function of the mass as σ (m) = σ 0 (1 − e −αm ) and which will be referred to as "variable-width Gaussian" in the following. The degree of the polynomial is chosen to be the lowest one that satisfactorily fits the data, according to a statistical criterion based on the F-test, which can be briefly described as follows. The correlated continuum is first fitted with two polynomials of degree n and n + 1 respectively. The null hypothesis is that the continuum is equivalently described by the lower and the higher degree functions. This hypothesis is tested by applying the F-test to the results of the fit. If the resulting p-value is larger than 5%, the lowest degree polynomial is used, otherwise the null hypothesis is rejected and the higher degree polynomial is chosen. The procedure starts from n = 3 and is iterated until the F-test fails. The degree of the polynomial determined by means of the F-test is typically n = 4. The F-test criterion is also applied, properly adapted, to choose among the single or double exponential options.
The reconstructed opposite-sign dimuon mass spectrum is then fitted with a superposition of the hadronic cocktail and the regularized continuum discussed above: the procedure is applied independently for each of the three options considered for the empirical function describing the continuum. The free parameters of the fit are the normalization of the continuum and the η → μμγ , ω → μμ, and φ → μμ contributions, while the other processes are fixed according to the relative branching ratios or cross sections known from existing measurements [14,[19][20][21]. The mass shapes of the processes included in the hadronic cocktail are extracted through full Monte Carlo (MC) simulations that include the resolution effects induced by the apparatus. The raw number of φ for each ( p T , y) interval is then calculated as the mean of the results from the fits with the three different descriptions of the continuum. The typical reduced χ 2 of these fits is around unity. In Fig. 1 the raw dimuon invariant mass spectrum after combinatorial background subtraction is shown for √ s = 5.02, 8 and 13 TeV, in their respective p T ranges. The components of the fit are also shown. Additional details on the signal extraction procedure are reported in Ref. [18].
The differential cross section for the φ meson is calculated as where N raw φ→μμ ( p T , y) is the raw number of dimuons in the φ → μμ decay channel in a given p T , y interval as obtained from the fit procedure described above, [A × ε]( p T , y) is the corresponding product of the geometrical acceptance and the reconstruction efficiency, B R φ→μμ is the branching ratio for the φ → μμ decay, and L int is the integrated luminosity of the analyzed data sample. The [A × ε]( p T , y) factor is evaluated by means of MC simulations, where the generation of the φ → μμ process is based on a parametric generator that takes as input p T and y distributions iteratively tuned to the results of the present analyses. In detail, a first set of p T and y distributions, corresponding to the results of the measurements at √ s = 7 TeV, is taken as an input to the calculation. The resulting [A × ε]( p T , y) values are then used to correct the raw distributions obtained from the fits of the invariant mass spectra at the different energies. The corrected distributions are then used as input for another [A × ε]( p T , y) calculation, until convergence is reached. For the branching ratio of the φ → μμ decay, the value measured for the dielectron channel B R φ→ee = (2.954 ± 0.030) × 10 −4 [22] is used instead of the one of the dimuon channel, assuming lepton universality (i.e. electroweak interaction coupling to all leptons with the same strength), because the latter is affected by a larger uncertainty. The integrated luminosity is evaluated for each data set as L int = N μμ × F norm /σ MB , where N μμ is the number of analyzed opposite-sign dimuon triggered events, F norm is the inverse probability to obtain an opposite-sign dimuon trigger in a MB-triggered event, and σ MB is the MB cross section measured using the van der Meer scan method [23]. The resulting values are L int (5.02 TeV) = 1.19±0.03 pb −1 [24], L int (8 TeV) = 2.32±0.06 pb −1 [25], and L int (13 TeV) = 7.35±0.40 pb −1 [26], where the quoted uncertainties are the systematic ones, as the statistical uncertainties are negligible.

Systematic uncertainties
The systematic uncertainties on the measurement of the φ meson cross section include the following contributions: signal extraction, A × ε factor, branching ratio and integrated luminosity.
Three sources of systematic uncertainty have been considered for the evaluation of N raw φ→μμ : the choice of the function used to describe the correlated continuum, the choice of the fit range, and the uncertainty on the relative branching ratios or cross sections used to adjust some of the processes contributing to the hadronic cocktail. For the description of the correlated continuum, the three empirical functions described above were considered, providing descriptions of the data of equivalent quality.
The stability of the results under the choice of the fit range was studied by modifying the default range 2m μ < M μμ < 2 GeV/c 2 . Two alternative upper limits were considered: 2m μ < M μμ < 1.8 GeV/c 2 and 2m μ < M μμ < 2.2 GeV/c 2 .
The third contribution to the systematic uncertainty on the signal extraction was determined by varying the normalization of the ρ → μμ, η → μμγ , η → μμ and ω → μμπ 0 processes relative to the main contributions to the hadronic cocktail (η → μμγ , ω → μμ and φ → μμ), by modi-fying the relative branching ratios or cross sections. To test the sensitivity of the fit results on the normalization of the first two processes, the σ ρ /σ ω ratio was varied by 10% [21], while the σ η /σ η ratio was varied by 50%, according to the differences between σ η /σ φ obtained with calculations performed with PHOJET and with the PYTHIA 6 tunes D6T [27], ATLAS-CSC [28], Perugia 0 and Perugia 11 [29]. To account for the variation of the other two processes, the relative branching ratios of the 2-body and Dalitz decay channels of η and ω mesons were varied by one standard deviation, taking into account the values and the uncertainties reported in Ref. [22]. Table 1 Sources of systematic uncertainties in the measurement of the φ yield at various energies for rapidity-integrated (y-int.) and doubledifferential (2-diff.) analysis √ s  5 TeV  8 TeV  13 TeV The total systematic uncertainty on N raw φ→μμ is evaluated as the RMS of the values resulting from the above tests. At √ s = 5.02 and 8 TeV the variations of this contribution as a function of p T are found to be small, so their mean value is assumed as final systematic uncertainty for the signal extraction. At √ s = 13 TeV, the systematic uncertainty takes larger values in the lowest region of the measured p T range, mainly due to a progressive worsening in the description of the correlated continuum above the φ meson mass.
A potential source of systematic uncertainty on the correction for geometrical acceptance and reconstruction efficiency is associated to the choice of the input kinematic distributions used to generate the φ meson in the MC simulations. However, at √ s = 5.02 and 13 TeV these distributions are tuned to the measured data and the corrections are performed in sufficiently small p T and rapidity intervals. The resulting uncertainty, evaluated by varying the parameters of the φ meson kinematic distributions, is found to be negligible. At √ s = 8 TeV only the p T distribution was tuned to the data. To evaluate the uncertainty due to the input rapidity distribution, alternative MC simulations at √ s = 8 TeV were performed using as input the rapidity distributions measured at √ s = 5.02 and 13 TeV, and the corresponding A × ε factors were calculated. The half-difference between the two results is taken as the uncertainty due to the input rapidity distribution at √ s = 8 TeV and amounts to 3%. In addition to that, three specific sources of systematic uncertainty were considered for the reconstruction efficiency: tracking efficiency, trigger efficiency and matching efficiency. The tracking efficiency was evaluated both from data and MC simulations. To this purpose, the MC was tuned to the detector condition during data taking, in order to reproduce the uncertainties arising from correlated and anticorrelated dead areas in the muon tracker. As the same tracking algorithm is applied on MC and data, the difference observed between the two estimates is assumed as the systematic uncertainty on the tracking efficiency. The uncertainty on the trigger efficiency is mainly related to the imperfections in the description of two effects in the MC simulations: the interaction of the muons with the hadron absorber and the muon filter, and the occupancy of the trigger chambers. The uncertainty corresponding to the first effect was estimated as the difference between A × ε obtained from simulations implementing GEANT3 or GEANT4 as alternative transport codes. The uncertainty on the occupancy of the trigger chambers was evaluated comparing A × ε resulting from simulations where the MC signal was simulated either as generated or by embedding it in the environment of a real event. The uncertainty on the matching efficiency is related to the choice of the χ 2 cut used to define the matching between the tracks reconstructed in the tracking system and the track segments reconstructed in the trigger chambers, and amounts to 1% for all data samples.
The remaining contributions to the systematic uncertainty are the ones due to the branching ratio of the φ → ee decay channel (∼1%) [22] and to the integrated luminosity (2.1, 2.4 and 5.0% at √ s = 5.02, 8 and 13 TeV respectively). Two terms contribute to the uncertainty on the luminosity: the uncertainty on the visible cross section evaluated with the van der Meer scan technique and the difference between the luminosity measured with the T0 and the V0 detectors [24][25][26], while the uncertainty on the normalization factor F norm , evaluated by calculating it with two different methods, amounts to 1% for all data samples.
The systematic uncertainties listed above depend on both transverse momentum and rapidity, with the exception of the ones on the branching ratio and integrated luminosity. The sources of systematic uncertainty characterized by a dependence on transverse momentum and rapidity are classified into bin-to-bin uncorrelated (signal extraction) and bin-to-bin correlated (tracking and trigger efficiency, A × ε estimation, matching efficiency).
The contributions from the different sources of systematic uncertainties are reported in Table 1, for both the rapidityintegrated and the double-differential analyses.

Results
The differential φ meson production cross section at √ s = 5.02, 8 and 13 TeV, measured in the rapidity range 2.5 < y < 4, is shown in Fig. 2 as a function of the transverse momentum. The cross sections are fitted with a Levy-Tsallis function [30,31]: where m T = p 2 T + m 2 φ is the transverse mass and N 0 , n and T are the free parameters of the fit. While the total systematic uncertainty is shown in the spectra of Fig. 2, only the contribution coming from the signal extraction (added in quadrature to the statistical uncertainty) is considered when performing the fit, since the signal extraction is the only source of systematic uncertainty resulting in fully bin-to-bin uncorrelated fluctuations of the measured points.
The results of the fits are summarized in Table 2. The average p T , calculated using the fit functions, increases by about 20% when increasing the center-of-mass energy from 5.02 to 13 TeV.
In the same figure, the spectra at midrapidity measured by ALICE [32-34], normalized using the inelastic cross sections measured in [35], are also reported for comparison. In the bottom panels, the ratio between the fits with the Levy-Tsallis functions at forward rapidity and the data at midrapidity is shown. The cross section in the rapidity range covered by this analysis is approximately one half of the one measured at midrapidity. The p T spectra are harder at midrapidity. The difference between the slopes at forward and midrapidity is In Fig. 3, the differential cross sections are compared with the previously published results at √ s = 2.76 TeV [13] and √ s = 7 TeV [14]. The ratio to the measurement at √ s = 13 TeV (right) in the rapidity interval 2.5 < y < 4, compared with EPOS 3 [8][9][10], PHOJET [6,7] and the Monash 2013 tune of PYTHIA 8.1 [1]. The boxes represent the systematic uncertainties, the error bars the statistical uncertainties. Bottom: ratio between the measured cross section and the calculations √ s = 13 TeV is also reported in the bottom panel for a direct comparison. A hardening of the p T spectra is observed when increasing the center-of-mass energy. The values of the ratio at p T ∼ 5 GeV/c change from ∼ 0.2 for √ s = 2.76 TeV to ∼ 0.65 for √ s = 8 TeV. In Fig. 4 data are compared with the calculations performed with the models EPOS 3 [8][9][10], PHOJET [6,7] and the Monash 2013 tune of PYTHIA 8.1 [1]. At all collision energies, EPOS 3 underestimates the cross section at low p T , while it describes the data for p T > 4 GeV/c. Vice versa, PHOJET reproduces the lowp T region up to p T ∼ 2 GeV/c, but does not describe the shape of the spectra, which is predicted to be harder by the model. PYTHIA 8.1 with the Monash 2013 tune better reproduces the shape of the differential cross section at all energies. However, it still underestimates the measurement at √ s = 5.02 and 8 TeV, and reproduces the results at √ s = 13 TeV at high p T only, while it underestimates the data by about 20% at low p T . At √ s = 5.02 and 13 TeV, the p T dependence of the differential cross section was also measured in several rapidity intervals. Results are shown in Fig. 5. The p T coverage depends on the rapidity interval: in fact, at low p T , A × ε significantly increases with rapidity; on the other side, at high p T , high rapidity dimuons are more affected by statistical limitations than low rapidity ones. A significant dependence on p T and rapidity is also observed for the systematic uncertainty, namely for the contributions coming from the signal extraction and the trigger efficiency. Since this effect is related to the S/B and the data taking conditions, the impact on the results depends on the considered data sample: this explains the slight difference in the p T and rapidity coverage of the measurements at √ s = 5 and 13 TeV. The p T −differential cross sections are fitted with a Levy-Tsallis function, fixing the T parameter to the value obtained from the fit in the full range 2.5 < y < 4. Fit results are reported in Table 2, together with the average p T calculated using the fit functions. For both energies, √ s = 5.02 and 13 TeV, a moderate decrease of p T as a function of rapidity is observed. However, the relatively large uncertainties do not allow to draw any firm conclusion on this trend.
To cross check the consistency between the results shown in Figs. 2 and 5, the latter were integrated over the rapidity range 2.5 < y < 4. The differences between the two methods in the common p T region amounts to about 5%. As a comparison, the systematic uncertainty on the cross section due to signal extraction is about 4%.
The differential cross sections measured at √ s = 5.02 and 13 TeV are shown in Fig. 6 as a function of rapidity for several p T intervals, together with the corresponding values at midrapidity [32][33][34]. The calculations performed with EPOS 3, PHOJET and PYTHIA 8.1 are also plotted. None of the considered models manages to reproduce the measured rapidity spectra in all the p T ranges, neither at 5 TeV nor at 13 TeV. √ s = 13 TeV (right) in several p T intervals, compared with EPOS 3 [8][9][10], PHOJET [6,7] and the Monash 2013 tune of PYTHIA 8.1 [1]. The ratio of the data to the lowest p T interval is shown in the bottom panel. The boxes represent the systematic uncertainties, the error bars the statistical uncertainties In the bottom panel of Fig. 6, the ratio of the rapidity distributions in a given p T interval to the lowest is shown, scaled such that the ratio at midrapidity is set to unity. A moderate narrowing of the rapidity distribution is observed when increasing the transverse momentum. This effect depends on the collision energy, being stronger for the lowest √ s. The φ meson production cross section integrated in the p T range 1.5 < p T < 5 GeV/c, common to the five energies, is plotted as a function of √ s in Fig. 7. Results are compared with EPOS 3, PHOJET and PYTHIA 8.1 with the Monash 2013 tune. In the p T interval considered for this study, the evolution of the cross section with the collision energy is well described by PHOJET. Both EPOS 3 and PYTHIA 8.1/Monash 2013 underestimate the absolute values by a factor ranging from about 1.2 to 1.7, while reproducing the trend as a function of the center-of-mass energy. The differences between the measurement and the calculations are mainly due to the overestimation of the cross section at the lowest p T accessible to the measurement.

Conclusions
The φ meson production cross section was measured in pp collisions at the center-of-mass energies √ s = 5.02, 8, 13 TeV in the forward rapidity region 2.5 < y < 4, complementing the previously published results at √ s = 2.76 and 7 TeV. The p T spectra are well described by a Levy-Tsallis function. A  Fig. 7 φ meson production cross section as a function of √ s for 1.5 < p T < 5 GeV/c and 2.5 < y < 4, compared with EPOS 3 [8][9][10], PHOJET [6,7] and the Monash 2013 tune of PYTHIA 8.1 [1]. The boxes represent the systematic uncertainties, the error bars the statistical uncertainties hardening of the p T -differential cross section with the collision energy is observed, as is evinced from the comparison between the average values of the transverse momentum or from the ratios between the differential cross sections as a function of p T . At each energy, the p T spectra at midrapidity are harder than the corresponding ones at forward rapidity. Results were compared with the predictions from PYTHIA 8.1-Monash 2013, EPOS 3 and PHOJET. PHOJET reproduces the cross section at low p T , while EPOS 3 better approaches the data for p T > 4 GeV/c. PYTHIA 8.1 with the Monash 2013 tune fairly describes the shape of the p Tdifferential cross section at all energies and reproduces the results at √ s = 13 TeV, but underestimates the measurement at lower energies.
At √ s = 5.02 and 13 TeV, a double differential study of the φ meson production cross section was performed as a function of p T and rapidity. None of the calculations manages to reproduce the rapidity spectra both at low and high p T . A small decrease of the p T value was observed with increasing rapidity, although with relatively large uncertainties. Analogously, the dependence of the cross section on rapidity appears to be slightly narrower when going towards higher p T values, thus showing that the correlation between p T and rapidity cannot be neglected in this rapidity range.
Acknowledgements The ALICE Collaboration would like to thank all its engineers and technicians for their invaluable contributions to the construction of the experiment and the CERN accelerator teams for the outstanding performance of the LHC complex. The ALICE Collaboration gratefully acknowledges the resources and support provided by all Grid centres and the Worldwide LHC Computing Grid (WLCG) collaboration. The ALICE Collaboration acknowledges the following funding agencies for their support in building and running the ALICE detector: A. I. Alikhanyan National Science Laboratory  Data Availability Statement This manuscript has associated data in a data repository. [Authors' comment: Manuscript has associated data in a HEPData repository at https://www.hepdata.net/.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .