First measurement of large area jet transverse momentum spectra in heavy-ion collisions

Jet production in lead-lead (PbPb) and proton-proton (pp) collisions at a nucleon-nucleon center-of-mass energy of 5.02 TeV is studied with the CMS detector at the LHC, using PbPb and pp data samples corresponding to integrated luminosities of 404 μb−1 and 27.4 pb−1, respectively. Jets with different areas are reconstructed using the anti-kT algorithm by varying the distance parameter R. The measurements are performed using jets with transverse momenta (pT) greater than 200 GeV and in a pseudorapidity range of |η| < 2. To reveal the medium modification of the jet spectra in PbPb collisions, the properly normalized ratio of spectra from PbPb and pp data is used to extract jet nuclear modification factors as functions of the PbPb collision centrality, pT and, for the first time, as a function of R up to 1.0. For the most central collisions, a strong suppression is observed for high-pT jets reconstructed with all distance parameters, implying that a significant amount of jet energy is scattered to large angles. The dependence of jet suppression on R is expected to be sensitive to both the jet energy loss mechanism and the medium response, and so the data are compared to several modern event generators and analytic calculations. The models considered do not fully reproduce the data.


JHEP05(2021)284
p T imbalance [28][29][30], modifications of the jet yield in the medium [31][32][33][34][35][36][37], electroweak boson-jet p T imbalance [38,39], jet fragmentation functions [40][41][42], missing p T in dijet systems [28,29,43], jet-track correlations [44], and the radial p T profile of tracks within jets [45][46][47][48][49] have been studied. Complementary to these measurements, inclusive jet spectra reconstructed using different distance parameters R in the anti-k T algorithm are of great interest because they are less sensitive to hadronization effects than observables involving individual final-state hadrons. The value of R defines the area of the reconstructed jet. By varying R, different fractions of energy from the quenched jet and the medium response will be included in the reconstructed jet. A differential study of the suppression versus R provides new sensitivity to the QGP properties [50] and to the underlying jet quenching mechanism. In particular, theoretical models and generators based on perturbative QCD [50][51][52] and anti-de Sitter/conformal field theory correspondence [53] predict different dependences of the jet suppression on R.
Modifications to jet production can be quantified by the ratio of the inclusive jet yields per event in nucleus-nucleus (AA) collisions (N AA ) to those in proton-proton (pp) collisions (N pp ), scaled by the mean number of binary nucleon-nucleon (NN) collisions ( N coll ) [54]. This ratio is called the nuclear modification factor R AA and is defined as where p jet T is the transverse momentum of the jet. The R AA is typically measured in bins of centrality, which characterizes the degree of overlap of the two colliding lead nuclei [29,55]. The nuclear overlap function T AA is defined as the ratio of N coll to the total inelastic pp cross section, T AA = N coll /σ pp inel , and can be calculated from a Glauber model of the nuclear collision geometry [55]. If the ratio is less than one, it indicates a transfer of energy out of the jet cone. Measurements of the dependence of jet spectra and nuclear modification factors on the jet distance parameter R can help differentiate between competing models of parton energy loss mechanisms [56].
In studies of jet suppression from LHC Run 1 with lead-lead (PbPb) collisions at a nucleon-nucleon center-of-mass energy of √ s NN = 2.76 TeV, it was shown that the level of suppression is roughly independent of p jet T in the range p jet T = 200-400 GeV [33]. This suggests that the shape of the spectra is not significantly modified, and the modifications are predominantly through the overall number of jets. However, these initial measurements were statistically limited. At √ s NN = 5.02 TeV, this measurement can be extended to higher p T . Furthermore, at this higher center-of-mass energy, partons traverse a medium of higher density and temperature.
In this paper, measurements of jet R AA at p jet T > 200 GeV using PbPb collisions at √ s NN = 5.02 TeV are reported. The jets are reconstructed using the anti-k T algorithm [27] with R varying between 0.2 and 1.0. The results are presented as a function of p jet T in bins of PbPb event centrality.

The CMS apparatus
The central feature of the CMS detector is a superconducting solenoid of 6 m internal diameter, providing a magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead tungstate crystal electromagnetic calorimeter (ECAL), and a brass and scintillator hadron calorimeter (HCAL), each composed of a barrel and two endcap sections. Hadron forward (HF) calorimeters extend the pseudorapidity coverage up to |η| = 5.2 and are used for event selection. In the case of PbPb events, the HF signals are also used to determine the centrality class of the event. In the barrel section of the ECAL, an energy resolution of about 1% is achieved for unconverted or late-converting photons that have energies in the range of tens of GeV. The remaining barrel photons have a resolution of about 1.3% up to |η| = 1, rising to about 2.5% at |η| = 1.4. In the endcaps, the resolution of unconverted or late-converting photons is about 2.5%, while the remaining endcap photons have a resolution between 3 and 4% [57]. When combining information from the entire detector, the jet energy resolution amounts typically to 15% at 10 GeV, 8% at 100 GeV, and 4% at 1 TeV, to be compared to about 40, 12, and 5% obtained when the ECAL and HCAL calorimeters alone are used [58]. Muons are detected in gas-ionization chambers embedded in the steel flux-return yoke outside the solenoid. The silicon tracker measures charged particles within |η| < 2.5. It consists of 1440 silicon pixel and 15 148 silicon strip detector modules. For nonisolated particles of 1 < p T < 10 GeV and |η| < 1.4, the track resolutions are typically 1.5% in p T and 25-90  µm in the transverse (longitudinal) impact parameter [59]. Events of interest are selected using a twotiered trigger system [60]. The first level, composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events at a rate of around 100 kHz within a time interval of less than 4 µs. The second level, known as the high-level trigger, consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing, and reduces the event rate to around 1 kHz before data storage. A more detailed description of the CMS detector, together with a definition of the coordinate system used and the relevant kinematic variables, can be found in ref. [61].

Event selection
The event samples are recorded with dedicated jet triggers with different p jet T thresholds, the smallest of which is p jet T > 80 GeV [39]. The efficiencies of the triggering algorithms are evaluated in data and are found to reach unity in both pp and PbPb collisions for jets considered in this paper (p jet T > 200 GeV). A number of requirements are made to the events to remove non-collision events (e.g., beam-gas interactions) and to select only inelastic hadronic collisions [39,62]. Both pp and PbPb events are required to have at least one reconstructed primary interaction vertex with a distance from the center of the nominal interaction region of less than 15 cm along the beam direction. In addition, in PbPb collisions the shapes of the clusters in the pixel detector have to be compatible with those produced by a genuine collision [63]. The PbPb collision events are also required to -3 -  [54].

JHEP05(2021)284
have at least three towers in each of the HF detectors with energy deposits of more than 3 GeV per tower. These criteria select 99% of inelastic hadronic PbPb collisions [29].
The collision centrality for PbPb events is determined using the total sum of transverse energy from the calorimeter towers in the HF region. The transverse energy distribution is used to divide the event sample into bins of percentage of the total hadronic interaction cross section [29]. The results in this paper are presented in four centrality intervals, where 0% corresponds to a full overlap of the two nuclei: 0-10, 10-30, 30-50, and 50-90%. The corresponding T AA and N coll values used in this paper for the centrality intervals are listed in table 1.

Monte Carlo simulations
Several Monte Carlo (MC) simulated jet event samples are used to evaluate background components, efficiencies, misreconstructed jet rates (arising from upward fluctuations of the underlying event (UE) without a corresponding hard parton), jet energy corrections and jet energy resolutions (JER). Proton-proton collisions are generated with pythia 8.212 [64], with the UE tune CUETP8M1 [65], as well as with pythia6 [66], with the UE tune Z2 [67] with PDF set CTEQ6L1 [68]. For the PbPb MC samples, each pythia signal event is embedded into a PbPb collision event generated with hydjet v1.8 [69], which is tuned to reproduce global event properties such as the charged-hadron p T spectrum and particle multiplicity. The detailed simulation of the CMS detector response is performed using the Geant4 package [70].

Jet reconstruction and underlying event subtraction
Particle candidates are reconstructed with the particle-flow (PF) algorithm [58], where information from different parts of the detector are combined to form an optimized description of the event. Jets are clustered from the PF candidates using the anti-k T algorithm with distance parameters of R = 0.2, 0.3, 0.4, 0.6, 0.8, and 1.0, as implemented in the FastJet framework [27,71].

JHEP05(2021)284
One of the main challenges to reconstructing jets in heavy-ion collisions is the additional soft UE coming from the QGP. In order to subtract the soft UE in PbPb collisions on an event-by-event basis, an iterative algorithm [72] is employed. The mean value, E PF , and dispersion, σ(E PF ), of the transverse energies from the PF candidates are calculated in a number of η bins [29,38,73] for each event. Then, a two-step procedure is employed to account for the azimuthal modulation of background activity arising from the bulk properties of the QGP. In the first step, the so-called event plane angles (Φ EP,2 , Φ EP,3 ) for the second-and third-order harmonics of the azimuthal distribution are derived from the HF calorimeters (3 < |η| < 5) [74]. This method of estimating the UE gives underlying energy estimations that are consistent with a previous analysis of photon-and Z-tagged jets in which event plane mixing was used [75]. The event plane angles are not corrected for detector effects since the only goal of this procedure is to obtain a better description of the modulation of the background level. For the second step, a fit to the azimuthal angle (φ, in radians) distribution of charged-hadron PF candidates with 0.3 < p T < 3.0 GeV and |η| < 1 is performed. No explicit exclusion of regions close to the jet is performed, since their effect on the extracted parameters is negligible. The functional form of the fit is as follows: where N 0 is the magnitude of average UE activity. The parameters v 2 and v 3 quantify the strengths of the collective behaviors of the UE known as "elliptic" and "triangular" flow, respectively. The event plane angles Φ EP,2 and Φ EP,3 are fixed to the result from the first step. A fit is performed per event to extract the parameters N 0 , v 2 , and v 3 . The fit is discarded if the minimum required number of candidates (at least 10 entries in each bin) are not met, or if the reduced χ 2 of the fit is greater than 2. In this case, the background energy density is estimated as a flat distribution in φ, without flow modulations.
An example of this procedure is shown for data in figure 1. The left plot shows the fit in the extraction region, along with a breakdown of the components of the fit. The right plot takes parameters extracted from mid-rapidity (|η| < 1) and renormalizes the function to data at forward-rapidity (1 < |η| < 2). A good agreement of the shape for background modulations in the two η ranges is observed.
Finally, the UE subtraction in PbPb collisions is performed using a constituent subtraction method [76]. This is a particle-by-particle approach that corrects the energy of each jet constituent based on the local average UE density ρ(η, φ). This density is assumed to factorize in η and φ according to the form Here ρ(η) encodes the variation of the UE density as a function of η, and the flow parameters are determined in the previous fit. The average UE density ρ(η) is calculated as the average energy in given η bins, excluding regions overlapping with jets. In pp collisions, where the UE level is negligible, jets are reconstructed without UE subtraction.  The v 2 (blue curve) and v 3 (yellow curve) of the flow components are shown, together with the total modulation used in the analysis to account for the background (red curve). The flow coefficients are extracted from the left plot and overlaid in the right plot.

Jet energy scale and resolution
Jet energy corrections are derived from simulation separately for pp and PbPb data following methods outlined in ref. [77]. The energy scales are verified with an energy balance method applied to dijet and photon+jet events in pp data. For this study, jets with |η jet | < 2 and (corrected) p jet T > 160 GeV are selected. The p jet T binning of the analysis is chosen based on the jet energy resolution (JER) for each cone size and centrality. For pp events, the JER varies by less than 10% for different values of R. These variations reflect how the probability for energy to move into or out of the jet cone changes with cone sizes. Figure 2 shows the PbPb jet energy scale (JES, upper), defined as the reconstructed p jet T divided by the generated p jet T , and JER (lower), for R = 0.2 (left) and R = 1.0 (right) as functions of the generated p jet T . The JES is rather flat vs. p jet T while JER decreases with p jet T . As expected, the resolution is worse for more central events and for larger values of R, because of the larger UE contribution that must be subtracted. For R ≤ 0.4 the difference found in both pp and PbPb simulations between the JES of generated and reconstructed p jet T is below 2% at mid-rapidity (|η| < 1) and of order 4% for (1 < |η| < 2).
A small nonclosure of up to 2% was observed for all values of R in the peripheral 70-100% PbPb bin, where the nonclosure is defined as the deviation of the corrected JES from unity. The UE in this bin is most comparable to that in pp collisions, and it is used to evaluate the performance of the jet algorithm with heavy ion reconstruction and subtraction in the absence of UE. This is necessary, as the difference in tracking and the subtraction of an UE in PbPb, compared to pp, results in modest performance changes even without a significant UE contribution. The φ-modulation of ρ shown in eq. (5.2) improves the jet resolution without introducing any biases to the energy scale. However, as can be seen in figure 2, there is evidence of over-subtraction for the largest values of R at low p jet T . This is because of uncertainties in the estimation of ρ. Errors in the estimation of ρ tend to be handled much better for small R, as the subtraction scales with the multiplicative area. This over-subtraction causes the nonclosure to reach up to 4% for R = 1.0 jets, as seen in figure 2. For smaller values of R the nonclosure is below 2%.
Another source of over-subtraction is caused by the flow-modulated subtraction. The minimum candidate count requirement for a good UE shape estimation does not account for the fact that jets could bias the fit. The over-subtraction occurs when a jet biases the flow modulation fit. While the fitting for φ modulation is turned off for events with a small number of tracks, events close to this threshold could still be affected by these biases, resulting in a source of nonclosure.

JHEP05(2021)284
Finally, for the most central events, where ρ(η, φ) ranges from 200-300 GeV per unit area and the fluctuations are the largest, there is a global underestimation of the true UE, particularly in the forward region (|η| > 1). This occurs because towers within jets are nominally excluded in the estimation of ρ to avoid bias from the hard process. In the most central events, upward fluctuations of the UE may cause some towers to be included in the jet and excluded from the UE. If too many towers are excluded, ρ is underestimated. This underestimation of ρ results in the largest uncertainty in the final R AA and spectra for the most central bins. It is mitigated by setting an upper limit on the number of excluded towers, with a cutoff that is tuned to achieve the best performance.

Unfolding
Raw spectra are unfolded according to response matrices constructed using pythia+hydjet MC for PbPb and pure pythia for pp results, in matched bins of p jet T , η jet , and for PbPb only, event centrality. The matrices are constructed with an N coll distribution that matches the expectations from data. The unfolding is done with the d'Agostini iteration with early stopping [78] as implemented in the RooUnfold package [79]. Examples of response matrices are shown in figure 3 for pp and 0-10% PbPb collisions with R = 0.2 and 1.0. Underflow bins are shown to account for bin migration. As expected, the matrices are more tightly diagonal for pp events than for PbPb, and for R = 0.2 than for R = 1.0. The unfolded pp and PbPb spectra are then used to construct the R AA distribution.

Systematic uncertainties
The systematic uncertainties in the spectra are estimated by varying analysis parameters one at a time within a reasonable range, propagating the change through the full analysis chain, and then considering the deviation from the nominal results. For R AA , any correlation between the uncertainties in the pp and PbPb spectra is accounted for by simultaneously changing the same parameter in the pp and PbPb analyses, calculating a new R AA and taking the difference from the nominal result. This procedure produces a significant reduction in the uncertainty from data-simulation differences that impact JES and JER since the pp and PbPb were taken in run periods separated by just a few days. For ratios of R AA between different jet radii, the luminosity and the T AA uncertainties cancel.
Finally, in the R AA ratio between different radii, and the pp ratios of spectra between radii, there are statistical cancellations as the same jet may contribute to multiple R spectra. These are accounted for by comparing ratios of spectra in pseudo-experiments generated independently from the spectra and those generated with the correlation between different R taken from the data. Figure 4 shows the principal systematic uncertainties as a function of p jet T for R = 0.2 and 1.0, and for pp and PbPb collisions as a function of centrality. The dominant uncertainty arises from the JES. This tends to increase with p jet T and centrality but does not have a strong dependence upon R. The unfolding and JER uncertainties tend to decrease with p jet T and increase with centrality and R. The T AA uncertainty decreases from peripheral to central events and is independent of p jet T and R. The origins of these uncertainties are listed below in order of importance for R AA .
1. Jet energy scale. The uncertainty ranges from 15 to 20% and is dominated by the data-simulation difference. It consists of several components, summed in quadrature: (a) Nonclosure in simulation. This uncertainty is evaluated as a function of centrality and η but independently of p jet T . It is estimated by varying data by the observed nonclosure in simulation, see figure 2, and then propagating this change through the analysis chain. In pp and peripheral (50-90%) PbPb collisions, a 1% variation is made for all η. For 10-50% centrality, the variation is 1% within |η jet | < 1 and 2% for 1 < |η jet | < 2. For the most central (0-10%) events, a 2% variation is made for jets within |η jet | < 1 and a 4% variation for 1 < |η jet | < 2.
(b) Data-simulation differences. A flat 2% variation is performed in all bins following the procedure used for the nonclosure uncertainties above. This uncertainty is dominant in the pp spectra, and comparable to the nonclosure uncertainty in semicentral and semiperipheral PbPb bins. (c) Differences from the UE description between data and simulation. These differences are extracted by comparing random cone mean/widths between data and simulation, and the full difference is taken as a systematic uncertainty. As this is the centrality-dependent component of the JES, it does not cancel in the ratios between pp and PbPb data, and only cancels partially in R-dependent ratios of R AA .
(a) The JER uncertainty is extracted from simulation. This is subdominant compared to the data-simulation differences for spectra, but does not cancel in R AA .
(b) Jet energy resolution from data-simulation differences. The resolution in data is found to be 10 to 15% worse than that in simulation. To propagate this uncertainty, the simulation is first smeared by 10%, such that central values are closer to those in data. The systematic uncertainty is estimated by applying an additional smearing on top of these new central values such that the resolution is increased by 10% in all bins. The effect is subdominant in part because the p jet T binning was chosen to minimize bin migration. Furthermore, there is partial cancellation in R AA , coming from the constant and stochastic terms of the resolution, which are partially shared between the pp and PbPb data. 4. Integrated luminosity and T AA . The uncertainty in the integrated luminosity for pp collisions is 2.3% [81]. For the T AA , the relative uncertainties vary between 3% for the 0-10% bin, to 11% in the most peripheral 50-90% bin [54]. The absolute uncertainties for each of the four values are listed in table 1.

5.
Misreconstructed jets which arise from fluctuations of the UE. The contamination from these jets is evaluated from simulation, and it is found to be negligible in the considered kinematic range.

Results
The unfolded jet spectra as functions of p jet T for R = 0.2 and 1.0 for both pp and PbPb collisions of various centralities are shown in figure 5. The lower bound of p jet T is chosen based on the observed noise level for each centrality class, and the upper bound is driven by the amount of statistics. The upper plot of figure 6 shows the ratio of spectra of jets with different radii in pp collisions, normalized to the spectrum for R = 1.0. The number of jets with a given p jet T increases with the size of the jet cone. The increase of jet yield with R becomes weaker at higher values of p jet T suggesting that jets become narrower as p jet T increases. Figure 6 also shows predictions using the pythia6 and pythia8 MC generators. Both generators capture the trends of the data but pythia8 is closer to the scale of the data. The lower plot of figure 6 shows the ratios of the jet spectra from pythia to the data spectrum for R = 0.2 and 0.4. For pythia6 the ratio rises with p jet T for both values of R. The pythia8 ratios show little dependence on p jet T and are generally closer to unity than those of pythia6. The R AA factors compare PbPb data to the scaled pp reference. Figure 7 shows R AA , the ratio of PbPb data to a scaled pp reference, as functions of p jet T , jet radius, and For all values of R, R AA for the most peripheral collisions (50-90%) is independent of p jet T and consistent with unity after considering the T AA uncertainty. In the most central bin, a strong suppression of the PbPb data (≈0.6-0.7) is observed, which is well outside the systematic uncertainties. However, there are hints of an increasing R AA with p jet T for the smaller values of R in the central bins, with values up to 0.8 for jets with p jet T > 500 GeV. To highlight the jet radius dependence of the jet R AA , the ratios of R AA for a given R with respect to R = 0.2 are presented in figure 8. This observable is particularly sensitive to the recovery of the quenched energy and the presence of the medium response [56]. For 400 < p jet T < 500 GeV, the R AA ratios are above unity and increase with p jet T in both the 0-10% and 10-30% centrality intervals. On the other hand, for p jet T > 500 GeV, the AA is close to unity or slightly below it for the 0-10% and 10-30% centrality intervals, respectively.  Jewel increases greatly as R increases. For R = 1.0 the predictions without recoil are a factor of four below the default mode with recoil. The Jewel predictions with recoil are significantly below the data for R = 0.2 but come increasingly close to the data as R increases.
Predictions from pyquen are shown with (the default, shown in teal) and without (turquoise) medium-induced wide-angle radiation. The default pyquen generator overpredicts R AA particularly for smaller values of R and p jet T . The inclusion of wide-angle radiation lowers the predictions for R AA particularly for smaller R sizes and brings the pyquen predictions closer to the data, showing the importance of the medium effects.  Figure 11 shows a comparison of several models to R AA as functions of p jet T and R. The Hybrid model [56] combines a perturbative description of the weakly coupled physics of jet production and evolution, with a gauge/gravity duality description of the strongly coupled dynamics of the medium, and the soft-gluon exchanges between the jet and medium. As the jet passes through and deposits energy into the hydrodynamic medium, a wake is left behind the jet. The Hybrid model (dark orange) tends to under-predict R AA at high p jet T . Calculations without a wake (brown) and with only the positive contribution of the wake (yellow) are also shown. These two are not physical and are included here only for better understanding of the effect of the wake contribution. The effect of the wake is more important at large R and lower p jet T . In the Linear Boltzmann Transport (lbt) model [83], the effects of recoil thermal partons and their propagation in the dense medium are described by a 3+1D viscous relativistic hydrodynamic model. Predictions from lbt are shown in figure 11 with and without the medium response. It is clear that the medium response becomes more and more dominant as the size of the jet cone increases. A similar effect is seen for the jetcoupled fluid model [52,84,85] Ccnu. Although predictions are only available for a limited p jet T range, it is clear from comparing the blue and violet points in figure 11 that the hydrodynamic component of Ccnu becomes increasingly important with increasing R.
The predictions from Martini [86] (Modular Algorithm for Relativistic Treatment of Heavy IoN Interactions) are shown as purple boxes in figure 11. The model follows a hybrid approach where it embeds the high energy parton into an evolving hydrodynamic medium, and the shower evolution of the jet is modified following the McGill-AMY formalism [87][88][89][90][91]. The Martini generator predicts a larger increase of jet R AA ratio as a function of R than what is observed in data.  The same data in figures 11 and 12 are also compared to additional models. Figure 13 shows R AA vs. p jet T for several values of R and for the top 0-10% centrality class as well as several predictions from generators and analytic calculations. The gray boxes in figure 13 are predictions from a jet factorization model based on a phenomenological approach to establish QCD factorization of jet cross sections in heavy ion collisions [92]. Mediummodified jet functions are extracted from jet nuclear modification factors at smaller jet distance parameter values (R = 0.2 and 0.4) and predictions are made for larger distance parameter values. At R < 0.4, the data are described reasonably well by the factorization model. However, the model tends to underpredict R AA at larger R values. The data are also compared to the coherent antenna bdmps calculations [93] (orange), which is an analytical approach that resums multiple emissions to leading-logarithmic accuracy including both radiative energy loss and color coherence effects [94][95][96]. The predictions are in general agreement with the R AA data. Finally, calculations based on a soft collinear effective theory with Glauber gluon interactions Scet [50], are also compared to the data. The Scet calculations with collisional energy loss [97,98] (navy blue) are slightly below the R AA measurements while those without collisional energy loss (sky blue) are consistent with the data. Figure 14 shows R R AA /R R=0.2 AA vs. R for several values of p jet T together with predictions from the Scet, bdmps and jet factorization models. The bdmps (orange) and Scet predictions (sky blue and navy blue) are consistent with the data but the factorization calculations (gray) decrease with R in contrast to the data. central PbPb collisions, a strong suppression is observed for jets with high transverse momentum reconstructed with all distance parameters. Predictions from quenched jet event generators, theoretical models, and analytical calculations are compared to these results. The new data place further constraints on the underlying jet quenching mechanisms. While state of the art models have made important progress, significant tension remains in view of the large area jet data presented here.

JHEP05(2021)284
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

JHEP05(2021)284
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.