Measurement of $\phi$-meson production at forward rapidity in $p$$+$$p$ collisions at $\sqrt{s}$=510 GeV and energy dependence of $\sigma_\phi$ from $\sqrt{s}$=200 GeV to 7 TeV

The PHENIX experiment at the Relativistic Heavy Ion Collider has measured the differential cross section of $\phi$(1020)-meson production at forward rapidity in $p$$+$$p$ collisions at $\sqrt{s}=$510 GeV via the dimuon decay channel. The partial cross section in the rapidity and $p_T$ ranges $1.2<|y|<2.2$ and $2<p_T<7$ GeV/$c$ is $\sigma_\phi=[2.28 \pm 0.09\,{\rm (stat)} \pm 0.14\,{\rm (syst)} \pm 0.27\, {\rm (norm)}] \times 10^{-2}$~mb. The energy dependence of $\sigma_\phi$ ($1.2<|y|<2.2, \; 2<p_T<5$ GeV/$c$) is studied using the PHENIX measurements at $\sqrt{s}=$200 and 510 GeV and the Large-Hadron-Collider measurements at $\sqrt{s}=$2.76 and 7 TeV. The experimental results are compared to various event generator predictions ({\sc pythia6, pythia8, phojet, ampt, epos3,} and {\sc epos-lhc}).


I. INTRODUCTION
The φ(1020)-vector-meson production in p+p collisions was intensively studied by various experiments at different colliding energies and in different rapidity ranges [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]. It is the lightest bound state of s ands quarks and is considered a good probe to study strangeness production in p+p collisions. Production of φ mesons from an initial nonstrange colliding system, such as p+p collisions, is substantially suppressed in comparison to ω and ρ vector mesons due to the Okubo-Zweig-Iizuka (OZI) rule [19][20][21]. The φ-meson production at low transverse momentum is dominated by soft processes and is sensitive to the hadronization mechanism, while hard processes become dominant at higher transverse momentum. In p+p collisions, the production of strangeness is in general not well described by generators such as pythia, which tend to underestimate the production of strange particles [10,[22][23][24]. The study of φ-meson production in p+p collisions is an important tool to study quantum chromodynamics (QCD), providing data to tune phenomenological QCD models in which an interplay is mandatory between perturbative QCD calculations, used in particular for hard parton production dominant at higher p T , and phenomenological QCD models, needed to describe the nonperturbative hadronization into strange hadrons like the φ meson.
In addition, recently, a long-range near-side angular correlation was observed in p+p collisions at Large-Hadron-Collider (LHC) energies [25][26][27], which led to the observation of collectivity in p+p collisions [28]. This observation generated various explanations [29], including those based on the color-glass-condensate (CGC) model [30], and collective hydrodynamic flow [31] or color reconnection [32,33]. Being the heaviest easily accessible meson made of light quarks, φ-meson production provides the largest lever arm accessible to study effects that scale with mass, as should be the case for collective effects [34]. * Deceased † PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov The study of φ-meson production in p+p collisions can be an important tool to gain insight into new phenomena, such as long-range angular correlations, that would have a direct impact in the field of relativistic heavy-ion collisions. The φ-meson production is an excellent observable to probe the strangeness enhancement in the quark-gluon plasma created in heavy-ion collisions [35][36][37].
We report the φ-meson-production cross section measured in p+p collisions at √ s = 510 GeV. The analysis uses a data sample of 144.6 pb −1 of integrated luminosity obtained by the PHENIX experiment in 2013. The cross section is averaged over the rapidity (y) interval 1.2 < |y| < 2.2 and reported in several bins of transverse momentum (p T ) in the range 2 < p T < 7 GeV/c. The results are compared to several model predictions [24,34,[38][39][40][41] and to the measurements previously reported by the PHENIX experiment at √ s = 200 GeV [15] and by the LHC experiments measuring the φ-mesonproduction cross section at forward rapidity at √ s = 2.76 and 7 TeV [10][11][12][13]17]. Measurements from experiments at the Relativistic Heavy Ion Collider (RHIC) and the LHC allow extracting the energy dependence of the φ-meson-production cross section in the rapidity range 1.2 < y < 2.2, which provide information to further constrain model predictions.

II. EXPERIMENTAL SETUP
A complete description of the PHENIX detector can be found in Ref. [42]. The results presented here are obtained by measuring the φ meson via its µ + µ − decay channel using both PHENIX muon spectrometers covering forward and backward pseudorapidities, 1.2 < |η| < 2.2, and the full azimuth.
Each muon arm spectrometer comprises hadron absorbers, a muon tracker (MuTr), which resides in a radial field magnet, and a muon identifier (MuID). The absorbers are situated in front of the MuTr to provide hadron (mostly pion and kaon) rejection and are built of 19 cm of copper, 60 cm of iron, and 36.2 cm of stainless steel. The MuTr comprises three sets of cathode strip chambers in a radial magnetic field with an inte-grated bending power of 0.8 T.m. The final component is the MuID, which has five alternating steel absorbers and Iarocci tubes to further reduce the number of punchthrough hadrons misidentified as muons. Muon candidates are identified by reconstructed tracks in the MuTr matched to MuID tracks that penetrate through to the last MuID plane.
Another detector system relevant to this analysis is the beam-beam counter (BBC), comprising two arrays of 64Čerenkov counters, located on both sides of the interaction point and covering the pseudorapidity 3.1 < |η| < 3.9. The BBC system is used to measure the p+p collision vertex position along the beam axis (z vtx ) with 2 cm resolution and to provide the minimum bias (MB) trigger.

III. DATA ANALYSIS
The results presented here are based on the data sample collected by PHENIX during the 2013 p+p run at √ s = 510 GeV. The BBC counters provide the MB trigger, which requires at least one hit in each of the BBCs. Events, in coincidence with the MB trigger, containing a muon pair within the acceptance of the spectrometer are selected by the level-1 dimuon trigger requiring that at least two tracks penetrate through the MuID to its last layer. A total of 5.3 × 10 8 dimuon triggered events are recorded, which corresponds to a sampled integrated luminosity of 144.6 pb −1 .

A. Raw yield extraction
A set of quality assurance cuts is applied to the data to select p+p events and muon candidates as well as to improve the signal-to-background ratio. Good p+p events are selected by requiring that the collision occurs in the fiducial interaction region |z vtx | < 30 cm as measured by the BBC. No selection is made on the event's charged particle multiplicity. The MuTr tracks are matched to the MuID tracks at the first MuID layer in both position and angle. In addition, the track is required to have more than a minimum number of possible hits in the MuTr (12 out of the maximum 16) and MuID (6 out of the maximum 10), and cuts on the individual track χ 2 values are applied. Furthermore, there is a minimum allowed single muon momentum along the beam axis, p z , which is reconstructed and energy-loss corrected at the collision vertex, of 2.4 GeV/c corresponding to the momentum cut effectively imposed by the absorbers. Finally, a cut on the χ 2 of the fit to the common vertex of the two candidate tracks near the interaction point is made.
The invariant mass distribution is formed by combining muon candidate tracks of opposite charges. This unlike-sign dimuon spectrum is composed of correlated and uncorrelated pairs. In the low-mass region (below ≈1.5 GeV/c 2 ) the correlated pairs arise from the two- body and Dalitz decays of the light neutral mesons η, ρ, ω, η and φ as well as semi-muonic decays of correlated charmed hadrons (and beauty in a negligible contribution). The uncorrelated pairs are mainly coming from semi-muonic decays of pions and kaons and punchthrough hadrons, and form the so-called combinatorial background. The ratio of φ-meson signal over combinatorial background is of the order of 0.7. This combinatorial background is estimated using two methods: the first one derives the combinatorial background from the distribution formed within the same event by the muon candidates of the same sign (like-sign pairs); and the second one derives the combinatorial background from the pairs formed by muon candidates of opposite charges (unlike-sign pairs) coming from different events (mixedevent). The normalization of the mass distribution of the combinatorial background using the same-event like-sign dimuon distributions (N ++ and N −− ) is calculated as: The mixed-event like-sign dimuon mass distribution is normalized to the same-event like-sign combinatorial background distribution in the invariant mass range 0.2 − 2.5 GeV/c 2 . This factor is then used to normalize the mixed-event unlike-sign dimuon mass distribution. Figure 1 shows the unlike-sign dimuon spectrum together with the combinatorial background estimated by both methods that agree within 15% in the invariant mass range of interest (0.8 < M µµ < 1.3 GeV/c 2 ).
The signal invariant mass spectrum is extracted by first subtracting the uncorrelated combinatorial background spectra from the unlike-sign spectra. The signal spectra are then fitted to extract the φ contribution. The mass resolution of both muon spectrometers is estimated using Monte Carlo simulation to be 93 (94) MeV/c 2 for the lowest p T bin (2 < p T < 2.5 GeV/c) and up to 114 (111) MeV/c 2 for the highest p T bin (5 < p T < 7 GeV/c) for the negative (positive) pseudorapidity muon spectrometer. Those resolutions being greater than the natural widths of the φ and ω, the two-body decay of φ  and ω contributions are described by Gaussians while the two-body decay of the ρ-meson contribution is described by a Breit-Wigner distribution convoluted with a Gaussian. The contribution from ρ dimuon decay is fixed by the assumption that the production cross section of ρ and ω are related such as σ ρ = 1.15 × σ ω , as measured in Ref. [12] and used in previous PHENIX analysis related to φ-meson production in the dimuon decay channel [15,43,44]. To evaluate the shape of the correlated background, a pythia [45] MB simulation followed by geant3 [46] transport and detector response simulation of the PHENIX detector is performed. The correlated background distribution is found to be well described by an exponential plus a polynomial of first order (χ 2 /ndf ≤ 1). To summarize, eight free parameters are needed to describe the signal spectrum: two parameters for the φ and (ω + ρ) signal normalizations, two parameters to describe relative changes of Gaussian widths and central masses with respect to simulation estimates and four parameters to describe the correlated background distribution and its normalization. The starting values of the free parameters describing the shapes of the different distributions are taken to be the ones from the Monte Carlo simulation. Figure 2 shows the fit results for the entire p T range at backward rapidity. Extracted peak positions and widths are found to be in good agreement with Monte Carlo simulations.

B. Detector acceptance and reconstruction efficiency
The product of detector acceptance and reconstruction efficiency, A rec , of dimuon decays of φ mesons is determined by the full event reconstruction of the φmeson signal run through a full geant3 simulation of the 2013 PHENIX detector setup, and embedded in MB real-data. The p T distribution of the simulated φ-meson signal is iteratively re-weighted to match the data p T dis- and backward (−2.2 < y < −1.2) muon spectrometers (a) in the pT -rapidity plane and (b) integrated in rapidity per spectrometer for each pT bin considered in the analysis.
tribution, the initial p T distribution being obtained from pythia6 [45] using tune ATLAS CSC [38]. The embedded simulated events are then reconstructed in the same manner as data with the same cuts applied as in the real data analysis. The A rec factor is extracted from the simulation as the ratio of reconstructed φ distribution over the generated one in the same kinematic range. Figure 3 shows the A rec as a function of φ-meson p T and rapidity. The main sources of the relative difference between both spectrometers A rec are different detection efficiencies of the MuTr and MuID systems and different amount of absorber material. Br φ→µ + µ − 6.6%

C. Differential cross section extraction
The p T -dependent differential cross section is calculated according to: [47]. N raw is the extracted φ raw yield for each p T bin, N BBC MB = 4.16×10 12 is the number of sampled MB events. The BBC trigger samples a cross section of σ BBC pp = 32.5±3.2 mb in p+p collisions, according to Vernier scans, however, it samples a larger fraction of the cross section when the collision includes a hard scattering process [48]. Studies with high p T π 0 yields show an increase of the luminosity scanned by the BBC by a factor of 1/ BBC , BBC = 0.91±0.04 [49]. The inelastic cross sections given by pythia8 [50] for √ s = 500 and 510 GeV p+p collisions differ by 0.3%, therefore no correction or additional systematic uncertainty is added.

D. Systematic uncertainties
The main source of systematic uncertainties in the signal extraction comes from the uncorrelated and correlated background distributions used. To estimate this uncertainty, the extracted φ raw yields are compared using the following two methods; (1) the mixing and like-sign pair methods are separately used for subtraction of uncorrelated background and (2) the correlated background is fit by an exponential plus first-order polynomial and by an exponential plus second-order polynomial. The extracted φ raw yields are consistent among all different fit trials. The quadratic mean of the raw yields extracted from the trials is used as the central value, and the uncertainty on the central value is the quadratic mean of the uncertainties of all the trials. Table I summarizes the  systematic uncertainties. Type-A is a point-to-point uncorrelated uncertainty which allows the data points to move independently with respect to one another and are added in quadrature with statistical uncertainties. A systematic uncertainty equal to the difference between the central and the extreme values of the extracted yields accounts for the systematic uncertainty related to the background description as a whole. The systematic uncertainty associated with the signal extraction method ranges from 3 to 23%, depending on the p T bin and the muon spectrometer considered (negative/positive rapidity).
Type-B is a point-to-point correlated uncertainty which allows the data points to move coherently. To evaluate the A rec systematic uncertainty, different p T and rapidity input distributions of the simulated φ mesons are used. The p T distribution is allowed to vary over the range of the data statistical uncertainty (statistical plus Type-A systematics uncertainties added in quadrature, see above), yielding an up to 8% uncertainty. The rapidity distribution shapes given by five generator models (pythia6, pythia8, phojet, epos3 and epos-lhc) are used as input rapidity distributions of the simulated φ mesons, resulting in up to 5% uncertainty. The relative systematic uncertainty of acceptance caused by the fluctuation of vertex width is estimated to be 3.5% [51]. A 4% uncertainty from the measured MuID tube efficiency and a 2% uncertainty from MuTr chamber efficiency are assigned [15]. Simulation parameters are adjusted in order to reproduce the tracking efficiency observed in the data. While the relative tracking efficiency is validated using J/ψ → µµ data, data-driven evaluation of the absolute tracking efficiency are not available. Therefore, we assign 10% uncertainty for the absolute tracking efficiency as a conservative value [51].
Finally, Type-C is an overall normalization uncertainty, which allows the data points to move together by a common multiplicative factor. Type-C is composed of 10% uncertainty assigned for the BBC cross section and efficiency uncertainties and a 6.6% uncertainty from the measurement of BR φ→µ + µ − .

IV. RESULTS
The p T -differential cross section is calculated independently for each muon arm, then the results are combined using the best-linear-unbiased-estimate method [52]. The p T integrated (2 < p T < 7 GeV/c) cross section dσ φ /dy is given in Table II. Results obtained using the two muon spectrometers are consistent within uncertainties. Combining both arm results, the integrated cross section in the kinematic range 2 < p T < 7 GeV/c and 1.2 < |y| < 2.2 is σ φ = 2.28 ± 0.09 (stat) ± 0.14 (syst) × 10 −2 mb, to which a 12% normalization uncertainty applies.
The φ-meson-differential cross section as a function of p T measured in p+p collisions at √ s = 510 GeV is shown in Fig. 4 and listed in Table III. The data points are bin shifted in p T using the Lafferty and Wyatt method [53] TABLE II. The φ-meson-production cross section dσ φ /dy in p+p collisions at √ s = 510 GeV integrated in the transverse momentum range 2 < pT < 7 GeV/c. The first uncertainty represents the statistical and Type-A systematic uncertainties, while the second is the systematic uncertainty of Type-B and the third one is the additional ±12% Type-C normalization systematic uncertainty.
y range dσ φ /dy (mb) 1.2 < y < 2.2 (2.13 ± 0.14 ± 0. 16 (2.28 ± 0.09 ± 0.14 ± 0.27) × 10 −2 TABLE III. The φ-meson-differential-production cross section d 2 σ φ /dpT dy for 1.2 < |y| < 2.2 in p+p collisions at √ s = 510 GeV.pT is the pT at which the data point is plotted (see text for details). The first uncertainty represents the statistical and Type-A systematic uncertainties, while the second is the systematic uncertainty of Type-B and the third one is the additional ±12% Type-C normalization systematic uncertainty. to correct for the finite width of the p T bins. The data are fitted by a Tsallis function [54] with a resulting χ 2 /ndf = 0.66. The results are compared to calculations performed using six different generator models: pythia6 [45] using tune ATLAS CSC [38], pythia8.210 [50] using tune Monash2013 [24], phojet 1.12 [39], epos3.117 [34], epos-lhc [40] and ampt v1.26 [41]. Data and models are compared as the ratio of the model prediction over the Tsallis fit of the data.
The ampt simulation is done with the default ampt model version 1.26 (without string melting), where the initial conditions are determined by hijing [55]. Parton scattering is done using the Zhang's parton-cascade (ZPC) model [56]. The hadronization is accomplished using the Lund string fragmentation model [57,58]. The final state hadronic interactions are based on the "a relativistic transport" (ART) model [59]. We used the set of parameters tabulated in Ref. [60] describing both the charged particle distribution and elliptic flow measured in Au+Au collisions at RHIC. The Lund string fragmentation parameters are a = 0.5 and b = 0.9 GeV −2 , the QCD coupling constant is α s = 0.33, and the screening mass is µ = 3.2 fm −1 , leading to a parton-scattering cross section of 1.5 mb. Besides their production from the fragmenta- tion of excited strings in the initial collisions, φ mesons can also be produced and absorbed from hadronic matter via various hadronic reactions (baryon-baryon, mesonbaryon and meson-meson scatterings) [41].
The epos3 model includes, in addition to the description of the initial scattering based on a Gribov-Regge approach [61], a viscous hydrodynamic expansion of the created system followed by a hadronization phase and a final state hadronic cascade using the urqmd model [62,63]. In epos3, the hydrodynamic evolution and the hadronic cascade can be turned on or off, separately. The socalled "Full" version of epos3 includes hydrodynamic expansion of the created system followed by a final state hadronic cascade. The epos3 "No-Casc" version does not include the final state hadronic cascade and "No-Hydro/No-Casc" has both hydrodynamic and the final state hadronic cascade turned off. The epos-lhc calcu-    lation presented in Fig. 4 is performed including a parameterized viscous hydrodynamic expansion of the created partonic system. As shown in panels (b) and (c) of Fig. 4, the experimental data are better reproduced by the ampt model and by epos3 without the hadronic cascade. The epos3 "Full" and epos-lhc overestimate the φ-meson production, and phojet and pythia models tend to underestimate it by a factor of two. A previous study of Monash2013 tune of pythia8 showed that the calculated transversemomentum spectra of φ mesons is overestimating the experimental data at very soft momenta (below ∼500 MeV/c) and underestimates it at higher momenta, the overall yield of φ mesons being correctly reproduced [24].
Additional calculations using the ampt model with string melting (version 2.26) were performed. The φmeson-production yield was found to be a factor of two higher than the one extracted using the default ampt model with approximately the same p T dependence. For clarity, those calculations are not shown in Fig. 4.

V. ENERGY DEPENDENCE OF φ-MESON PRODUCTION
The PHENIX experiment previously measured the φmeson cross section at forward rapidity and for 1 < p T < 7 GeV/c in p+p collisions at √ s = 200 GeV [15]. At the LHC, the ALICE experiment measured the φ-mesonproduction cross section via its dimuon decay channel in p+p collisions at forward rapidity 2.5 < y < 4.0 and for 1 < p T < 5 GeV/c at √ s = 2.76 TeV [17] and 7 TeV [12]. Measurement of the φ-meson production was also performed via the K + K − decay channel at midrapidity |y| < 0.5 and for 0.4 < p T < 6 GeV/c at √ s = 7 TeV [13]. The LHCb experiment measured the inclusive φ-meson-production cross section in the K + K − decay channel in the kinematic range 2.44 < y < 4.06 and 0.6 < p T < 5 GeV/c in p+p collisions at √ s = 7 TeV [10]. Figures 5-8 show comparisons between d 2 σ φ /dp T dy measurements at forward rapidities done by PHENIX at √ s = 200 GeV [15], by ALICE at √ s = 2.76 TeV [17] and 7 TeV [12] and by LHCb at √ s = 7 TeV [10], respectively, along with model predictions. The ampt model is in good agreement with the measured cross sections at both RHIC energies, but overestimates the production cross section at LHC energies, especially at 7 TeV. The pythia6 and phojet calculations at LHC energies are in better agreement with the data than at RHIC energies, where the models underestimate the measured production cross section. The pythia8 prediction underestimates the cross section for all four energies.
Panel (c) of Figs. 4-8 show the comparison between the measurements fitted by a Tsallis function and epos3 using three different model settings (see above for details). The comparison of those results reveals the effect of the hydrodynamic expansion of the partonic system created in p+p collisions and of the final state hadronic cascade on the φ-meson production. The hydrodynamic evolution does not impact the φ-meson production at RHIC energies ("No-Casc" and "No-Hydro/No-Casc" curves are almost identical on panel (c) of Figs 4-8. A significant effect appears at √ s = 2.76 TeV and becomes stronger at 7 TeV where the φ-meson-production cross section increases by a factor of two for the p T range 1-3 GeV/c when turning on the hydrodynamic evolution. The same behavior was already observed for the production of Λ 0 , K s and Ξ ± in p+p collisions at 7 TeV [34], showing that the flow effects increase with the mass of the particle. The final state hadronic cascade using the urqmd model enhances the φ-meson-production cross section in the entire p T range and for all collision energies. The epos3 "No-Casc" is the best configuration to reproduce the experimental data over the full collision energy range, while the addition of the urqmd hadronic cascade overestimates the φ-meson production compared to the experimental data.
In the following, the φ-meson cross sections in the forward rapidity range 1.2 < y < 2.2 at the different measured energies (0.2, 0.51, 2.76 and 7 TeV) are presented. The p T range is fixed to 2 < p T < 5 GeV/c which is the common range of all experimental measurements.
The cross sections measured by PHENIX in the kinematic range 1.2 < y < 2.2 and 2 < p T < 5 GeV/c are: where the uncertainties correspond to the quadratic sums of the statistical and systematic uncertainties.
The rapidity domains of the LHC measurements are different from those of PHENIX. Accordingly, to compare with PHENIX measurements the LHC measurements are extrapolated to the same rapidity coverage (i.e. 1.2 < y < 2.2). The procedure followed here is to fit the LHC data points using the dσ φ /dy shapes obtained using the different models mentioned above, the only free parameter being the normalization of the simulated dσ φ /dy distributions. Figure 9 shows the LHC p T integrated data points overlaid on the dσ φ /dy distributions obtained using pythia6, pythia8, phojet, epos3, epos-lhc and ampt models at √ s = 7 TeV (a) before and (b) after the minimization procedure.
The LHC dσ φ /dy at 1.2 < y < 2.2 is calculated as the quadratic mean of the dσ φ /dy from each of the model fits. The difference between the mean and the extreme value is taken as a systematic uncertainty, due to the rapidity shifting procedure, and added in quadrature to the experimental uncertainties. This uncertainty is 22.1% for the 2.76 TeV measurement and 15.5% at 7 TeV. The obtained cross sections in 1.2 < y < 2.2 and 2 < p T < 5 GeV/c at LHC energies are: • σ φ (2.76 TeV) = (1.15 ± 0.28) × 10 −1 mb, • σ φ (7 TeV) = (2.23 ± 0.35) × 10 −1 mb.  [10,12,13] and results of five simulations using pythia6, pythia8, phojet, ampt, epos3 and epos-lhc generator models. (b) Results of the fit of the measurements using the simulated dσ φ /dy shapes with normalization as the only free parameter. Figure 10 shows the energy dependence of the partialφ-meson-production cross section integrated in 1.2 < y < 2.2 and 2 < p T < 5 GeV/c in p+p collisions compared to pythia6, pythia8, phojet, ampt, epos3 and eposlhc model predictions.
The experimental measurements follow a power-law versus the colliding energy defined as σ φ (s) ∝ s n , with n = 0.43±0.03 (black dotted line in Fig. 10). The χ 2 /ndf of the power-law fit is 0.19.
The phojet generator reproduces the partial φ-meson cross section correctly for LHC energies, but completely fails at RHIC energies. On the other hand, the ampt model performs well at lower energies but overshoots the experimental data at 7 TeV. pythia6 shows an energy dependence following a power law with exponent n = 0.43, comparable to that of the data, but underestimates the cross section by ∼30%. Accounting for hydrodynamic evolution of the partonic system makes epos3 qualitatively and quantitatively better consistent with the data from both RHIC and LHC. The increasing effect of the hydrodynamic evolution of the system on the φ-meson production as the energy increases can clearly be seen in Fig. 11. Also, the φ-meson enhancement caused by the hadronic cascade is approximately constant over the whole energy range, ≈20-30%.  . Partial-φ-meson-production cross section in 1.2 < y < 2.2 and 2 < pT < 5 GeV/c in p+p collisions versus the center-of-mass energy √ s compared to epos3 model using different options. The LHC data points are interpolated at the PHENIX forward rapidity, see text for details.
In epos3, when the hydrodynamic evolution is turned off the hadrons are produced via string decays. On the other hand, when hydrodynamic calculation is included, the various string segments originating from the initial Pomerons are separated into two collections named "core" and "corona". The "core" part will experience the hydrodynamic evolution while the segments in the "corona" will leave the bulk matter and decay to hadrons. String segments are placed in the "core" or "corona" depending on their transverse momenta and on the local string density [34]. After its hydrodynamical evolution, the "core" hadronizes following the Cooper-Frye freeze-out procedure. Figure 12 shows the "core" and the "corona" contributions to the production of φ mesons in p+p collisions for the four energies studied in this work. The contribution of the "core" part increases with the colliding energy, being negligible compared to the "corona" contribution at RHIC energies and of the same order of magnitude at LHC energies for 1 < p T < 3 GeV/c. The difference in the shape of the p T distributions between the "core" and the "corona" part (shift from low to intermediate p T ) is due to the fact that in the "core" the φ mesons are produced from "fluid cells characterized by a radial flow velocities" [34]. The heavier the particle is the more transverse momentum it receives from this mechanism.

VI. SUMMARY AND CONCLUSIONS
In summary, the φ-meson-production differential cross section is measured in p+p collisions at √ s =510 GeV in the kinematic range 1.2 < |y| < 2.2 and 2 < p T < 7 GeV/c. The cross section integrated in p T and averaged over positive and negative rapidities is σ φ = [2.28 ± 0.09 (stat) ± 0.14 (syst) ± 0.27 (norm)] × 10 −2 mb. The measured p T -differential cross section is compared to various model predictions based on pythia6, pythia8, phojet, ampt, epos3 and epos-lhc generators. The default ampt model and the epos3 model without hadronic cascade provide the best description of the data.
The energy dependence of the φ-meson-production cross section is studied in the kinematic range 1.2 < y < 2.2 and 2 < p T < 5 GeV/c, shifting LHC measurements to the same rapidity range as PHENIX measurements. The epos3 model shows that the addition of the hydrodynamic evolution of the system induces an enhancement of the φ-meson production at the LHC energies for 1 < p T < 3 GeV/c whereas no effect is seen for RHIC energies. The LHC measurements tend to favor the scenario with hydrodynamic evolution of the system included in epos3 showing a possible hint of collective effects in p+p collisions at high energy.
The epos3 model shows that the hydrodynamic flow induces a shift from low to intermediate p T of the produced φ mesons. A similar effect is obtained from tuning the color reconnection mechanism in pythia8 [32,33]. The study of the p T as a function of the charged particle multiplicity produced in p+p collisions and its evolution versus the colliding energy would be a relevant observable of such effect, and would allow to discriminate between alternative models. In addition to the already published data at √ s = 2.76 and 7 TeV regarding the production of φ mesons at forward rapidity, the LHC experiments took data in p+p collisions at 5, 8, and recently 13 TeV where the effect should be even larger.