1 Introduction

The primary goals of high-energy heavy-ion (A–A) collisions are to create a system of deconfined quarks and gluons known as quark–gluon plasma (QGP) and to study its properties [1,2,3,4]. Asymmetric collision systems like proton-nucleus (p–A) and deuteron-nucleus (d–A) can be considered as control experiments where the formation of an extended QGP phase is not expected. These collision systems are used as baseline measurements to study the possible effects of cold nuclear matter and disentangle the same from hot dense matter effects produced in heavy-ion collisions [5,6,7,8,9,10,11,12,13,14]. In addition, p–A collisions at Large Hadron Collider (LHC) energies enable probing the parton distribution functions in nuclei at very small values of the Bjorken x variable, where gluon saturation effects may occur [15,16,17]. Recent measurements in high-multiplicity pp, \(\text{ p }{-}\text{ Pb } \), p–Au, d–Au, and \(^{3}\)He–Au collisions at different energies have shown features such as anisotropies in particle emission azimuthal angles, strangeness enhancement, and long-range structures in two-particle angular correlations on the near and away side, which previously have been observed in nucleus-nucleus collisions [18,19,20,21,22,23,24,25,26,27,28,29]. The origin of these phenomena in small systems is not yet fully understood. A systematic study of multiplicity and rapidity dependence of hadron production allows us to investigate the mechanism of particle production and shed light on the physics processes that contribute to the particle production [15]. Similar studies have been reported by the experiments at the LHC [9, 10, 16, 17] and Relativistic Heavy Ion Collider (RHIC) [6, 7, 19,20,21]. The mechanism of hadron production may be influenced by different effects such as nuclear modification of the parton distribution functions (nuclear shadowing) and possible parton saturation, multiple scattering, and radial flow [16, 30,31,32]. These effects are expected to depend on the rapidity of the produced particles. In \(\text{ p }{-}\text{ Pb } \) collisions, one can expect that the production mechanism may be sensitive to different effects at forward (p-going) and backward (Pb-going) rapidities [9, 10, 16, 17, 32, 33]. The partons of the incident proton are expected to undergo multiple scattering while traversing the Pb-nucleus. It is thus interesting to study the ratio of particle yields between Pb- and p-going directions, represented by the rapidity asymmetry \((Y_{\textrm{asym}})\) defined as:

$$\begin{aligned} Y_{\textrm{asym}}(p_{\textrm{T}}) =\frac{\left. \frac{{\textrm{d}}^{2}N}{{\textrm{d}}p_{\textrm{T}} {\textrm{d}}y}\right| _{-0.3< y< 0}}{\left. \frac{{\textrm{d}}^{2}N}{{\textrm{d}}p_{\textrm{T}}{\textrm{d}}y}\right| _{0< y < 0.3}} \end{aligned}$$
(1)

where \(\mathrm{{d}}^{2}{N}/ \mathrm{{d}}{p}_{\mathrm{{T}}}\mathrm{{d}}{y}|_{-0.3< y < 0}\) is the particle yield in the rapidity (y) interval \(-0.3< y < 0,\) considered as the Pb-going direction, and \(\mathrm{{d}}^{2}{N}/\mathrm{{d}}{p}_{\mathrm{{T}}}\mathrm{{d}}{y}|_{0< y < 0.3}\) is the particle yield in the rapidity interval \(0< y < 0.3,\) corresponding to the p-going direction. From the experimental point of view, the \(Y_{\textrm{asym}}\) is a powerful observable because systematic uncertainties cancel out in the ratio and hence it can better discriminate rapidity-dependent effects among models [16,17,18]. Gluon saturation effects at low Bjorken x values [7, 18] may affect the transverse momentum distribution of hadron production at large rapidities in the p-going direction in \(\text{ p }{-}\text{ Pb } \) collisions at LHC energies. The gluon saturation effects depend on the colliding nuclei and rapidity as A\(^{1/3}e^{\lambda y}\), where A represents the mass number [18], and \(\lambda \) is a parameter whose value lies between 0.2 and 0.3, and is obtained from fits to the HERA measurements [8]. The effect of rapidity dependence on particle production is tested by measuring the ratios of integrated yield (dN/dy) and mean transverse momentum \(\left( \langle p_{\textrm{T}} \rangle \right) \) at given y to the values at \(y=0\), i.e. denoted as (dN/dy)/(dN/d\(y)_{y=0}\) and \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}} \rangle _{y=0}\). It is also important to study the variation with rapidity of the nuclear modification factor between central and non-central collisions. This factor \(\left( Q_{\textrm{CP}} (p_{\textrm{T}})\right) \) is defined as

$$\begin{aligned} {Q}_{\textrm{CP}} (p_{\textrm{T}})= \left. \frac{\frac{{\textrm{d}}^{2}N}{{\textrm{d}}p_{\textrm{T}}{\textrm{d}}y}}{\langle N_{\textrm{coll}}\rangle }\right| _{\text {HM}} \Bigg / \left. \frac{\frac{{\textrm{d}}^{2}N}{{\textrm{d}}p_{\textrm{T}}{\textrm{d}}y}}{\langle N_{\textrm{coll}} \rangle }\right| _{\text {LM}}, \end{aligned}$$
(2)

where \(\langle N_{\textrm{coll}} \rangle \) is the average number of nucleon–nucleon collisions in low-multiplicity (LM) and high-multiplicity (HM) events, respectively. The multiplicity dependence of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) meson production at midrapidity was studied in pp, \(\text{ p }{-}\text{ Pb },\) and \(\text{ Pb }{-}\text{ Pb } \) collisions at LHC energies and reported in Refs. [14, 26, 27, 34]. The lifetime of the \({\textrm{K}}^{\mathrm{{* 0}}}\) meson is about 4 fm/c, which is comparable to the lifetime of the hadronic phase. In contrast, the lifetime of the \(\phi \) meson is 10 times higher. The \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons are thus useful probes of the late-stage evolution of high-energy hadronic collisions. As they have similar mass but differ in their strangeness content by one unit, they are suitable candidates to understand the modification of particle production due to rescattering effects and the role of the strangeness content, as discussed in Refs. [26, 28, 29, 34].

This article reports the first measurements of the rapidity dependence of K\(^{*}(892)^{0}~\)and \(\mathrm {\phi (1020)}~\)mesons production in \(\text{ p }{-}\text{ Pb } \) collisions at center-of-mass energy per nucleon–nucleon collisions at \(\sqrt{s_{\textrm{NN}}} = {5.02}\) TeV by the ALICE experiment at the LHC. The large size of the data sample and the excellent particle identification (PID) provide opportunities to extend these measurements in a wider rapidity interval and multiplicity classes compared to earlier measurements [23, 24, 26,27,28,29]. This enables the investigation of the nuclear effects on the particle production in \(\text{ p }{-}\text{ Pb } \) collisions. The \(p_{\textrm{T}}\) spectra, \(Y_{\textrm{asym}}(p_{\textrm{T}})\) and \( {Q}_{\textrm{CP}} (p_{\textrm{T}})\) are studied in the rapidity range \(-1.2< y < 0.3\) and three multiplicity classes along with a measurement on the multiplicity-integrated sample. Similar measurements for charged particles and strange hadrons at RHIC and the LHC energies were reported in Refs. [7, 16, 18].

The measurements presented here are compared with various model predictions such as EPOS-LHC [35], EPOS3 with and without UrQMD [36,37,38,39], DPMJET [40], HIJING [41] and PYTHIA8/Angantyr [42]. EPOS-LHC is an event generator for minimum-bias hadronic interactions that incorporates a parameterization of flow based on LHC data [35]. It is an event generator based on multiple partonic scatterings described using Gribov’s Reggeon field theory formalism, supplemented with collective hadronization and the core-corona mechanism from pp to A–A collisions [35]. EPOS3 is an event generator based on 3\(+\)1D viscous hydrodynamical evolution in the Gribov–Regge multiple scattering framework, which is used to understand hadronic resonances and their interactions in the partonic and hadronic medium. The UrQMD model introduces the description of rescattering and regeneration effects in the hadronic phase [36, 38, 39]. DPMJET is a QCD-inspired dual parton model based on the formalism of the Gribov–Glauber approach that treats the soft- and hard-scattering interactions differently [40]. HIJING is used to study jet and the associated particle production in high energy collisions and is based on QCD-inspired models, including multiple minijet production, soft excitation, nuclear shadowing of parton distribution functions, and jet interaction in dense matter. Due to nuclear shadowing, the parton distribution functions of quarks and gluons are expected to differ from the simple superposition of their distribution in a nucleon. The initial nuclear shadowing effect on particle production is tested using two shadowing depth constant values 0.1 and 0 that represent the HIJING model predictions with and without shadowing [41]. PYTHIA8/Angantyr is an event generator based on the Fritiof model [43], which includes the features of multi-parton interactions and diffractive excitation in each nucleon–nucleon (NN) sub-collision. It acts as a baseline for understanding the non-collective background in the observables sensitive to collective behavior, as PYTHIA8 does not include a collective expansion stage in its description of pp collisions [42]. The comparison of data with the results from these phenomenological models helps to understand the relative contribution of the nuclear effects on the particle production in \(\text{ p }{-}\text{ Pb } \) collisions.

For the results presented here, K\(^{*}(892)^{0}\) and \({{\overline{\textrm{K}}}}^{*}(892)^{0}\) are averaged and denoted by the symbol \({\textrm{K}}^{\mathrm{{* 0}}}\), while \(\mathrm {\phi (1020)}~\)is denoted by \(\phi \). The article is organized as follows. In Sect. 2, the data sample, event and track selection criteria, the analysis techniques, the procedure of extraction of the yields, and the study of the systematic uncertainties are discussed. In Sect. 3, the results on the \(p_{\textrm{T}} \) spectra, the dN/dy, the \(\langle p_{\mathrm{{T}}}\rangle \), the \(Y_{\textrm{asym}}\) and the \({Q}_{\textrm{CP}}\) in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}} = 5.02\) TeV are presented. Finally, the results are summarized in Sect. 4.

2 Data analysis

Measurements of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) meson production are carried out on the data sample collected in 2016 during the second LHC run with \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  \(=\) 5.02 TeV. The resonances are reconstructed from their decay products by using the invariant-mass method. The considered decay channels are \({\textrm{K}}^{\mathrm{{* 0}}}\) \(\rightarrow \textrm{K}^{+} \pi ^{-}\) and its charge conjugate, and \(\phi \) \(\rightarrow \textrm{K}^{+}\textrm{K}^{-}\) with respective branching ratios (BR) of 66.6\(\%\) and 49.2\(\%\) [26, 27]. In the \(\text{ p }{-}\text{ Pb } \) configuration, the \({}^{208}\)Pb beam with energy of 1.58 TeV per nucleon collides with a proton beam with an energy of 4 TeV resulting in collisions at a nucleon–nucleon center-of-mass energy \(\sqrt{s_{\textrm{NN}}}\)  \(=\) 5.02 TeV [26]. It leads to the rapidity in the center-of-mass frame being shifted by \(\Delta y = -0.465\) in the direction of the proton beam with respect to the laboratory frame. The measurements are performed in the rapidity range \(-1.2< y < 0.3\) for five rapidity intervals with width of 0.3 units and three multiplicity classes along with a multiplicity-integrated class. The details of the ALICE detector setup and its performance can be found in Refs. [44, 45]. The measurements are carried out with the ALICE central barrel detectors, which are utilized for tracking, PID, and primary vertex reconstruction and are housed inside a solenoidal magnet with a magnetic field of 0.5 T. The main detectors that are used for the analyses presented here are the Inner Tracking System (ITS) [46], the Time Projection Chamber (TPC) [47], and the TOF (Time-Of-Flight) [48] detectors. These detectors have full azimuthal coverage and have a common pseudorapidity coverage of \(|\eta |< 0.9.\)

2.1 Event and track selection and particle identification

The trigger and event selection criteria are the same as those discussed in previous publications [26, 27]. The events are selected with a minimum-bias trigger based on the coincidence of signals in two arrays of 32 scintillator detectors covering full azimuth and the pseudorapidity regions \(2.8<\eta < 5.1\) (V0A) and \(-3.7<\eta < -1.7\) (V0C) [49]. The primary vertex of the collision is determined using the charged tracks reconstructed in the ITS and the TPC. Events are selected whose reconstructed primary vertex position lies within ± 10 cm from the center of the detector along the beam direction. The Silicon Pixel Detector (SPD) which is the innermost detector of the ITS, is used to reject events in which multiple collision vertices are found (pile-up) [46]. In this work, approximately 540 million events are selected with the criteria described above. The minimum-bias events are further divided into three multiplicity classes, which are expressed in percentiles according to the total charge deposited in the V0A detector [49]. The yield of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons is measured in five rapidity regions \(-1.2< y< -0.9, -0.9< y < -0.6,\) \( -0.6< y < -0.3,\) \( -0.3< y < 0\) and \(0< y < 0.3\) for the multiplicity classes 0–10\(\%\), 10–40\(\%\), 40–100\(\%\) in addition to the multiplicity-integrated (0–100\(\%\)) measurement, corresponding to all minimum-bias events. The 10\(\%\) of the events with the highest multiplicity of charged particles correspond to the 0–10\(\%\) class and similarly, the 40–100\(\%\) class corresponds to the lowest multiplicity. The \(\langle N_{\textrm{coll}} \rangle \) values are estimated from a Glauber model analysis [50] of the charged particle multiplicity distribution in the V0A detector, and they are 13.8 ± 3.8, 10.5 ± 3.9 and 4.0 ± 2.6, respectively for 0–10\(\%\), 10–40\(\%\), and 40–100\(\%\) multiplicity classes taken from Ref. [51]. Charged-particle tracks reconstructed in the TPC with \(p_{\textrm{T}} > 0.15\) GeV/c and pseudorapidity \(|\eta |< 0.8\) are selected for the analysis. The selected charged tracks should have crossed at least 70 out of 159 readout-pad rows of the TPC. The distance of closest approach of the track to the primary vertex in the longitudinal direction (DCA\(_{z}\)) is required to be less than 2 cm. In the transverse plane (xy) a \(p_{\textrm{T}}\)-dependent selection of DCA\(_{xy}(p_{\textrm{T}})< 0.0105 + 0.035\) \(p_{\textrm{T}} ^{-1.1}\) cm is applied. The \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons are reconstructed from their decay daughters (pions and kaons), which are identified by measuring the specific ionization energy loss (dE/dx) in the TPC [47] and their time-of-flight information using the TOF [48]. For the selection of pions and kaons, the measured \(\langle \text {d}E/\text {d}x \rangle \) is required to be within \(n\sigma _{\textrm{TPC}}\) from the expected \(\langle \text {d}E/\text {d}x \rangle \) values for a given mass hypothesis, where \(\sigma _{\textrm{TPC}}\) is the TPC \(\langle \text {d}E/\text {d}x \rangle \) resolution. The values of n are momentum-dependent (p) and are set to 6\(\sigma _{\textrm{TPC}}\), 3\(\sigma _{\textrm{TPC}}\), and 2\(\sigma _{\textrm{TPC}}\) in the momentum intervals \(p < 0.3\,\) GeV/c, \(0.3< p < 0.5\,\) GeV/c and \(p > 0.5\,\) GeV/c, respectively. If the TOF information is available for the considered tracks, it is used for pion and kaon identification in addition to the TPC one by requiring the time-of-flight of the particle to be within 3\(\sigma _{\textrm{TOF}}\) from the expected value for the considered mass hypothesis, where \(\sigma _{\textrm{TOF}}\) is the time-of-flight resolution of the TOF.

2.2 Yield extraction

The \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) resonances are reconstructed from their decay products using the invariant-mass reconstruction technique described in Refs. [26, 27]. The invariant-mass distributions of K\(^{\pm }\pi ^{\mp }\) and K\(^+\)K\(^-\) pairs in the same event are reconstructed. The shape of uncorrelated background is estimated using two techniques, namely mixed-event and like-sign methods. In the mixed-event method, the shape of the uncorrelated-background distribution for \({\textrm{K}}^{\mathrm{{* 0}}}\) (\(\phi \)) is obtained by combining pions (kaons) from a given event with opposite-sign kaons from other events. Each event is mixed with five different events to reduce the statistical uncertainties of the estimated uncorrelated-background distribution. The events which are mixed are required to have similar characteristics, i.e. the longitudinal position of the primary vertices should differ by less than 1 cm, and the multiplicity percentiles, computed from the V0A amplitude, should differ by less than 5\(\%\). The mixed-event distributions for \({\textrm{K}}^{\mathrm{{* 0}}}\) (\(\phi \)) candidates are normalized in the invariant mass interval \(1.1< M_{{\textrm{K}}\pi } < 1.15\,\) GeV/\(c^{2}\) \((1.06< M_{\textrm{KK}} < 1.09\,\) GeV/\(c^{2}),\) which is well separated from signal peak.

Fig. 1
figure 1

Invariant-mass distributions after combinatorial background subtraction for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) candidates in the multiplicity class 0–10\(\%\) and transverse momentum range \(2.2 \le p_{\textrm{T}} < 3.0\,\) GeV/c in the rapidity interval \(-0.3< y < 0\) (panels a and b) and \(0< y < 0.3\) (panels c and d). The \({\textrm{K}}^{\mathrm{{* 0}}}\) peak is described by a Breit–Wigner function whereas the \(\phi \) peak is fitted with a Voigtian function. The residual background is described by a polynomial function of order 2

In the like-sign method, tracks with the same charge from the same event are paired to estimate the uncorrelated background contribution. The invariant-mass distribution for the uncorrelated background is obtained as the geometric mean 2\(\sqrt{n^{++}\times n^{--}}\), where \(n^{++}\) and \(n^{--}\) are the number of positive-positive and negative-negative pairs in each invariant-mass interval, respectively. The mixed-event technique is used as the default method to extract the yields for both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons, while the difference with respect to the yield obtained using the combinatorial background from the like-sign method is included in the estimation of the systematic uncertainty. After the subtraction of the combinatorial background, the invariant-mass distribution consists of a resonance peak sitting on the top of a residual background of correlated pairs. The residual background originates from correlated pairs from jets, misidentification of pions and kaons from \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) meson decays, and partially reconstructed decays of higher-mass particles [26]. Figure 1 shows the \({K^{\pm }\mathrm {\pi ^{\mp } }}\) and \(\mathrm{{K^{+}K^{-}}}\) invariant-mass distributions after subtraction of mixed-event background in the transverse momentum interval \(2.2 \le p_{\textrm{T}} < 3.0\,\) GeV/c for the rapidity intervals \(-0.3< y < 0\) (panels (a) and (b)) and \(0< y < 0.3\) (panels (c) and (d)) in the 0–10\(\%\) multiplicity class.

The signal peak is fitted with a Breit–Wigner and a Voigtian function (convolution of Breit–Wigner and Gaussian functions) for \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) resonances, respectively. For the \({\textrm{K}}^{\mathrm{{* 0}}} \), a pure Breit–Wigner is used because the invariant mass resolution is negligible with respect to the natural width of the resonance peak. A second-order polynomial function is used to describe the shape of the residual background for both resonances. The fit to the invariant-mass distribution is performed in the interval \(0.75< M_{\mathrm{{K}}\pi } < 1.15\,\) GeV/\(c^{2}\) \((0.99< M_{\textrm{KK}} < 1.07\,\) GeV/\(c^{2})\) for \({\textrm{K}}^{\mathrm{{* 0}}} \)(\(\phi \)). The widths of \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) peaks are fixed to their known widths \(\Gamma \)(K*\(^{0})=47.4\pm 0.6\,\text {MeV}/c^{2}\), \(\Gamma (\phi )=4.26 \pm 0.04\,\text {MeV}/c^{2}\) [52], whereas the resolution parameter of the Voigtian function for \(\phi \) is kept as a free parameter. In the estimation of the systematic uncertainties, the width of the Breit–Wigner is taken as a free parameter. The mass and width values extracted from the fit have similar magnitude and trend with \(p_{\textrm{T}}\) as reported in previous publications [26, 27, 34, 53]. In the present study, it is found that the mass and width obtained from the fit are independent of rapidity and multiplicity for both \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) mesons. The sensitivity of the systematic uncertainty to the choice of the fitting range, normalization interval of the mixed-event background, shape of the residual background function, width, and resolution parameters have been studied by varying the fit configuration, as described in Sect. 2.3. The raw yields of \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) mesons are extracted in the transverse momentum range from 0.8 to 16  GeV/c for various rapidity intervals and multiplicity classes.

To obtain the transverse momentum spectra, the raw yields are normalized by the number of accepted non-single-diffractive (NSD) events and corrected for the branching ratio and the detector acceptance (A) times the reconstruction efficiency \((\epsilon _{\textrm{rec}}).\) The \(A \times \epsilon _{\textrm{rec}}\) is obtained from a Monte Carlo (MC) simulation based on the DPMJET [40] event generator and the GEANT3 package to model the transport of the generated particles through the ALICE detector [54]. The \(A \times \epsilon _{\textrm{rec}}\) is defined as the ratio of the reconstructed \(p_{\textrm{T}}\) spectra of \({\textrm{K}}^{\mathrm{{* 0}}} \) \((\phi )\) mesons in a given rapidity interval to the generated ones in the same rapidity interval. The track and PID selection criteria applied to the decay products of resonances in the MC are identical to those used in the data. Since the efficiency depends on \(p_{\textrm{T}}\) and the \(p_{\textrm{T}}\) distributions of \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) mesons from DPMJET are different from the real data, a re-weighting procedure is applied to match the generated \(p_{\textrm{T}}\) shapes to the measured ones.

The effect of the re-weighting on \(A\times \epsilon _{\textrm{rec}}\) depends on \(p_{\textrm{T}} \) and amounts to \(\sim \) 5–17\(\%\) at \(p_{\textrm{T}} \) \(<1.5\) GeV/c. At higher \(p_{\textrm{T}} \), the effect is negligible. The effect of re-weighting also depends on rapidity at low \(p_{\textrm{T}}\). The re-weighted \(A \times \epsilon _{\textrm{rec}}\) is used to correct the raw \(p_{\textrm{T}} \) distribution. The \(A \times \epsilon _{\textrm{rec}}\) is calculated for each rapidity interval and multiplicity class considered in the analysis. The \(A \times \epsilon _{\textrm{rec}}\) as a function of \(p_{\textrm{T}}\) shows a rapidity dependence for a given multiplicity class, however, no significant multiplicity dependence of \(A \times \epsilon _{\textrm{rec}}\) is observed for a given rapidity interval.

2.3 Systematic uncertainties

The procedure to estimate the systematic uncertainties is similar to the one adopted in previous analyses [26, 27]. The sources of systematic uncertainties on the measured yield of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons are signal extraction, track selection criteria, PID, global tracking efficiency, uncertainty on the material budget of the ALICE detector, and the hadronic interaction cross section in the detector material. A summary of the systematic uncertainties on the \(p_{\textrm{T}}\) spectra is given in Table 1. The uncertainty due to signal extraction is estimated from the variation of the yields when varying the invariant mass fit range, the treatment of the Breit–Wigner width in the fits, the mixed-event background normalization interval, the choice of residual background function, and the method to determine the combinatorial background. The fitting range is varied by 50 MeV\(/c^{2}\) for \({\textrm{K}}^{\mathrm{{* 0}}}\) and 10 MeV\(/c^{2}\) for \(\phi \). The normalization interval of the mixed-event background is varied by 150 (50) MeV\(/c^{2}\) with respect to the default value for \({\textrm{K}}^{\mathrm{{* 0}}}\) (\(\phi \)). The width of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) resonances is left as a free parameter in the fit, instead of fixing it to the world-average value. For \(\phi \) resonances, the effect on the yield due to the variation of resolution parameter (\(\sigma \) of the Gaussian function in the Voigtian distribution) is also considered. The residual background is parameterized using first-order and third-order polynomial functions for estimating its contributions to the systematic uncertainties. The combinatorial background from the like-sign method is used instead of the one from event mixing. The estimated systematic uncertainties due to the yield extraction is 5.2\(\%\) for \({\textrm{K}}^{\mathrm{{* 0}}} \) and 3.3\(\%\) for \(\phi \). The systematic effects due to charged track selection have been studied by varying the selection criteria on the number of crossed rows in the TPC, the ratio of the numbers of TPC crossed rows to findable clusters, and the DCA to the primary vertex of the collisions. The estimated uncertainties due to the track selection is 2.5\(\%\) for \({\textrm{K}}^{\mathrm{{* 0}}} \) and about 5\(\%\) for the \(\phi \) mesons. To estimate the systematic uncertainty due to the PID, the selections on the dE/dx and time-of-flight of the pions and kaons are varied. Two momentum-independent selections: 2\(\sigma _{\mathrm{{TPC}}}\) with 3\(\sigma _{\mathrm{{TOF}}}\) and 2\(\sigma _{\mathrm{{TPC}}}\) only, are used for both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \). The estimated systematic uncertainties are 3\(\%\) for \({\textrm{K}}^{\mathrm{{* 0}}} \) and 1.7\(\%\) for \(\phi \). The uncertainty due to the global tracking efficiency, description of the detector material budget in the simulation, and the cross sections for hadronic interactions in the material are taken from Ref. [26]. The total systematic uncertainty is taken as the quadratic sum of all contributions leading to 7.5\(\%\) for \({\textrm{K}}^{\mathrm{{* 0}}} \) and 7.3\(\%\) for \(\phi \) mesons. No multiplicity and rapidity dependence of the systematic uncertainties is observed. Therefore, the systematic uncertainties on the \(p_{\textrm{T}} \) spectra determined for minimum-bias events in the rapidity interval \(0< y < 0.3\) are assigned in all rapidity intervals and multiplicity classes.

Table 1 Relative systematic uncertainties for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) yields in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  \(=\) 5.02 TeV. The quoted relative uncertainties are averaged over \(p_{\textrm{T}} \) in the range 0.8–16 GeV/c. The total systematic uncertainty is the sum in quadrature of the uncertainties due to each source

The systematic uncertainties on \(Y_{\mathrm{{asym}}}\) are estimated by considering the same approaches and variations as for the corrected yields. The systematic uncertainties due to signal extraction and PID are uncorrelated among different rapidity intervals whereas the other sources of systematic uncertainties such as track selections, global tracking uncertainties, material budget and hadronic interactions are correlated and cancel out in the \(Y_{\mathrm{{asym}}}\) ratio. For the uncorrelated sources of uncertainty, the same variations considered for the yields were studied by estimating their effects on the \(Y_{\mathrm{{asym}}}\) ratio. The resulting uncertainty was estimated to be about 2.5\(\%\) \((2\%)\) for \({\textrm{K}}^{\mathrm{{* 0}}}\) \((\phi )\) mesons. No multiplicity and rapidity dependence of the uncertainties is observed for \(Y_{\mathrm{{asym}}}\). Therefore, the systematic uncertainties determined for minimum-bias events are assigned to the ratios in the different rapidity intervals and multiplicity classes. The systematic uncertainties on the ratios (dN/dy)/(dN/d\(y)_{y=0}\) and \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}} \rangle _{y=0}\) as a function of rapidity are calculated in a similar way as for \(Y_{\mathrm{{asym}}}\). The systematic uncertainties on (dN/dy)/(dN/d\(y)_{y=0}\) and \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}} \rangle _{y=0}\) are 2.2\(\%\) (2\(\%\)) and 1.2\(\% \) (1\(\%\)) for \({\textrm{K}}^{\mathrm{{* 0}}}\) (\(\phi )\), respectively.

3 Results and discussion

The rapidity and multiplicity dependence results on the \(p_{\textrm{T}} \) spectra, the dN/dy, the \(\langle p_{\textrm{T}} \rangle \), the \(Y_{\textrm{asym}}\), and the \({Q}_{\textrm{CP}}\) in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}} = 5.02\) TeV are discussed. The measurements are also compared with various model predictions.

3.1 Transverse momentum spectra

Figures 2 and 3 show the \(p_{\textrm{T}}\) spectra of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  \(=\) 5.02 TeV for five rapidity intervals within \(-1.2< y < 0.3\) and for two multiplicity classes 0–10\(\%\) and 40–100\(\%\), respectively. The ratios of the \(p_{\textrm{T}}\) spectra in different rapidity intervals to that in the interval \(0< y < 0.3\) are presented in the bottom panels of Figs. 2 and 3. The measured \(p_{\textrm{T}}\) spectra of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons in the 0–10\(\%\) multiplicity class show a rapidity dependence at low \(p_{\textrm{T}} \) (< 5 GeV/c) indicating that the production of these resonances is higher in the Pb-going direction \((y < 0)\) than in the p-going direction \((y > 0).\) For high \(p_{\textrm{T}}\), no rapidity dependences are observed.

Fig. 2
figure 2

Top panels: The transverse momentum spectra of \({\textrm{K}}^{\mathrm{{* 0}}}\) for five rapidity intervals within \(-1.2< y < 0.3\) and for two multiplicity classes (0–10\(\%\), 40–100\(\%\)) in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 5.02 TeV. The data for different rapidity intervals are scaled for better visibility. Bottom panels: The ratios of the \(p_{\textrm{T}}\) spectra in various rapidity intervals to that in the interval \(0< y < 0.3\) for a given multiplicity class. The statistical and systematic uncertainties are shown as bars and boxes around the data points, respectively

Fig. 3
figure 3

Top panels: The transverse momentum spectra of \(\phi \) for five rapidity intervals within \(-1.2< y < 0.3\) and for two multiplicity classes (0–10\(\%\), 40–100\(\%\)) in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 5.02 TeV. The data for different rapidity intervals are scaled for better visibility. Bottom panels: The ratios of the \(p_{\textrm{T}}\) spectra in various rapidity intervals to that in the interval \(0< y < 0.3\) for a given multiplicity class. The statistical and systematic uncertainties are shown as bars and boxes around the data points, respectively

3.2 Integrated particle yield and mean transverse momentum

The dN/dy and the \(\langle p_{\textrm{T}} \rangle \) are obtained from the transverse momentum spectra in the measured \(p_{\textrm{T}}\) interval and using a fit function to account for the contribution of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons in unmeasured regions. The spectra are fitted with a Lévy–Tsallis function [55] and the fit function is extrapolated to unmeasured regions at low \(p_{\textrm{T}}\) (< 0.8 GeV/c). The integral of the fit function in the extrapolated region accounts for 33\(\%\) (39\(\%\)) of the total yield in the 0–10\(\%\) (40–100\(\%\)) multiplicity class for both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons. The contribution of the extrapolated yield at low \(p_{\textrm{T}}\) is the same for all rapidity intervals. The contribution of the yield in the unmeasured region at high \(p_{\textrm{T}}\) (> 16 GeV/c) is negligible for both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons. The extrapolated yield contribution at low \(p_{\textrm{T}} \) obtained with different fitting functions (i.e., \(m_{\textrm{T}}\)-exponential, Bose–Einstein and Boltzmann–Gibbs Blast-Wave function [56]) and that obtained with the default Lévy–Tsallis function is 5\(\%\) (8\(\%\)) for the 0–10\(\%\) (40–100\(\%\)) included as the systematic uncertainties in the dN/dy and it varies by 2–5\(\%\) for the \(\langle p_{\textrm{T}} \rangle \). In Fig. 4 the dN/dy (top panels) and \(\langle p_{\textrm{T}} \rangle \) (bottom panels) of \({\textrm{K}}^{\mathrm{{* 0}}}\) (left) and \(\phi \) (right) mesons are shown as a function of y for minimum-bias \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}} \) = 5.02 TeV. The central values of the dN/dy of both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons decrease slightly from the rapidity interval \(-1.2< y < -0.9\) to \(0< y < 0.3\) even though within the systematic uncertainties all the data points are compatible among each other. Nevertheless, considering that the systematic uncertainties are mostly correlated among the rapidity intervals, the measured dN/dy values suggest a decreasing trend with increasing y in the rapidity interval covered by the measurement. The \(\langle p_{\textrm{T}} \rangle \) is constant as a function of rapidity for both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) resonances.

Fig. 4
figure 4

The \(p_{\textrm{T}}\) integrated yield (d\(N/{\textrm{d}}y\)) (top panels) and mean transverse momentum (\(\langle p_{\textrm{T}} \rangle \)) (bottom panels) for \({\textrm{K}}^{\mathrm{{* 0}}}\) (left) and \(\phi \) (right) mesons as a function of y measured for the multiplicity class 0–100\(\%\) in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  \(=\) 5.02 TeV. The predictions from EPOS-LHC [35], EPOS3 with and without UrQMD [38, 39], DPMJET [40], HIJING [41], and PYTHIA8/Angantyr [42] are shown as different curves. The statistical uncertainties are represented as bars whereas the boxes indicate total systematic uncertainties

The model predictions from EPOS-LHC [35], EPOS3 with and without UrQMD [38, 39], DPMJET [40], HIJING [41], and PYTHIA8/Angantyr [42] are also shown in the Fig. 4. In general, the models show a similar trend with rapidity as the data except EPOS3 with and without UrQMD for \(\langle p_{\textrm{T}} \rangle \), which shows a pronounced decreasing trend with rapidity. All the model predictions shown in Fig. 4 underestimate the \(\langle p_{\textrm{T}} \rangle \) of both meson species. For the dN/dy, HIJING and EPOS3 with and without UrQMD overpredict the measured values for both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \), while PYTHIA8/Angantyr overpredicts the \({\textrm{K}}^{\mathrm{{* 0}}}\) and underpredicts the \(\phi \) yield. EPOS-LHC provides the best overall description of the dN/dy and \(\langle p_{\textrm{T}} \rangle \) measurements for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons. The \(\langle p_{\textrm{T}} \rangle \) also shows a flat behavior as a function of rapidity for all the considered multiplicity classes as it can be seen in Fig. 10 of Appendix A.

A similar behavior in the average transverse kinetic energy as a function of rapidity for strange hadrons was reported in Ref. [16]. The rapidity dependence of dN/dy and \(\langle p_{\textrm{T}} \rangle \) for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons in the multiplicity class 0–100\(\%\) is further studied by dividing the dN/dy and \(\langle p_{\textrm{T}} \rangle \) values in a given rapidity interval by the corresponding values at \(y=0\), as shown in Fig. 5. The dN/dy and \(\langle p_{\textrm{T}} \rangle \) value at \(y=0\) is computed from the \(p_{\textrm{T}} \) spectrum measured in the rapidity interval \(-0.3< y < 0.3.\) The systematic uncertainties on these ratios are estimated by studying the effects of the variations directly on the ratios as discussed in Sect. 2.3. This procedure takes into account the correlation of the systematic uncertainties across rapidity bins: as a result, these ratios have smaller systematic uncertainties than those on the dN/dy and \(\langle p_{\textrm{T}} \rangle \), and allow for a better insight into the y dependence. The ratio (dN/dy)/(dN/d\(y)_{y=0}\) decreases with rapidity, whereas \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}}\rangle _{y=0}\) shows a flat behavior within uncertainties as a function of rapidity for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons. The measurements are compared with various model predictions. The predictions from HIJING qualitatively reproduce the trend and are the closest to the data for both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \). The predictions from PYTHIA8/Angantyr, DPMJET, EPOS-LHC, EPOS3 with and without UrQMD show a decreasing trend of (dN/dy)/(dN/d\(y)_{y=0}\) with increasing y, but the rapidity dependence is less pronounced than the one in data, as it can be seen by the fact that they all tend to underestimate the measured yield ratios in the lowest rapidity intervals, especially for \({\textrm{K}}^{\mathrm{{* 0}}}\) meson.

For \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}}\rangle _{y=0}\) as a function of y, also shown in Fig. 5, EPOS3 with and without UrQMD overestimate the measurements at low y and predict a marked decreasing trend of \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}}\rangle _{y=0}\) with rapidity, which is not supported by the data. From the other models, less pronounced trends are expected, which are consistent with the data. In particular, HIJING predicts a slightly decreasing \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}}\rangle _{y=0}\) with increasing rapidity, while PYTHIA8/Angantyr, DPMJET, and EPOS-LHC predict a slightly increasing trend.

Similar studies of the ratio \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}}\rangle _{y=0}\) of charged hadrons in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 5.02 TeV compared with the predictions of hydrodynamics and color-glass condensate (CGC) model were reported in [33]. Predictions from hydrodynamic calculations show a decrease in \(\langle p_{\textrm{T}} \rangle \) with rapidity, whereas CGC predicts an increase in \(\langle p_{\textrm{T}} \rangle \) with rapidity [33], while the data are flat within uncertainties. The dN/dy and \(\langle p_{\textrm{T}} \rangle \) increase with multiplicity at midrapidity as observed for light-flavor hadrons and resonances in pp and \(\text{ p }{-}\text{ Pb } \) collisions [26, 27, 56]. A similar behavior is observed in this article for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) in all the different rapidity intervals shown in Fig. 10 in Appendix A.

Fig. 5
figure 5

The \(p_{\textrm{T}}\) integrated yield (d\(N/{\textrm{d}}y\)) (upper panels) and mean transverse momentum \(\left( \langle p_{\textrm{T}} \rangle \right) \) (bottom panels) for \({\textrm{K}}^{\mathrm{{* 0}}}\) (left) and \(\phi \) (right) mesons as a function of y, divided by the d\(N/{\textrm{d}}y\) and \(\langle p_{\textrm{T}} \rangle \) at \(y = 0\) for the multiplicity class 0–100\(\%\) in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  \(=\) 5.02 TeV. The predictions from EPOS-LHC [35], EPOS3 with and without UrQMD [38, 39], DPMJET [40], HIJING [41], and PYTHIA8/Angantyr [42] are shown as different curves. The statistical uncertainties are represented as bars whereas the boxes indicate total systematic uncertainties

Fig. 6
figure 6

Rapidity asymmetry (\(Y_{\textrm{asym}}\)) of K\(^{*0}\) (red circles) and \(\phi \) (blue squares) meson production as a function of \(p_{\textrm{T}}\) in the rapidity range \(0< |y| < 0.3\) for various multiplicity classes in p–Pb collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 5.02 TeV. The statistical uncertainties are shown as bars whereas the boxes represent the systematic uncertainties on the measurements

Fig. 7
figure 7

The comparison of experimental results of \(Y_{\textrm{asym}}\) for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) meson production as a function of \(p_{\textrm{T}}\) in the rapidity range \(0< |y| < 0.3\) with the model predictions from EPOS-LHC [35], EPOS3 with and without UrQMD [38, 39], DPMJET [40], HIJING [41], and PYTHIA8/Angantyr [42]. Data points are shown with blue markers, and model predictions are shown by different color bands, where bands represent the statistical uncertainity of the model. The statistical uncertainties on the data points are represented as bars whereas the boxes indicate total systematic uncertainties

Fig. 8
figure 8

The \({Q}_{\textrm{CP}}\) of \({\textrm{K}}^{\mathrm{{* 0}}}\) (red circles) and \(\phi \) (blue squares) mesons as a function of \(p_{\textrm{T}} \) for 0–10\(\%\)/40–100\(\%\) (top panels) and 10–40\(\%\)/40–100\(\%\) (bottom panels) in various rapidity intervals within the range \(-1.2< y < 0.3\) in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 5.02 TeV. The statistical and systematic uncertainties are represented by vertical bars and boxes, respectively

3.3 Rapidity asymmetry

The rapidity asymmetry \((Y_{\textrm{asym}})\) is calculated from \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons yields in \(-0.3< y < 0\) and 0\(< y < 0.3,\) as defined by Eq. 1. Figure 6 shows the \(Y_{\textrm{asym}}\) of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons in the measured \(p_{\textrm{T}}\) intervals for various multiplicity classes in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 5.02 TeV.

The \(Y_{\textrm{asym}}\) values for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) as a function of \(p_{\textrm{T}}\) are consistent within uncertainties for all multiplicity classes. The \(Y_{\textrm{asym}}\) values deviate from unity at low \(p_{\textrm{T}} \) (< 5 GeV/c), suggesting the presence of a rapidity dependence in the nuclear effects. The deviations are more significant for events with high multiplicity. The \(Y_{\textrm{asym}}\) values are consistent with unity at high \(p_{\textrm{T}} \) \((> 5\,\) GeV/c) for all multiplicity classes, suggesting the absence of nuclear effects at high \(p_{\textrm{T}} \) for the production of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons in \(\text{ p }{-}\text{ Pb } \) collisions. Similar results have been reported for charged hadrons, pions, protons in d–Au collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 200  GeV by the STAR Collaboration [18] and for charged hadrons and multi-strange hadrons in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 5.02 TeV by the CMS Collaboration as discussed in Refs. [16, 17]. Figure 7 shows the comparison of the measured \(Y_{\textrm{asym}}\) for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons as a function of \(p_{\textrm{T}}\) in minimum-bias events (0–100\(\%\)) with the model predictions from EPOS-LHC, HIJING with and without shadowing, DPMJET, PYTHIA8/Angantyr, and EPOS3 with and without UrQMD.

Fig. 9
figure 9

The \(Q_{\textrm{CP}}\) of \({\textrm{K}}^{\mathrm{{* 0}}}\) (red circles) and \(\phi \) (blue squares) mesons as a function of rapidity for 0–10\(\%\)/40–100\(\%\) (solid markers) and 10–40\(\%\)/40–100\(\%\) (open markers) in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  \(=\) 5.02 TeV. The solid and doted lines represents the linear fit to data. The statistical and systematic uncertainties are represented by vertical bars and boxes on the measurements, respectively

HIJING with and without shadowing, and EPOS3 with and without UrQMD describe the measured \(Y_{\textrm{asym}}\) at low \(p_{\textrm{T}}\) within uncertainties, but they significantly overestimate the data at high \(p_{\textrm{T}}\), predicting an increasing trend with \(p_{\textrm{T}}\) (more pronounced for \({\textrm{K}}^{\mathrm{{* 0}}}\) than for \(\phi \)) that is not supported by the measurements, which are consistent with a flat or decreasing trend for both meson species. Model predictions from EPOS-LHC, PYTHIA8/Angantyr, and DPMJET for \({\textrm{K}}^{\mathrm{{* 0}}}\) and DPMJET for \(\phi \) at high \(p_{\textrm{T}}\) are in agreement with the data within uncertainties.

3.4 Nuclear modification factor

The nuclear modification factor \({Q}_{\textrm{CP}}\) is calculated from the \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) yields normalized to \(\langle N_{\textrm{coll}} \rangle \) in high multiplicity (central) and low multiplicity (peripheral) collisions, as defined by Eq. 2. Figure  8 shows the \({Q}_{\textrm{CP}}\) of \({\textrm{K}}^{\mathrm{{* 0}}}\) (red circles) and \(\phi \) (blue squares) mesons as a function of \(p_{\textrm{T}} \) for 0–10\(\%\)/40–100\(\%\) (top panels) and 10–40\(\%\)/40–100\(\%\) (bottom panels) in various rapidity intervals within the range \(-1.2< y < 0.3\) for \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}}\)  = 5.02 TeV. The \({Q}_{\textrm{CP}}\) of \(\phi \) mesons seems to be slightly higher than the \({\textrm{K}}^{\mathrm{{* 0}}}\) one for the ratio of 0–10\(\%\)/40–100\(\%\), however, the results for the two meson species are consistent within uncertainties for the ratio 10–40\(\%\)/40–100\(\%\) for all measured rapidity intervals. An enhancement at intermediate \(p_{\textrm{T}}\) \((2.2< p_{\textrm{T}} < 5.0\,\) GeV/c),  reminiscent of the Cronin effect, is seen for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons in the \({Q}_{\textrm{CP}}\). This enhancement is more pronounced at high negative rapidity, i.e., in the Pb-going direction, and for high multiplicity events. The more pronounced Cronin-like enhancement for the 0–10\(\%\) multiplicity class suggests that multiple scattering effects are more relevant for high multiplicity (central) collisions. At high \(p_{\textrm{T}}\) \((> 5\,\) GeV/c),  the \(Q_{\textrm{CP}}\) values are greater than unity, which is a known feature of \(Q_{\textrm{pPb}}\)Footnote 1 and \(Q_{\textrm{CP}}\) when the centrality or multiplicity classes are defined with the V0 detector, and it is interpreted as a selection bias due to the multiplicity estimator [51]. The results for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons are consistent between each other within uncertainties.

To quantify the rapidity dependence of the nuclear modification factor, the \({Q}_{\textrm{CP}}\) values of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons for intermediate \(p_{\textrm{T}} \) \((2.2< p_{\textrm{T}} < 5.0\) GeV/c) are shown as a function of rapidity in Fig. 9. The values of \({Q}_{\textrm{CP}}\) at intermediate \(p_{\textrm{T}} \) show a faster decrease from the rapidity interval \(-1.2< y < - 0.9\) to \(0< y < 0.3\) for 0–10\(\%\)/40–100\(\%\) than for 10–40\(\%\)/40–100\(\%\), indicating a stronger rapidity dependence of the Cronin-like enhancement in events with high multiplicity. The stronger rapidity dependence for 0–10\(\%\)/40–100\(\%\) can be inferred from the slope parameter \((\alpha )\) of the linear function fit to the \({Q}_{\textrm{CP}}\) of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons reported in Fig. 9. The slope of the \(\phi \) meson \({Q}_{\textrm{CP}}\) is slightly larger than the \({\textrm{K}}^{\mathrm{{* 0}}}\) one. A similar conclusion on the \(\eta \) dependence of nuclear modification factors of charged hadrons was reported by the BRAHMS Collaboration [7].

4 Summary

The transverse momentum differential yields of \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) mesons have been measured in the rapidity interval \(-1.2< y < 0.3\) for various multiplicity classes over the transverse momentum range \(0.8< p_{\textrm{T}} <16\,\) GeV/c in \(\text{ p }{-}\text{ Pb } \) collisions at \(\sqrt{s_{\textrm{NN}}} = 5.02\) TeV with the ALICE detector. The \(p_{\textrm{T}}\) spectra of \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) mesons show a multiplicity and rapidity dependence at low \(p_{\textrm{T}}\), whereas the spectral shapes are similar for all multiplicity classes and rapidity intervals at high \(p_{\textrm{T}}\) (> 5  GeV/c). This suggests that nuclear effects influence \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) meson production at low \(p_{\textrm{T}}\). The (dN/dy)/(dN/d\(y)_{y=0}\) ratios decreases with increasing rapidity in the measured interval \(-1.2< y < 0.3,\) whereas the average transverse momentum \(\left( \langle p_{\textrm{T}} \rangle \right) \) and the \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}} \rangle _{y=0}\) ratios show a flat behavior for both \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons. The rapidity dependence of \({\textrm{d}}N/{\textrm{d}}y\), \(\langle p_{\textrm{T}} \rangle \) and their ratios with respect to the corresponding values at \(y =0\) are compared with model predictions for minimum-bias events. The EPOS-LHC model, which includes parameterized flow, provides the best description for the magnitudes of \({\textrm{K}}^{\mathrm{{* 0}}} \) and \(\phi \) \({\textrm{d}}N/{\textrm{d}}y\) and \(\langle p_{\textrm{T}} \rangle \), whereas HIJING predictions are in closest agreement with the measured rapidity dependence, which is studied via the ratios (dN/dy)/(dN/d\(y)_{y=0}\) and \(\langle p_{\textrm{T}} \rangle /\langle p_{\textrm{T}} \rangle _{y=0}\). The \(Y_{\textrm{asym}}\) ratios for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons as a function of \(p_{\textrm{T}}\) show deviations from unity at low \(p_{\textrm{T}}\) for high multiplicity events, while, their values are consistent with unity within uncertainties at high \(p_{\textrm{T}}\) in the measured multiplicity and rapidity intervals. The \(Y_{\textrm{asym}}\) ratios of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons are found to be consistent between each other within uncertainties in the measured kinematic region. The measured deviations of \(Y_{\textrm{asym}}\) from unity at low \(p_{\textrm{T}}\) suggest the presence of rapidity dependent nuclear effects such as multiple scattering, nuclear shadowing, parton saturation, and energy loss in cold nuclear matter. None of the models presented here is able to describe the \(Y_{\textrm{asym}}\) of \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons at low \(p_{\textrm{T}}\). The nuclear modification factors between the central and peripheral collisions \({Q}_{\textrm{CP}}\) for \({\textrm{K}}^{\mathrm{{* 0}}}\) and \(\phi \) mesons as a function of \(p_{\textrm{T}} \) show a bump, with a maximum around \(p_{\textrm{T}}\)  = 3 GeV/c, suggestive of the Cronin effect. This Cronin-like enhancement is more pronounced for large negative rapidities (in the Pb-going direction) and for more central (higher multiplicity) collisions. The measurements reported in this paper confirm that nuclear effects play an important role in particle production in p–Pb collisions at the LHC energies. They will contribute, along with previous and upcoming measurements of other hadron species, to constrain models and event generators.