High-energy heavy ion collisions are useful to study the properties of the quark-gluon plasma (QGP), a strongly coupled medium of deconfined quarks and gluons. It has long been proposed that the yields of quarkonium states are suppressed because of interactions in the QGP [1][2][3][4][5][6][7]. These in-medium effects have been theoretically studied with calculations based on lattice quantum chromodynamics and effective field theories [8][9][10][11]. One of the fundamental aspects of these interactions with the medium is that the amount of suppression for different quarkonium states is expected to be stronger for those with smaller binding energies. On the other hand, quarkonia can also be produced by recombination processes [7][8][9][10][12][13][14][15]. Studies of ϒ production are particularly interesting, because the number of quark-antiquark pairs in a single lead-lead (Pb-Pb) collision is much smaller for bottom than for charm quarks, so that the recombination of independently produced quarks and antiquarks can be neglected [16].
Experimentally, the dynamics of quarkonium production in heavy ion collisions are commonly studied using the nuclear modification factor (R AA ), which is defined as the ratio of particle yields in nucleus-nucleus (AA) collisions to those in proton-proton (pp) collisions scaled by the average number of binary nucleon-nucleon (NN) interactions in the AA events.
Measurements of R AA have been performed at the SPS, RHIC, and LHC accelerators, in both the charmonium and bottomonium families [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. For ϒ states, R AA results have been reported by CMS [27], ALICE [28], and ATLAS [29] at the LHC and by STAR [30,31] and PHENIX [32] at RHIC. These results show a significant suppression of ϒð1SÞ mesons in heavy ion collisions, and ϒð2SÞ mesons are even more suppressed. For the ϒð3SÞ meson, only upper limits have been reported in AA collisions, by CMS [27].
In this Letter, the first observation of the ϒð3SÞ meson in Pb-Pb collisions is reported. Nuclear modification factors for ϒð2SÞ and ϒð3SÞ mesons are found using Pb-Pb and pp data collected with the CMS detector at a nucleonnucleon center-of-mass energy of ffiffiffiffiffiffiffi ffi s NN p ¼ 5.02 TeV in 2018 and 2017, respectively. The double ratios, obtained by dividing the ϒð3SÞ over ϒð2SÞ yield ratios in Pb-Pb by those for pp collisions, are also reported. The ϒ states are identified using their decay into two oppositely charged muons. The results are presented as functions of Pb-Pb collision centrality, as well as ϒ meson transverse momentum (p T ) in the rapidity range of jyj < 2.4. Centrality is related to the overlap of the two lead nuclei and is quantified as the percentage of the total inelastic nucleus-nucleus hadronic cross section, with 0% representing the largest overlap [33]. Tabulated results are provided in the HEPData record for this analysis [34].
The CMS apparatus [35] is a multipurpose, nearly hermetic detector, designed to trigger on [36,37] and identify electrons, muons, photons, and hadrons [38][39][40][41]. A superconducting solenoid of 6 m internal diameter provides a magnetic field of 3.8 T. Within the solenoid volume are the silicon pixel and strip tracker, a crystal electromagnetic calorimeter, and a brass-scintillator hadron calorimeter. The pseudorapidity (η) coverage of these calorimeters is extended by the forward hadron (HF) calorimeters, located at 3 < jηj < 5. Muons are measured in the range jηj < 2.4 using gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid. Events are filtered using a two-tiered trigger system [37]. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors [36]. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software. Centrality is determined using the sum of the total transverse energy deposited in both of the HF calorimeters [42].
The events are selected online using an L1 trigger requiring two muons in a single bunch crossing without explicit requirements on the muon momentum. Additional criteria on the single-muon quality and dimuon mass selection are applied at the HLT in Pb-Pb collisions [43]. The collected pp (Pb-Pb) sample corresponds to an integrated luminosity of 300 pb -1 (1.61 nb -1 ) [44][45][46].
In order to reject beam-gas interactions and nonhadronic collisions, an offline event selection is applied for pp and Pb-Pb collisions [47]. For both systems, events are required to have at least one reconstructed primary vertex (PV), which is the vertex corresponding to the hardest scattering in the event, evaluated using tracking information alone [48]. In addition, more than 25% of the tracks have to pass a tight track-quality selection in pp collisions [41,47]. For Pb-Pb collisions, selections are applied on the silicon pixel detector cluster widths [47,49] and on the HF information [43]. The reconstructed muons are selected using criteria that optimize the muon identification [43]. In addition, to ensure high reconstruction efficiency for ϒ mesons, the individual muons are required to have p T > 3.5 GeV=c and jηj < 2.4.
Monte Carlo (MC) samples for each ϒ meson are simulated using the "CP5" tune [50] of PYTHIA8.212 [51], with the assumption that they are produced unpolarized [52-55]. To reproduce the background environment in Pb-Pb collisions, each PYTHIA ϒ event is embedded into a Pb-Pb event simulated using HYDJET1.9 [56]. The MC events are then weighted to match the measured ϒ p T spectra for either pp or Pb-Pb collisions [43]. A full simulation of the CMS detector using GEANT4 [57] is performed to determine the acceptance and reconstruction efficiency. The feed-down contributions, i.e., decays from heavier quarkonium states, are not considered in this analysis. The effect of such contributions on the kinematic distributions of the simulated ϒ states is, to a large extent, accounted for by the weighting procedure.
The dimuon invariant mass distribution is studied in the 8-14 GeV=c 2 region. To reduce the large amount of background in Pb-Pb collisions, signal-enriched dimuon candidates are selected using boosted decision trees (BDTs), using the TMVA package [58]. For the BDT training, MC samples and dimuons from the invariant mass spectra in data except where the ϒ signals are present (i.e., 8.8-10.8 GeV=c 2 ) are used for the signal and background, respectively. In both cases, the dimuons must satisfy the selection criteria previously mentioned. The training variables for the BDT include the χ 2 probability of the dimuon vertex fit; the distance of closest approach of the ϒ meson momentum vector relative to the PV; the distance between the PV and the dimuon vertex, both in three dimensions and projected onto the transverse plane; and the pointing angle, defined as the angle between the line segment connecting the PV and decay vertex and the momentum vector of the reconstructed particle candidates, again in three dimensions and projected onto the transverse plane. To avoid potential biases, the algorithms and selections are optimized using a quarter of the Pb-Pb data for dimuon p T < 30 GeV=c, jyj < 2.4, and centrality 0%-90%. A threshold is set on the resulting BDT variable to maximize the signal significance, S= ffiffiffiffiffiffiffiffiffiffiffiffi S þ B p , where S and B represent the yields for signal and background, respectively.
The yields of the ϒ states are extracted by an extended unbinned maximum likelihood fit to the dimuon invariant mass distributions. The line shape of each ϒ state is modeled by a sum of three Crystal Ball (CB) functions [59], with common means and tail parameters but independent widths. The means and widths for the ϒð2SÞ and ϒð3SÞ are found by multiplying the fitted values of those for the ϒð1SÞ by the ratio of their world-average masses [60] over that for the ϒð1SÞ. All other fit parameters, except the yields for the excited states, are fixed to those found for the ground state. The shape parameters of the signal fit model, with the exception of the ϒð1SÞ mean, are constrained by Gaussian probability density functions (PDFs), of means and widths determined, for each p T bin, from the central values and uncertainties obtained from fits to the MC samples. The background of the mass distribution is modeled using three different functional forms: an error function multiplied by an exponential, a simple exponential, and a Chebyshev polynomial. An Akaike information criterion (AIC) test [61], which estimates the relative quality of each function, is performed to determine, for each kinematic region, the nominal function for the background PDF. The order of the Chebyshev polynomial is chosen based on a log likelihood ratio test [62].
Figure 1 shows the dimuon invariant mass distribution in Pb-Pb collisions, integrated over the full kinematic range, p T < 30 GeV=c and jyj < 2.4, and centrality 0%-90%.
Acceptance and efficiency correction factors are applied to the extracted number of ϒ mesons to compute the R AA values. The acceptance is computed with simulated samples, as the fraction of generated ϒ mesons that decay to muons of p μ T > 3.5 GeV=c and jη μ j < 2.4. The efficiency is also evaluated with simulated events, as the fraction of accepted dimuons that are reconstructed and pass the trigger and analysis selection criteria. To take into account possible differences between data and MC simulations, the individual components of the dimuon efficiency (reconstruction, identification, and triggering) are measured using single muons from J=ψ meson decays in both data and simulation, with the tag-and-probe (T&P) method [63]. The data over simulation ratios, for each of the three components, are applied as event-by-event weights to each dimuon. The acceptance and efficiency values for the full kinematic region are 40% (43%) and 37% (38%) for ϒð2SÞ [ϒð3SÞ] in Pb-Pb collisions.
The systematic uncertainties are analyzed for various sources and summarized in Table I. For each source, the difference of the signal yields in the variation compared to the nominal is taken as the systematic uncertainty. Three sources are considered for the uncertainty in the signal extraction: choice of signal PDF parameters, choice of signal PDF, and choice of background PDF. For the parametrization of the signal PDF, the fit results using either the ϒð2SÞ or ϒð3SÞ MC samples are used to determine different values for the mean and width of the Gaussian function for each signal parameter. The signal PDF is modified from a sum of three CB functions to a sum of two CB functions and a Gaussian function. For the background PDF systematic study, functions rejected by the AIC test in the nominal background PDF determination, for each kinematic region, are used as alternatives.
The acceptance and efficiency uncertainties are evaluated by varying the p T weights (used to match the simulated and measured ϒ p T spectra) by their fit uncertainties. The uncertainty on the T&P correction of the single muon efficiencies is propagated to the dimuon efficiencies.
For the Pb-Pb analysis, additional uncertainties arise from the BDT training and centrality calibration. To evaluate the uncertainty reflecting the BDT training, the event samples used as signal or background are split in two. Each trained algorithm is applied to the other half of the samples, and the average of the two BDT variable values is used as the nominal. Alternatively, the two individual values of each BDT variable are used, and the largest difference of the corrected signal yields compared to the nominal is taken as the systematic uncertainty. The centrality calibration is varied by changing the boundaries of the centrality intervals [43] considering the inefficiency in the event selection.
The total uncertainty, dominated by the uncertainties from the BDT training and background PDF, is the quadratic sum of the uncertainties from all the different sources. In addition, global uncertainties reflect the integrated luminosity of the pp dataset (1.9%) [44] and the number of minimum bias Pb-Pb collision events (N MB ) sampled by the trigger (1.3%) [64]. The overlap function
TABLE I. Systematic uncertainties from various sources in pp and Pb-Pb collisions listed in percentage. The global uncertainties described in the text are not included in the total uncertainties. ϒð2SÞ (%) ϒð3SÞ (%) ϒð3SÞ=ϒð2SÞ (%) Source pp Pb-Pb pp Pb-Pb pp Pb-Pb BDT selection Á Á Á 0.3-9.0 Á Á Á 1.5-18.6 Á Á Á 1.2-22.8 Background PDF 0.1-1.4 0.3-11.7 0.2-1.6 1.4-21.4 <0.5 0.6-17.6 Signal PDF 0.1-1.1 0.5-2.6 0.4-1.1 0.1-2.5 0.3-0.6 0.1-3.0 Signal parameter 0.1-1.2 0.0-3.8 0.1-1.6 0.3-3.7 0.05-1.4 0.1-0.9 Event selection Á Á Á 0.0-0.5 Á Á Á 0.2-13.1 Á Á Á 0.1-13.6 Correction factors <0.1 <0.5 <0.1 <0.4 <2.0 T&P 0.9-1.0 3.8-4.5 0.9-1.1 3.8-4.4 Á Á Á Total uncertainty 1.0-1.8 3.9-13.5 1.1-2.2 6.0-22.2 0.4-1.5 4.1-23.8 8 9 1 0 1 1 1 2 1 3 1 4 T AA for each centrality interval is the average number of binary NN collisions per Pb-Pb interaction divided by the inelastic NN cross section. Its uncertainty is estimated by varying the Glauber model parameters within their uncertainties [65] and is found to be in the range of 1.8%-5.4%.
The significance of the ϒð3SÞ in Pb-Pb collisions corresponds to 5.6 standard deviations, calculated using the fit likelihood ratio [66-68]. The nuclear modification factors are computed as
where N AA and σ pp are the efficiency-and acceptancecorrected normalized yields in Pb-Pb collisions and the pp cross sections, respectively, for ϒ mesons in a given kinematic range. The average value of T AA computed in each centrality bin is denoted by hT AA i.
Figure 2 presents the R AA values for ϒð2SÞ and ϒð3SÞ, together with the previous measurements for the ϒð1SÞ [27], as functions of ϒ meson p T and Pb-Pb collision centrality, with the latter represented using hN part i, the average number of participating nucleons in each of the centrality intervals [65]. The measurements as a function of p T are reported in five intervals for the ϒð2SÞ, 0-3, 3-6, 6-9, 9-15, and 15-30 GeV=c, and four for the ϒð3SÞ, 0-4, 4-9, 9-15, and 15-30 GeV=c. The centrality intervals and the corresponding hN part i values used for this analysis are listed in the HEPData record [34]. The results for the ϒð1SÞ are taken from Ref. [27], because the new data do not improve the statistical significance. Furthermore, the systematic uncertainties, which are dominant for the ϒð1SÞ results, are found to not improve when including the BDT training.
A gradual decrease of R AA is observed toward more central collisions (i.e., higher N part ) for both ϒð2SÞ and ϒð3SÞ. On the other hand, no significant dependence on p T is found for either of them. Both states are strongly suppressed in central Pb-Pb collisions, as well as over the entire p T region when averaged over centrality. Furthermore, the R AA of the ϒð3SÞ is smaller than that of the ϒð2SÞ, with values integrated over p T and centrality 0%-90% of 0.115 AE 0.008ðstatÞ AE 0.007ðsystÞ and 0.080 AE 0.014ðstatÞ AE 0.012ðsystÞ for the ϒð2SÞ and ϒð3SÞ, respectively. These results indicate, much more clearly than the previous measurements [27-29], that the sequential suppression pattern of the bottomonium states follows the ordering of their binding energies [69].
Figure 3 shows, as functions of hN part i and p T , the ϒð3SÞ to ϒð2SÞ double ratios, obtained by dividing the ϒð3SÞ over ϒð2SÞ yield ratios in Pb-Pb collisions by the corresponding ratios in pp collisions. The bins are identical to those used in the ϒð3SÞ measurement. The uncertainties in the ϒð3SÞ=ϒð2SÞ ratios are computed by propagating the ϒð2SÞ and ϒð3SÞ uncertainties, taking into account their correlation, whereas the global uncertainties on T AA , pp luminosity, and Pb-Pb N MB cancel out. The decreasing trend of the double ratios toward more central Pb-Pb collisions indicates a stronger suppression for the ϒð3SÞ than for the ϒð2SÞ, whereas no significant dependence on p T is seen. Comparisons to various theoretical calculations [9,10,70-73] and previous experimental results [24,[74][75][76][77] can be found in Supplemental Material [78].
In summary, data from Pb-Pb and pp collisions at a nucleon-nucleon center-of-mass energy of ffiffiffiffiffiffiffi ffi s NN p ¼5.02 TeV, collected with the CMS detector, were analyzed to measure the yields and nuclear modification factors, R AA , of the ϒð2SÞ and ϒð3SÞ mesons. The ϒð3SÞ meson is observed for the first time in Pb-Pb collisions, with a significance above 5 standard deviations. Dividing the ϒð3SÞ over ϒð2SÞ yield ratios in Pb-Pb by those in pp collisions gives the double ratios that quantify the relative modification of the two mesons. Results are shown as functions of ϒ transverse momentum and Pb-Pb collision centrality. Both the ϒð2SÞ and ϒð3SÞ mesons are suppressed (R AA < 1), with a stronger effect for the ϒð3SÞ. The suppression increases for more central Pb-Pb collisions, whereas no significant dependence on p T is seen. The ϒð3SÞ over ϒð2SÞ double ratios show a gradual decrease toward more central Pb-Pb collisions, indicating that the degree to which the suppression is stronger for the ϒð3SÞ meson increases for more central Pb-Pb collisions. Combined with previous measurements, these results indicate that the strength of the suppression increases in the sequence ϒð1SÞ, ϒð2SÞ, and ϒð3SÞ. These results provide new constraints on the understanding of the dynamics of quarkonium states in the QGP created in heavy ion collisions. We congratulate our colleagues in the CERN accelerator departments for the excellent performance of the LHC and thank the technical and administrative staffs at CERN and at other CMS institutes for their contributions to the success of the CMS effort. In addition, we gratefully acknowledge the computing centers and personnel of the Worldwide LHC Computing Grid and other centers for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies: BMBWF and FWF (Austria); FNRS and FWO (Belgium); CNPq, CAPES, FAPERJ, FAPERGS, and FAPESP (Brazil); MES and BNSF (Bulgaria); CERN; CAS, MoST, and NSFC (China); MINCIENCIAS (Colombia); MSES and CSF (Croatia); RIF (Cyprus); SENESCYT (Ecuador); MoER, ERC PUT, and ERDF (Estonia); Academy of Finland, MEC, and HIP (Finland); CEA and CNRS/IN2P3 (France); BMBF, DFG, and HGF (Germany); GSRI (Greece); NKFIH (Hungary); DAE and DST (India); IPM (Iran); SFI (Ireland); INFN (Italy); MSIP and NRF (Republic of Korea); MES (Latvia); LAS (Lithuania); MOE and UM (Malaysia); BUAP, CINVESTAV, CONACYT, LNS, SEP, and UASLP-FAI (Mexico); MOS (Montenegro); MBIE (New Zealand); PAEC (Pakistan); MES and NSC (Poland); FCT (Portugal); MESTD (Serbia); MCIN/AEI and PCTI (Spain); MOSTR (Sri Lanka); Swiss Funding Agencies (Switzerland); MST (Taipei); MHESI and NSTDA (Thailand); TUBITAK and TENMAK (Turkey); NASU (Ukraine); STFC (United Kingdom); and DOE and NSF (USA).
FIG. 3. The double ratios of ϒð3SÞ=ϒð2SÞ as functions of hN part i (upper), showing also the 0%-90% centrality interval, and of p T (lower). The vertical lines and boxes correspond to statistical and systematic uncertainties, respectively. The box at unity in the left plot shows the combined systematic and statistical uncertainties from pp data, which is common to all the points.
5.02 TeV Pb þ Pb collisions with the ATLAS experiment, Eur. Phys. J. C 78, 762 (2018). [26] CMS Collaboration, Measurement of prompt and nonprompt charmonium suppression in PbPb collisions at ffiffiffiffiffiffiffi ffi s NN p ¼ 5.02 TeV, Eur. Phys. J. C 78, 509 (2018). [27] CMS Collaboration, Measurement of nuclear modification factors of ϒð1SÞ, ϒð2SÞ, and ϒð3SÞ mesons in PbPb collisions at ffiffiffiffiffiffiffiffi s NN p ¼ 5.02 TeV, Phys. Lett. B 790, 270 (2019). [28] ALICE Collaboration, ϒ production and nuclear modification at forward rapidity in Pb-Pb collisions at ffiffiffiffiffiffiffi ffi s NN p ¼ 5.02 TeV, Phys. Lett. B 822, 136579 (2021). [29] ATLAS Collaboration, Production of ϒðnSÞ mesons in Pb þ Pb and pp collisions at ffiffiffiffiffiffiffi ffi s NN p ¼ 5.02 TeV, Phys. Rev. C 107, 054912 (2023). [30] L. Adamczyk et al. (STAR Collaboration), Suppression of ϒ production in d þ Au and Au þ Au collisions at ffiffiffiffiffiffiffi ffi s NN p ¼ 200 GeV, Phys. Lett. B 735, 127 (2014); 743, 537(E) (2015). [31] B. E. Aboona et al. (STAR Collaboration), Measurement of sequential ϒ suppression in Au þ Au collisions at ffiffiffiffiffiffiffi ffi s NN p ¼ 200 GeV with the STAR experiment, Phys. Rev. Lett. 130, 112301 (2023). [32] A. Adare et al. (PHENIX Collaboration), Measurement of ϒð1S þ 2S þ 3SÞ production in p þ p and Au þ Au collisions at ffiffiffiffiffiffiffiffi s NN p ¼ 200 GeV, Phys. Rev. C 91, 024913 (2015). [33] D. d'Enterria and C. Loizides, Progress in the glauber model at collider energies, Annu. Rev. Nucl. Part. Sci. 71, 315 (2021). [34] HEPData record for this analysis, 10.17182/hepdata.130959 (2023). [35] CMS Collaboration, The CMS experiment at the CERN LHC, J. Instrum. 3, S08004 (2008). [36] CMS Collaboration, Performance of the CMS Level-1 trigger in proton-proton collisions at ffiffi ffi s p ¼ 13 TeV, J. Instrum. 15, P10017 (2020). [37] CMS Collaboration, The CMS trigger system, J. Instrum. 12, P01020 (2017). [38] CMS Collaboration, Performance of electron reconstruction and selection with the CMS detector in proton-proton collisions at ffiffi ffi s p ¼ 8 TeV, J. Instrum. 10, P06005 (2015). [39] CMS Collaboration, Performance of the CMS muon detector and muon reconstruction with proton-proton collisions at ffiffi ffi s p ¼ 13 TeV, J. Instrum. 13, P06015 (2018). [40] CMS Collaboration, Performance of photon reconstruction and identification with the CMS detector in proton-proton collisions at ffiffi ffi s p ¼ 8 TeV, J. Instrum. 10, P08010 (2015). [41] CMS Collaboration, Description and performance of track and primary-vertex reconstruction with the CMS tracker, J. Instrum. 9, P10009 (2014). [42] CMS Collaboration, Observation and studies of jet quenching in Pb-Pb collisions at ffiffiffiffiffiffiffi ffi s NN p ¼ 2.76 TeV, Phys. Rev. C 84, 024906 (2011). [43] CMS Collaboration, Measurement of the azimuthal anisotropy of ϒð1SÞ and ϒð2SÞ mesons in PbPb collisions at ffiffiffiffiffiffiffiffi s NN p ¼ 5.02 TeV, Phys. Lett. B 819, 136385 (2021). [44] CMS Collaboration, Luminosity measurement in protonproton collisions at 5.02 TeV in 2017 at CMS, CMS Physics Analysis Summary, Report No. CMS-PAS-LUM-19-001, 2021, http://cds.cern.ch/record/2765655. [45] CMS Collaboration, CMS luminosity measurement using nucleus-nucleus collisions at ffiffiffiffiffiffiffiffi s NN p ¼ 5.02 TeV in 2018, CMS Physics Analysis Summary, Report No. CMS-PAS-LUM-18-001, 2022, https://cds.cern.ch/record/2809613. [46] CMS Collaboration, Precision luminosity measurement in proton-proton collisions at ffiffi ffi s p ¼ 13 TeV in 2015 and 2016 at CMS, Eur. Phys. J. C 81, 800 (2021). [47] CMS Collaboration, Charged-particle nuclear modification factors in PbPb and pPb collisions at ffiffiffiffiffiffiffiffi s NN p ¼ 5.02 TeV, J. High Energy Phys. 04 (2017) 039. [48] CMS Collaboration, Technical proposal for the Phase-II upgrade of the Compact Muon Solenoid, CMS Technical Proposal, Report No. CERN-LHCC-2015-010, 2015. [49] CMS Collaboration, Transverse momentum and pseudorapidity distributions of charged hadrons in pp collisions at ffiffi ffi s p ¼ 0.9 and 2.36 TeV, J. High Energy Phys. 02 (2010) 041. [50] CMS Collaboration, Extraction and validation of a new set of CMS PYTHIA8 tunes from underlying event measurements, Eur. Phys. J. C 80, 4 (2020). [51] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Skands, An introduction to PYTHIA8.2, Comput. Phys. Commun. 191, 159 (2015). [52] CMS Collaboration, Measurement of the ϒð1SÞ, ϒð2SÞ, and ϒð3SÞ polarizations in pp collisions at ffiffi ffi s p ¼ 7 TeV, Phys. Rev. Lett. 110, 081802 (2013). [53] CMS Collaboration, ϒðnSÞ polarizations versus particle multiplicity in pp collisions at ffiffi ffi s p ¼ 7 TeV, Phys. Lett. B 761, 31 (2016).