Midrapidity dip and limiting fragmentation of charged hadrons from RHIC to LHC energies

Introduction

The charged-hadron pseudorapidity density $dN_{\mathrm{ch}}/d\eta$ is one of the most precisely measured global observables in heavy-ion collisions. It is available with large acceptance from RHIC to the highest LHC energies [1, 2, 3, 4, 5, 6]. In central collisions it has a shallow depression at $\eta= 0$ that becomes more pronounced with energy. Pseudorapidity ignores the particle mass: at fixed transverse momentum the Jacobian $dy/d\eta= p/E$ is smallest at midrapidity, so a rapidity density $dN_{\mathrm{ch}}/dy$ that is flat near $y = 0$ produces a dip in $dN_{\mathrm{ch}}/d\eta$ by kinematics alone. How much of the observed dip this accounts for depends on the transverse-momentum spectra, the particle composition [7] and the flatness of $dN_{\mathrm{ch}}/dy$ near $y = 0$.

Earlier studies came to different conclusions, and they used different Jacobians. In the three-source relativistic diffusion model (RDM), an early LHC prediction with a single effective mass and transverse momentum found the size of the dip “essentially determined by the relative particle content in the three sources, not by the Jacobian” [8]. Later RDM analyses used an effective Jacobian determined from identified-particle data and found it “essential to generate the dip” [9]; near midrapidity “it generates the dip” [10]. Hoelck and Wolschin, solving a two-dimensional diffusion model for pions, kaons and protons with the full momentum- and mass-dependent Jacobian, found the dip “mainly due to the Jacobian transformation” [11]. ALICE converted its 2.76 TeV $dN_{\mathrm{ch}}/d\eta$ to $dN_{\mathrm{ch}}/dy$ with a Jacobian from its measured midrapidity $\pi$, $K$ and $p$ spectra and found the result consistent with a flat rapidity plateau of about $\pm1.5$ units around $y = 0$, and also with a wide Gaussian of width 3.86 [4]. Event generators provide the full momentum of every particle, but their spectra need not match the data: HIJING [12, 13], for example, has softer spectra than measured and therefore a deeper kinematic dip.

The question is addressed here without a model of particle production: do the measured $dN_{\mathrm{ch}}/d\eta$ require any midrapidity structure in $dN_{\mathrm{ch}}/dy$ once the exact $y \to\eta$ kinematics with measured identified spectra are applied? Candidate $dN_{\mathrm{ch}}/dy$ shapes are folded through a kernel built from the measured $\pi$, $K$ and $p$ spectra and their measured rapidity dependence. The null hypothesis is a family of smooth dip-free shapes; the alternative adds a dip of free depth and width. The test statistic is calibrated with pseudo-experiments, its sensitivity is quantified by injecting dips, and the decision rule was fixed before any $p$-value was computed, although it was not registered publicly; later changes are listed in Appendix § 10. The analysis is not a replacement for momentum-space models such as the RDM. It tests what the data require, and a $dN_{\mathrm{ch}}/dy$ without a dip in rapidity, such as a dominant midrapidity source with fragmentation sources, is one of the shapes it allows.

The same machinery bears on a second question. Limiting fragmentation [14], the energy independence of particle production near beam rapidity, has been observed at RHIC [1], and ALICE found its 2.76 TeV data consistent with it within the uncertainties of extrapolating them toward beam rapidity [4]. LHC data stop several units before beam rapidity, so statements at the LHC rest on an extrapolation. Double-Gaussian extrapolations of the ALICE data have been read as violating limiting fragmentation at 5.02 TeV when the Gaussians are fitted directly in $\eta$ [15], and as fulfilling it when they are fitted in rapidity and transformed with a Jacobian [9]. In the RDM the fragmentation tail is shaped by the Dirichlet boundary condition of the numerical solver, and at LHC energies that analysis of limiting fragmentation “remains a model-dependent prediction” [10]. The data alone are found to leave this extrapolation largely unconstrained, and the constraint provided by an energy sum rule calibrated on RHIC data is examined.

Section 2 lists the data and Sect. 3 describes the kernel and the test. Section 4 gives the results on the dip, Sect. 5 closure tests with event generators and Sect. 6 a comparison of identified and inclusive BRAHMS yields. Limiting fragmentation is discussed in Sect. 7, limitations in Sect. 8, and Sect. 9 summarizes.

Data

Pseudorapidity densities. For Au+Au collisions the data are the PHOBOS 0–6% central $dN_{\mathrm{ch}}/d\eta$ at $\sqrt{s_{NN}} =$ 19.6, 62.4, 130 and 200 GeV from the collaboration’s published data tables [1, 2, 3]. The tables give systematic uncertainties only, which PHOBOS quotes at 90% confidence level, with negligible statistical errors [1, 2]. We average the upper and lower uncertainties and divide by 1.645, the half-width of a central 90% interval of a normal distribution in units of its standard deviation, to obtain one standard deviation; this assumes Gaussian uncertainties (the effect of taking them as quoted is given in Sect. 7.4). The Pb+Pb data are ALICE 0–5% at 2.76 TeV [4] and 5.02 TeV [5], and CMS 0–5% at 5.36 TeV [6]. The LHC tables are published symmetrized in $\eta$, so each $|\eta|$ is used once. The data reach $|\eta| = 5.4$ (ALICE, 2.76 TeV), 4.9 (ALICE, 5.02 TeV) and 2.5 (CMS, 5.36 TeV). The mean numbers of participants used in Sect. 7, the collaborations’ Glauber-model estimates [16], are 337, 340 and 344 at 19.6, 130 and 200 GeV [1], and 382.8, 389.3 and 382.3 at 2.76, 5.02 and 5.36 TeV [4, 5, 6].

Identified spectra. The kernel uses midrapidity $\pi^\pm$, $K^\pm$, $p$ and $\bar{p}$ spectra for 0–5% central collisions (the PHOBOS $dN_{\mathrm{ch}}/d\eta$ are 0–6%): ALICE at 2.76 and 5.02 TeV [17, 18] (the 5.02 TeV spectra also serve at 5.36 TeV), PHENIX at 130 [19] and 200 GeV [20], and STAR at 19.6 GeV [21]. At 62.4 GeV the STAR spectra [22] are continued to higher $p_T$ with STAR [23] ($\pi$, $p$) and BRAHMS [24] ($K$) data (0–10%), scaled to the 0–5% spectra in the overlap so that only their shapes enter. The STAR spectra alone at 62.4 GeV [22] and 200 GeV [25] serve as a variant. The rapidity dependence of $\langle p_T\rangle$ and of the particle ratios is taken from BRAHMS at 200 GeV [26, 27].

Method

The $y \to\eta$ kernel

For a charged-hadron rapidity density $g(y) = dN_{\mathrm{ch}}/dy$,

with $y = \sinh^{-1}(p_T \sinh\eta/m_T)$, $m_T = \sqrt{p_T^2 + m_i^2}$, $p = p_T \cosh\eta$ and $E = \sqrt{p^2 + m_i^2}$ for a hadron of species $i \in\{\pi, K, p, \text{other}\}$ and mass $m_i$. Here $\phi_i$ is the normalized $p_T$ spectrum of species $i$ and $x_i$ its share of the charged yield, normalized to $\sum_i x_i = 1$ at every $y$. The measured spectra are used directly inside their $p_T$ range, and a Lévy–Tsallis fit [28, 29], $dN/dp_T \propto p_T[1 + (m_T - m_i)/(nT)]^{-n}$, outside it. Charged hyperons ($\Sigma^\pm$, $\Xi^-$, $\Omega^-$ and their antiparticles) are given 2% of the charged yield, an estimate from the measured $\Xi$ and $\Omega$ [30] and $\Lambda$ [31] yields at 2.76 TeV, with $\Sigma^{\pm}$ estimated from $\Lambda$ through isospin symmetry and an assumed $\Sigma^{0}/\Lambda$ ratio; they have the proton spectral shape and the $\Sigma$ mass, bracketed by variants with 0 and 4%. The rapidity dependence is modeled as $\phi_i(p_T;y)=\phi_i^0(p_T/\lambda_i(y))/\lambda_i(y)$, where $\phi_i^0$ is the midrapidity spectrum. The scale $\lambda_i(y)$ and the ratios $x_i(y)/x_{\pi}(y)$ are parametrized as $1+c_i y^2$, fitted to the BRAHMS points at 200 GeV in units of $y/y_{\mathrm{beam}}$, where $y_{\mathrm{beam}}=\cosh^{-1}(\sqrt{s_{NN}}/2m_p)$ is the beam rapidity, and held at their values at the largest measured $|y|$ beyond. For protons $\lambda_p(y)$ is not taken from the BRAHMS proton spectra, which cover four rapidity slices [27]; the kaon curve is used.

The kernel variants are: no rapidity dependence; the BRAHMS dependence in absolute $y$; an $m_T$-exponential low-$p_T$ extrapolation; a $\pm3.7\%$ global $p_T$ scale (the uncertainty of ALICE’s pion $\langle p_T\rangle$ [17]); the hyperon share; at RHIC a $\pm25\%$ proton share; and the STAR-only spectra at 62.4 and 200 GeV. At 2.76 TeV the kernel’s integrated yields ($\pi$ 1475, $K$ 218, $p$ 67) reproduce ALICE’s published values (1465, 218, 67) [17]. ALICE obtained its own $dN_{\mathrm{ch}}/dy$ from its $dN_{\mathrm{ch}}/d\eta$ with a Jacobian built from midrapidity identified spectra (at 5.02 TeV those measured at 2.76 TeV) [4, 5]. The ratio of the implied $dN_{\mathrm{ch}}/dy$ to it is 0.99–1.05 at 2.76 TeV and 0.95–1.04 at 5.02 TeV, a consistency check between two conversions rather than a comparison with a measured $dN_{\mathrm{ch}}/dy$. fig:1 compares the effective Jacobian of the kernel with those of the two event generators used in Sect. 5. At midrapidity EPOS4, which has radial flow, agrees with the measured-spectra kernel within −0.5 to +0.3% at 200 GeV and 2.76 TeV, whereas HIJING, with softer spectra, lies 4.4–4.9% lower and so produces a deeper kinematic dip.

Hypotheses

The null hypothesis H0 is a $dN_{\mathrm{ch}}/dy$ that is symmetric and non-increasing in $|y|$. Two families define it: a double Gaussian (dG), $A[G(y-t\sigma)+G(y+t\sigma)]$ with $G$ a Gaussian of width $\sigma$ and $0\leq t\leq1$, which keeps it unimodal, and a smeared box (box), $A[L_a(y-Y)-L_a(y+Y)]$ with $L_a(x)=1/(1+e^{x/a})$. A generalized Gaussian (gG), $A\exp[-(|y|/\sigma_g)^\beta]$, serves as a robustness check. The alternative H1 multiplies the H0 shape by $1-d\exp[-y^2/(2w^2)]$, a Gaussian dip of depth $0\leq d\leq1$ and width $w>0$. Apart from small lower limits on widths and on the exponent $\beta$ that keep the fits numerically defined, and $y_e\geq0$ below, the only parameter limits are those that define the hypotheses and positivity. Fits are multistart minimizations with MINUIT [32] through iminuit [33]; the lowest $\chi^2$ is kept, and H0 is refitted from the H1 solution, so that $\chi^2_{\mathrm{H0}}$ cannot be stuck above a minimum that H1 has found.

Each family also exists with a separate kinematic edge (dG+edge, box+edge, gG+edge), the shape multiplied by $L_a(|y|-y_e)$, so that the fall-off near beam rapidity is not tied to the central shape. These edge families, as opposed to the smooth families above, are used for the limiting-fragmentation test (Sect. 7), the closure tests (Sect. 5) and a check made after the results (Sect. 4); they are not part of the decision rule.

At RHIC the fit range of the dip test is $|\eta|\leq y_{\mathrm{beam}}$, with $y_{\mathrm{beam}}=3.04$, 4.20, 4.93 and 5.36 at 19.6, 62.4, 130 and 200 GeV. Near and beyond beam rapidity the kernel’s rapidity dependence, held constant beyond the BRAHMS range, is least reliable, so PHOBOS points beyond $y_{\mathrm{beam}}$ are not used in the dip test. The full range and $|\eta|\leq y_{\mathrm{beam}}-1$ are variants. The LHC data end well before beam rapidity and are used in full.

Uncertainties

The PHOBOS systematic uncertainties are not split into normalization, uncorrelated and shape parts. At RHIC the covariance is therefore

with $\sigma_j^{\mathrm{stat}}$ and $\sigma_j^{\mathrm{sys}}$ the statistical and systematic uncertainties of point $j$. Here $f$ is the fully correlated ends_mid=1 (normalization) share of the systematic variance; of the remainder, a share $u$ is uncorrelated from point to point and $1-u$ is correlated with length $\ell$ in $\eta$. 15 combinations are scanned: $f \in\{0, 0.5, 0.8\}$ with $u=1$, or with $u \in\{0.5, 0.2\}$ and $\ell\in\{0.5, 1.5\}$; $f=u=\ell=0.5$ is nominal. At the LHC the nominal covariance is built from the systematic components documented in each publication (Appendix A); its documented variants and the scan serve as variants. At nominal settings $\chi^2_{\mathrm{H0}}/\mathrm{ndf}$ of the decisive family (Sect. 3.4) is 0.05–0.31 at RHIC and 0.04–0.35 at the LHC: the published uncertainties are large compared with the point-to-point scatter of the data. At RHIC the more strongly correlated settings of the scan reach $\chi^2_{\mathrm{H0}}/\mathrm{ndf}=1.2$, and the dip-test results below are checked over the whole scan.

Test statistic, calibration and decision rule

The test statistic is $q=\chi^2_{\min}(\mathrm{H0})-\chi^2_{\min}(\mathrm{H1})\ge0$, with $\chi^2_{\mathrm{H0}}\equiv\chi^2_{\min}(\mathrm{H0})$. Its distribution under H0 is obtained from 200 pseudo-experiments per configuration, drawn from the best-fit H0 prediction with covariance $C$ and refitted exactly as the data; $p$ is the fraction of pseudo-experiments with $q\ge q_{\mathrm{obs}}$, and $q_{95}$, the 95th percentile, is 3.0–4.4 at nominal settings. The decision at each energy is drawn from the decisive family, the more flexible of the two H0 families (the one with the lower $\chi^2_{\mathrm{H0}}$ at nominal settings), because testing against a rigid null family would mistake its rigidity for evidence of a dip. A dip is required if $p<0.0027$ at nominal settings and $p<0.05$ in every kernel and covariance variant; it is not required if $p>0.05$ at nominal settings; anything else is inconclusive. With 200 pseudo-experiments the smallest resolvable $p$ is 0.005, above the threshold of 0.0027, so any configuration of the decisive family with $p<0.05$ was to be rerun with more pseudo-experiments before a decision. The $p$-values are those of the selected family; the selection is not repeated in the pseudo-experiments. At nominal settings both families give $p\ge0.065$ at every energy, so the decision does not depend on the selection.

Sensitivity

Dips of depth $d$ and width $w$ are injected into the best-fit H0 shape of the decisive family, on a grid of 7 widths ($w=0.3, 0.5, 0.75, 1, 1.5, 2, 3$) and 10 depths from 0.01 to 0.25. Pseudo-experiments are generated from each and analyzed with the same procedure, and $d_{\min}(w)$ is the depth at which half of them give $q>q_{95}$, i.e. $p<0.05$, interpolated linearly in $d$. $d_{\min}$ therefore refers to a truth that the null family describes exactly, folded with the nominal kernel; it does not include kernel or shape misspecification. With 200 pseudo-experiments per point, a rejection rate of 50% has a binomial uncertainty of 0.035. At the smallest injected depth, 1%, the rejection rate is 0.04–0.18 (0.05 without a dip, binomial uncertainty 0.015); it is largest for the narrowest widths at the LHC, where $d_{\min}$ is smallest, so a 1% dip is already partly detected there, but statistical scatter is not separated from this.

Results on the midrapidity dip

Figure 2 shows the data with the best dip-free fits and Figure 3 the implied rapidity densities. All three families describe the data at every energy. Table 1 summarizes the test, which covers 546 configurations with 200 pseudo-experiments each. At every energy the decision is “dip not required”. The smallest $p$ of the decisive family over all kernel, range and covariance variants is 0.23, and no configuration of a decisive family has $p<0.05$. The distributions of $q$ are shown in Figure 4.

$\sqrt{s_{NN}}$	family	$\chi^2_{H_0}$/ndf	$q_{\mathrm{obs}}$	$q_{95}$	$p$	$p_{\min}$	$d_{\min}$ (%), $w = 0.5/1/2
19.6 GeV	dG	1.4/27	0.21	3.9	0.75	0.23	15.8 / > 25 / > 25
62.4 GeV	box	11.3/39	0.00	3.9	1.00	0.81	6.7 / 12.5 / > 25
130 GeV	box	14.5/47	0.00	4.1	1.00	0.74	7.9 / 11.6 / > 25
200 GeV	box	7.6/51	0.00	4.4	1.00	0.75	7.4 / 9.6 / 24.4
2.76 TeV	dG	6.7/19	0.18	4.3	0.84	0.43	3.6 / 7.3 / 17.9
5.02 TeV	box	1.9/17	0.13	3.0	0.66	0.23	4.0 / 8.0 / > 25
5.36 TeV	box	0.4/10	0.32	3.6	0.54	0.24	6.1 / 8.2 / 21.6

“Not required” is not “excluded”. Figure 5 shows the sensitivity. For $w\le1$, $d_{\min}$ is 2.9–8.2% at the LHC and 5.9–12.5% at 62.4–200 GeV: deeper dips would have given $p<0.05$ in at least half of the pseudo-experiments. At 19.6 GeV, where the fit range is shorter and the uncertainties are larger, $d_{\min}$ is reached only for $w\le0.75$, with values of 11–21%. Broad depressions are only weakly constrained, because the dip-free families can absorb them: at $w=2$ a dip within the injected range is detectable only at 200 GeV, 2.76 TeV and 5.36 TeV, with $d_{\min}$ of 18–24%.

The only configurations with $p<0.05$ belong to the double-Gaussian family at 62.4–200 GeV, where it is not the decisive family, under covariance settings that leave little freedom for correlated shape errors: a normalization share $f=0.8$, or $f=0.5$ with an uncorrelated remainder ($u=1$); there $q$ reaches 13.4, with $p<0.005$ in the most extreme settings. Its fits sit at the unimodal boundary $t=1$, and its best “dip” is broader than the distribution itself. We interpret it as a reshaping of the whole distribution that the family cannot produce, rather than as a midrapidity structure. This is the situation the choice of the more flexible family was made for. As a check made after the results, the edge families (Sect. 3) were fitted to the RHIC data. They describe the data better than the smooth families do: at nominal settings $\chi^2_{\mathrm{H0}}$ of the smeared box falls from 7.6–14.5 to 3.0–4.2 at 62.4–200 GeV, suggesting that the data fall off faster near beam rapidity than the smooth shapes allow. They find no dip either: $q\le0.35$ at the four RHIC energies at nominal settings, and $q\le1.30$ at the covariance setting where the double Gaussian gives

its largest $q$. The test is one-sided: at RHIC the data have slightly less $\eta$-dip than the best dip-free prediction, within the uncertainties, and an excess at midrapidity is not tested.

Three further checks of the dip test were made after the results and are not part of the decision rule (Appendix § 10). First, PHOBOS corrections are presumably the same at $+\eta$ and $-\eta$, which the shape term of Sect. 3 does not assume; with the shape correlation taken as a function of $\lvert\lvert\eta_j\rvert-\lvert\eta_k\rvert\lvert$ instead, at nominal $f$, $u$ and $\ell$, the RHIC data give $p \ge0.77$, and $d_{\min}$ for $w \le1$ at 62.4–200 GeV becomes 6.6–15.3% (with 60 pseudo-experiments per depth), against 5.9–12.5% at nominal settings. Second, $d_{\min}$ refers to a truth that the null family describes. As a truth outside the families, the best fit of the smeared box with a kinematic edge at 200 GeV and 2.76 TeV, times a dip, was folded with the nominal kernel and fitted, without fluctuations, with the decisive family. Without a dip it gives $q \le0.57$, so the shape difference alone does not produce a false dip; the depth at which $q$ exceeds $q_{95}$, at $w = 0.5$ and 1, is 7.0–20.0% on the injection grid, 1.8–2.2 times $d_{\min}$. A dip in a truth that the families do not describe must therefore be about twice as deep to be detected. Third,

95% confidence-level upper limits on $d$ were computed with the CL$_{\mathrm{s}}$ method [34], with $q$ as the test statistic and the stored injection pseudo-experiments at nominal settings: $\mathrm{CL}_{\mathrm{s}}(d) = P(q \le q_{\mathrm{obs}} \mid d)/P(q \le q_{\mathrm{obs}} \mid0)$, with the limit where it falls to 0.05. They are 4.7–9.1% at $w = 0.5$ and 8.8–11.5% at $w = 1$ from 62.4 GeV to 5.36 TeV; at 19.6 GeV the limit is 18.3% at $w = 0.5$ and is not reached within the injected depths at $w = 1$. Like $d_{\min}$, they hold for a truth described by the null family and the nominal kernel.

Closure tests with event generators

The kernel relies on two approximations, the factorized rapidity dependence of the spectra and the extrapolation of the measured $p_{\mathrm{T}}$ ranges. Both, and the dip test itself, are tested in two event generators in which every final-state momentum is known: HIJING [12, 13], without collective flow, and EPOS4 [35, 36, 37, 38], whose core–corona procedure gives a viscous-hydrodynamic core with radial flow, followed by the UrQMD hadronic cascade [39, 40]. Both were run at 19.6, 62.4 and 200 GeV and at 2.76 and 5.02 TeV, for impact parameters $b \le3.3$ fm (Au+Au) and 3.5 fm (Pb+Pb), with $\pi^0$ and weak decays switched off and spectator nucleons removed. HIJING (through AliGenHijing of AliPhysics vAN-20260817-02-1; nuclear shadowing on, jet quenching off; spectators removed with AliGenHijing’s rule) has 2000 events per energy at 200 GeV and above and 1000 below; its final-state energy equals the participant energy $\langle N_{\mathrm{part}}\rangle\sqrt{s_{NN}}/2$ with a relative deviation below $7 \times10^{-4}$. EPOS4 (version 4.0.3, with the lead–lead example options of the EPOS4 authors: core–corona procedure, hydrodynamics with the x3ff equation of state, and the hadronic cascade) has at least 499 events per energy at 200 GeV and above and 300 below; the never-interacting spectator nucleons are identified in its event record, and the remaining final-state energy equals the participant energy within 0.74% per event (1 event with an energy violation was rejected). In these EPOS4 samples baryon number is not conserved exactly event by event; this does not enter the tests, because the energy rule takes its net baryons from data.

Folding each generator’s $dN_{\mathrm{ch}}/dy$ with its own $\eta$-versus-$y$ map reproduces its $dN_{\mathrm{ch}}/d\eta$ within 0.9% (HIJING) and 2.8% (EPOS4). The kernel is then rebuilt from the generator as from data: its midrapidity spectra, restricted to the $p_{\mathrm{T}}$ ranges measured at each energy and extrapolated in the same way, and its rapidity dependence fitted within the BRAHMS range and held constant beyond. At RHIC this is the generator’s own rapidity dependence at each energy; at the LHC the generator’s 200 GeV dependence is carried in $y/y_{\mathrm{beam}}$, as the data kernel carries BRAHMS’s. Folded with the generator’s true $dN_{\mathrm{ch}}/dy$, this kernel reproduces its $dN_{\mathrm{ch}}/d\eta$ at the $\eta$ points of the data within 1.8% (HIJING) and 1.1% (EPOS4) at the LHC (Table 2). At RHIC the deviations are largest at the most forward points, up to 10% in HIJING and 36% in EPOS4; within $|\eta| \le y_{\mathrm{beam}}-1$ they are 6–9% and 3–4%. In units of the data uncertainties the mean $\chi^2$ per point is at most 0.50 and 0.59. With full spectra and each generator’s exact rapidity dependence, which leaves only the factorized form $\phi_i(p_T/\lambda_i)/\lambda_i$, the deviations within $|\eta| \le y_{\mathrm{beam}}-1$ are 5–9% in HIJING and 3% in EPOS4, and up to 15% and 22% at all points. In HIJING the factorized $p_T$ scaling itself is therefore violated at this level; in EPOS4 part of the forward deviations comes from holding the rapidity dependence constant beyond the BRAHMS range.

generator	$\sqrt{s_{NN}}$	family	$\chi^2$/points	max. dev. (%)	$q_A$	$q_{95}$	measured range	$|y| \le y_{\mathrm{beam}}$
HIJING	19.6 GeV	dG	14.9/30	9.8	0.01	3.90	1.09	1.12
HIJING	62.4 GeV	dG	9.7/42	10.4	0.10	4.09	1.03	1.04
HIJING	200 GeV	box	11.5/54	10.0	12.71	4.39	0.99	1.00
HIJING	2.76 TeV	box	0.5/22	1.8	0.00	4.46	–	1.04
HIJING	5.02 TeV	box	0.1/20	0.7	0.00	2.98	–	1.05
EPOS4	19.6 GeV	box	17.7/30	19.1	0.22	3.94	1.22	1.24
EPOS4	62.4 GeV	box	15.6/42	36.2	0.00	3.95	1.15	1.17
EPOS4	200 GeV	dG	9.1/54	29.1	5.81	4.11	1.11	1.12
EPOS4	2.76 TeV	dG	0.4/22	1.1	0.01	4.27	–	1.28
EPOS4	5.02 TeV	box	0.3/20	0.9	0.24	2.98	–	1.31

Restricting the spectra to the narrow STAR $p_T$ ranges raises the extrapolated proton $\langle p_T\rangle$ at 200 GeV by 2% in HIJING and by 25% in EPOS4 (kaons 30%), whereas the nominal ranges change it by 1% and −2%. The generator with radial flow thus shows a bias of the same sign as the data, in which the STAR-only spectra give a proton $\langle p_T\rangle$ higher by 40% than the nominal spectra at 200 GeV (Sect. 8); HIJING, with softer spectra, shows it only weakly.

The dip test applied to noise-free generator pseudo-data with the uncertainties of the real data, whose $q$ (denoted $q_A$) approximates the median $q$, finds no dip at four of the five energies in both generators. At 200 GeV both give a false dip, $q_A$ above $q_{95}$: $q_A=12.7$ against 4.4 in HIJING (decisive family box) and $q_A=5.8$ against 4.1 in EPOS4 (dG; Table 2). In EPOS4 it is driven by the most forward points: restricting the fit to $|\eta| \le y_{\mathrm{beam}}-1$ removes it ($q_A=0.4$), as do the edge families (0.1). In HIJING, whose $dN_{\mathrm{ch}}/dy$ has a distinct forward component, neither does ($q_A=12.6$ and 12.7); only restricting the fit to $|\eta| \le3$ removes it ($q_A=1.9$; EPOS4 0.4). With each generator’s own $\eta$-versus-$y$ map instead of the rebuilt kernel, $q_A=10.5$ in HIJING, still above $q_{95}$, and 4.1 in EPOS4, just below it: the false dip in HIJING comes from the dip-free families’ failure to describe its forward $dN_{\mathrm{ch}}/dy$, and in EPOS4, where $q_A$ exceeds $q_{95}$ only marginally, partly from the kernel. HIJING’s true $dN_{\mathrm{ch}}/dy$ itself has a midrapidity dip of 1.4%, far below the sensitivity at this energy. The procedure can therefore produce a false dip at top RHIC energy, and a positive result would need checks with reduced fit ranges. In the data no dip is found with either the smooth or the edge families (Sect. 4).

Identified and inclusive BRAHMS yields at 200 GeV

BRAHMS measured both identified $\pi$, $K$ and $p$ rapidity densities and the charged $dN_{\mathrm{ch}}/d\eta$ at 200 GeV [26, 27, 41]. Folding the identified rapidity densities through the kernel (Figure 6) gives 0.97 of BRAHMS’s own $dN_{\mathrm{ch}}/d\eta$ at $\eta=0.15$, falling to 0.66 at forward pseudorapidity, for every kernel variant. The ratio of the BRAHMS and PHOBOS charged distributions is 0.88–1.10. The identified yields are thus narrower in rapidity than the inclusive charged distribution of the same experiment, and BRAHMS identified $dN_{\mathrm{ch}}/dy$ does not provide a reliable check of an $\eta\to y$ inversion beyond $|y| \approx1$. The origin is not resolved here; it may lie in the $p_T$ extrapolation of the forward identified spectra, in charged species other than $\pi$, $K$ and $p$, or in the forward systematics of the inclusive measurement. The kernel uses only the relative rapidity dependence of BRAHMS’s $\langle p_T\rangle$ and particle ratios, not the identified yields; its rapidity dependence is tested in two generators (Sect. 5) and varied between none, scaled and absolute.

Limiting fragmentation at the LHC

Limiting fragmentation is tested in the variable $\eta' = \eta-y_{\mathrm{beam}}$, in which it predicts $(2/\langle N_{\mathrm{part}}\rangle)\,dN_{\mathrm{ch}}/d\eta$ to be independent of energy. The most forward LHC points lie at $\eta' = -2.61$ (2.76 TeV), $-3.71$ (5.02 TeV) and $-6.15$ (5.36 TeV), so any statement at the LHC rests on an extrapolation. Only at 2.76 TeV do the data come close to the window $-2.5 \le\eta' \le0$ of the test below: at the most forward point, $\eta' = -2.61$, just outside it, the ratio to PHOBOS at 200 GeV, interpolated linearly in $\eta'$, is $1.09 \pm0.09$ (total uncertainties of both measurements in quadrature, without the $\langle N_{\mathrm{part}}\rangle$ uncertainties).

The data alone do not fix the extrapolation

Dip-free $dN_{\mathrm{ch}}/dy$ shapes that describe the 2.76 TeV data equally well give $(2/\langle N_{\mathrm{part}}\rangle)\,dN_{\mathrm{ch}}/d\eta= 0.64$–$0.76$ at $\eta' = 0$ for the three smooth families and 0.01 for the same families with a kinematic edge (medians of bootstrap fits, i.e. refits of pseudo-data drawn from the best fit with the data covariance, 100 per family). With the energy estimator $f_E$ of Sect. 7.2, the best fits of the smooth families give 1.4–1.5 (2.76 TeV) and 2.4–2.9 (5.02 TeV) times its RHIC calibration $\kappa$, and those of the edge families 0.6 times $\kappa$ at 2.76 TeV (Figure 7b, where $\kappa$ is the shaded band). A limiting-fragmentation conclusion drawn from such extrapolations reflects the assumed shape.

Energy sum rule

The estimated energy fraction

counts the produced charged hadrons with their neutral partners, $e_{\mathrm{prod}}(y)=\sum_i x_i(y)\langle m_T\rangle_i(y)\nu_i$, where $x_i$ and the mean transverse masses $\langle m_T\rangle_i$ follow from the kernel and $\nu_i$, the number of charged plus neutral hadrons per charged hadron, is $\nu_\pi=1.5$ and $\nu_K=\nu_p=\nu_{\mathrm{other}}=2$. The explicit net-baryon component $B(y)$ is normalized to $\langle N_{\mathrm{part}}\rangle$. Its profile and its mean transverse mass $\langle m_T\rangle_B$ come from the BRAHMS net-proton measurement at 200 GeV [27], with $c_B=2.03$ net baryons per net proton, and are carried to other energies at fixed rapidity loss, which BRAHMS found to be the same within uncertainties at 62.4 and 200 GeV [27, 42]. This net-baryon profile gives a mean rapidity loss of 2.04 (BRAHMS: $2.04 \pm0.10$) and an energy per net baryon of 27.7–33.2 GeV (BRAHMS: 30 $\pm$ 2 and 26 $\pm$ 5 GeV from two parametrizations) [27]. Besides this nominal rule, five variants are used: a flat rule ($\nu=1.5$ for all species and no net-baryon component), a species rule (the factors $\nu_i$, no net-baryon component), the net-baryon profile shifted away from beam rapidity by $-0.59$ and $+0.41$ units, toward the ends of the conservative range 1.45–2.45 of the rapidity loss quoted by BRAHMS [27] (mean loss 1.65 and 2.34 at 200 GeV after the cut at beam rapidity), and $\langle m_T\rangle_B$ falling to the proton mass at beam rapidity.

$f_E$ is an estimator, not the energy fraction itself: it depends on the composition and $\langle m_T\rangle$ near beam rapidity, where neither is measured. It is therefore calibrated where the data cover nearly the whole rapidity range and the constant-rapidity-loss profile applies, at 130 and 200 GeV; 62.4 and 19.6 GeV serve as checks. There the three dip-free families with a kinematic edge, the shapes used at the LHC, are fitted within $|\eta| \le y_{\rm beam}$ as in the dip test, and $f_E$ is integrated up to the most forward fitted point; at the LHC it is integrated up to beam rapidity. $\kappa$ is their mean with the nominal kernel. $\sigma_\kappa$ combines the larger of their standard deviation and half the largest deviation from $\kappa$ among the kernel’s rapidity-dependence variants, at least half the difference to the smooth families (which give 1.22), with the PHOBOS scale uncertainty: 5.3% of $\kappa-c$, from the normalization share of the systematic uncertainty and the $\langle N_{\rm part}\rangle$ uncertainty. The same scale also enters the reference window of Sect. 7.3; this correlation is neglected, which widens the intervals in $s$. This gives $\kappa= 1.19 \pm0.09$ for the nominal rule and $0.85 \pm0.09$ for the flat rule (Figure 7a). Without the net-baryon component (species rule) the calibration gives $0.96 \pm0.10$: the produced hadrons with their neutral partners account for about the participant energy, and the net-baryon component adds $c = 0.22$. Of the species-rule value, 47–52% comes from beyond the last BRAHMS point, where the composition and $\langle m_T\rangle$ are held at their last measured values; there the kernel’s proton yield (10.0 per unit $y$ at $y = 4$, 200 GeV) is below the net-proton density of the model (29.4). Why $\kappa$ exceeds unity is therefore not identified: the calibration absorbs it, and only its transfer to the LHC is taken from the generators. As checks, at 62.4 and 19.6 GeV, fitted and integrated in the same way, the nominal rule gives $f_E = 1.26$ and 1.25 (+0.9 and +0.7 $\sigma_\kappa$ from $\kappa$) and the flat rule 0.91 and 0.85 (+0.6 and −0.1 $\sigma_\kappa$), although near SPS energies the rapidity loss still grows with beam rapidity [27]. The generators, evaluated in the same way, give $R = 1.09$ at 19.6 GeV against 0.99 at 200 GeV in HIJING and 1.22 against 1.11 in EPOS4 (Table 2): their estimator bias changes more with energy than that of the data. Only the transfer of the calibration from RHIC to the LHC is therefore taken from the generators, and from both of them, below.

At the LHC the fits are repeated with the term $[(f_E-\kappa)/\sigma_\kappa]^2$ (the squared “pull”) added to $\chi^2$. Smooth shapes then describe the data poorly, and a separate fall-off below beam rapidity is needed, which the edge families of Sect. 3 provide: they describe both the data and the sum rule (Figure 7b). In the best fits of the test below (Table 3) the edge lies close to beam rapidity ($y_e = 7.4$–8.7, $a_e = 0.28$–0.81), except for box+edge and gG+edge at 5.02 TeV. There the central shape itself falls off sharply near beam rapidity, and the edge factor ($y_e = 4.5$, $a_e = 1.38$–1.39) acts as a broad suppression of the forward region.

$\sqrt{s_{NN}}$	family	best $s$ (95% interval)	$\Delta\chi^2(s=1)$	pull	best $s$	$\Delta\chi^2(s=1)$	best $s$	$\Delta\chi^2(s=1)$	best $s$
2.76 TeV	dG+edge	1.09 (0.94–1.28)	1.3	2.5	1.10	1.5	1.11	2.0	1.10
2.76 TeV	box+edge	1.09 (0.94–1.27)	1.3	2.5	1.09	1.5	1.10	1.9	1.10
2.76 TeV	gG+edge	1.09 (0.94–1.29)	1.3	2.5	1.10	1.6	1.12	2.3	1.11
5.02 TeV	dG+edge	1.03 (0.85–1.27)	0.1	2.9	1.06	0.4	1.11	1.3	1.07
5.02 TeV	box+edge	1.00 (0.83–1.23)	0.0	2.6	1.03	0.0	1.08	0.7	1.04
5.02 TeV	gG+edge	1.00 (0.83–1.23)	0.0	2.6	1.03	0.1	1.09	0.8	1.04
5.36 TeV	dG+edge	0.92 (0.77–1.12)	0.8	1.7	0.95	0.3	1.02	0.0	0.96
5.36 TeV	box+edge	0.84 (0.68–1.04)	2.6	1.2	0.88	1.4	0.96	0.1	0.90
5.36 TeV	gG+edge	0.90 (0.70–1.10)	1.1	1.6	0.94	0.4	1.00	0.0	0.95

Test of limiting fragmentation

The sum rule constrains the integral, not the shape, of the fragmentation region: the extrapolation itself falls below the PHOBOS data near $\eta' = 0$ (Figure 8). The RHIC shape is therefore imposed, and only its normalization is tested. The PHOBOS 200 GeV distribution $(2/\langle N_{\rm part}\rangle)\,dN_{\rm ch}/d\eta'$ with $\langle N_{\rm part}\rangle= 344 \pm11$, scaled by $s$, is added as pseudo-data at its own points in $-2.5 \le\eta' \le0$. Its covariance is the nominal RHIC covariance of these points plus the $\langle N_{\rm part}\rangle$ uncertainty as a common normalization; for $s$ times this distribution the uncertainties scale with $s$ and the covariance with $s^2$. This is equivalent to comparing the LHC prediction divided by $s$ with the PHOBOS points and their fixed covariance, so $\chi^2_{\rm pseudo}$ needs no determinant term. Only the nominal covariances of the LHC data and of the PHOBOS points are used here; the covariance variants of the dip test are not propagated. The total $\chi^2_{\rm total} = \chi^2_{\rm data} + \text{pull}^2 + \chi^2_{\rm pseudo}$ is minimized at each $s$ on a grid of step 0.02 for the three edge families, and limiting fragmentation corresponds to $s = 1$. Each scan point is refitted from its neighbors’ solutions until no point improves. Where two solution branches exist, the scan follows the lower one; the switch between them can produce small kinks in the scan curves. The intervals use the asymptotic thresholds $\Delta\chi^2 = 1$ and 3.84; their coverage was not checked with pseudo-experiments.

Results and the transfer of the calibration

Table 3 gives the results. At 2.76 TeV the ratio is $s = 1.09$ with $\Delta\chi^2(s=1) = 1.3$; the direct ratio at the most forward data point, $1.09 \pm0.09$, is consistent with it. At 5.02 TeV $s = 1.00$–1.03, with its minimum at $s \approx1$. At 5.36 TeV, where CMS covers only $|\eta| \le2.5$ and the result is an extrapolation by the sum rule, $s = 0.84$–0.92 with $\Delta\chi^2(s=1) = 0.8$–2.6. The 95% intervals, taken over the families, are wide, for example 0.94–1.29 at 2.76 TeV and 0.83–1.27 at 5.02 TeV, because the RHIC uncertainties in the window are large and partly correlated. The difference between 5.02 and 5.36 TeV, close in energy, is within these intervals and indicates the size of the extrapolation uncertainty at 5.36 TeV, where the data end furthest from the window. Over the six sum rules (coarse scan, step 0.1) the best $s$ changes by less than one grid step (the same grid point for every rule). The coarse grid underestimates $\Delta\chi^2(s=1)$ slightly; on it the largest value is 1.4 at 2.76 TeV and 2.6 at 5.36 TeV (flat rule, box+edge), which stays below the 95% threshold of 3.84 for every rule and family. The sum-rule pull at the best $s$ is 2.5 at 2.76 TeV, 2.6–2.9 at 5.02 TeV and

1.2–1.7 at 5.36 TeV: the LHC fits prefer more estimated energy than the RHIC calibration allows.

The covariances matter more. The scan was repeated on a grid of step 0.05, each point restarted from the nominal solutions, for the documented variants of the LHC covariances (and, at 5.02 TeV, the forward ALICE uncertainties correlated as at 2.76 TeV), for three settings of the covariance scan of Sect. 3 for the PHOBOS points in the window, $(f,u,\ell) = (0,1,-)$, $(0.8,0.2,1.5)$ and $(0.5,0.2,1.5)$, and for the kernel’s rapidity dependence at the LHC (none, in absolute $y$, and the BRAHMS dependence at a fixed distance from beam rapidity).

Over these variants the best $s$ is 1.00–1.23 at 2.76 TeV, 0.91–1.23 at 5.02 TeV and 0.75–1.14 at 5.36 TeV, and $s = 1$ leaves the 95% interval at 2.76 TeV, with $s > 1$, for all three families under strongly correlated PHOBOS uncertainties; at 5.02 TeV, with $s > 1$, for one family under strongly correlated PHOBOS uncertainties ($\Delta\chi^2(s=1) = 4.1$); at 5.36 TeV, with $s < 1$, for all three families under uncorrelated PHOBOS uncertainties, for two families under a kernel without rapidity dependence and for one family under the BRAHMS dependence anchored at beam rapidity ($\Delta\chi^2(s=1) = 4.0$). No single setting excludes $s = 1$ at all three energies, and interval over the families excludes it only at 2.76 TeV ($s>1$) for strongly correlated PHOBOS uncertainties and at 5.36 TeV ($s<1$) for uncorrelated PHOBOS uncertainties. The unpublished correlations of the PHOBOS uncertainties in the window move the best $s$ by up to 0.22, comparable to the half-width of the 95% interval. In 100 pseudo-experiments generated at $s=1$ for box+edge at 2.76 and 5.36 TeV, $\Delta\chi^2(s=1)$ exceeds 3.84 in a fraction 0.04–0.06 of them, consistent with the nominal 0.05 (at 5.36 TeV the minimum lies at the edge of the pseudo-experiments’ scan in up to 27 of them, for which $\Delta\chi^2(s=1)$ is a lower bound). Taking the PHOBOS uncertainties as quoted, at 90% confidence level, instead of dividing them by 1.645 (Sect. 2) leaves $\kappa$ unchanged, because the PHOBOS covariance is scaled as a whole; it raises $\sigma_\kappa$ to 0.10 through the scale uncertainty (7.6% instead of 5.3%) and widens the uncertainties of the window. The best $s$ then moves by at most 0.07, the 95% intervals over the families widen to 0.89–1.41, 0.79–1.54 and 0.63–1.23, $\Delta\chi^2(s=1)\leq1.4$, and the pull at 2.76 and 5.02 TeV falls to 1.6–1.9 $\sigma_\kappa$: the conversion sets the precision of the test, not its central values.

The calibration assumes that the estimator’s bias, from the composition and $\langle m_T\rangle$ near beam rapidity where neither is measured, is the same at RHIC and at the LHC. This is tested in both generators with the participant energy as truth (Table 2), with the kernel built as for the data: at the LHC the generator’s 200 GeV rapidity dependence is carried in $y/y_{\mathrm{beam}}$. The nominal rule gives $R=0.99$ (HIJING) and 1.11 (EPOS4) at 200 GeV within the measured range, and 1.04–1.05 and 1.28–1.31 at the LHC over $|y|\leq y_{\mathrm{beam}}$: the calibration transfers with a factor 1.05–1.06 in HIJING and 1.16–1.18 in EPOS4. With each generator’s own rapidity dependence at the LHC the factors would be 1.05–1.07 and 1.03, so the difference between the generators arises from how the rapidity dependence of the spectra changes with energy, which current data do not constrain. The absolute $R$ depends on the generator; what is used here is the transfer. With $\kappa\times1.06$ (HIJING), $s = 1.09$–$1.10$ at 2.76 TeV, $1.03$–$1.06$ at 5.02 TeV and $0.88$–$0.95$ at 5.36 TeV, with $\Delta\chi^2(s=1)$ at most 1.6, 0.4 and 1.4. With $\kappa\times1.17$ (EPOS4) the values are $s = 1.10$–$1.12$, $1.08$–$1.11$ and $0.96$–$1.02$, with $\Delta\chi^2(s=1)$ at most 2.3, 1.3 and 0.1. The calibration with the STAR-only spectra, $\kappa_{\mathrm{STAR}} = \kappa\times1.08$, gives nearly the same $s$ as the HIJING transfer (Table 3). At 2.76 TeV, where the data reach the edge of the window, the result depends only weakly on the calibration.

Under these assumptions limiting fragmentation is consistent with the data and the energy sum rule at all three LHC energies: $s = 1$ lies within the 95% interval of every family for the nominal rule and the nominal covariances, with or without either calibration transfer. Three caveats qualify this. With other plausible correlations of the PHOBOS uncertainties, or another rapidity dependence of the kernel, $s = 1$ lies outside the 95% interval of at least one family at each energy, although no single setting excludes it at all three energies and the interval over the families excludes it only at 2.76 TeV ($s > 1$) for strongly correlated PHOBOS uncertainties and at 5.36 TeV ($s < 1$) for uncorrelated PHOBOS uncertainties (above). At 2.76 and 5.02 TeV the LHC fits prefer 2.5–2.9 $\sigma_\kappa$ more estimated energy than the calibration, and at 5.36 TeV 1.2–1.7 $\sigma_\kappa$ (the pull above); the transfers indicated by the generators reduce this tension, to a pull of 2.0–2.5 (HIJING) and 1.1–1.9 (EPOS4) at 2.76 and 5.02 TeV. At 5.36 TeV a calibration lower by the same factor as the HIJING transfer, $\kappa\times0.94$, would place $s = 1$ outside the 95% interval for at least one of the three families ($\Delta\chi^2(s=1) \ge4.2$). All statements are conditional on the RHIC shape of the fragmentation region (Sect. 8); the data do not show that the LHC follows it.

Discussion

Relation to source models. The result is consistent with the RDM and diffusion-model analyses that attribute the dip mainly to the Jacobian [9, 10, 11]. It also suggests why a calculation that took the pion mass and an assumed $\langle p_T\rangle$ for its effective Jacobian [8] reached a different conclusion from later ones whose effective Jacobian was fixed by identified-particle data [9, 10]: whether a dip-free $\mathrm{d}N_{\mathrm{ch}}/\mathrm{d}y$ reproduces the $\eta$-dip depends on the flatness of $\mathrm{d}N_{\mathrm{ch}}/\mathrm{d}y$ near $y = 0$ and on the full species and $p_T$ dependence of the Jacobian. The test does not distinguish between production mechanisms that yield a dip-free $\mathrm{d}N_{\mathrm{ch}}/\mathrm{d}y$ of the allowed shapes, such as a dominant midrapidity source with fragmentation sources.

In p+Pb collisions the RDM finds the Jacobian insufficient to explain the dip [43]; the test does not address asymmetric systems. The limiting-fragmentation results are consistent with the RDM expectation that limiting fragmentation also holds at LHC energies [9, 10]; here the extrapolation is constrained by an energy sum rule calibrated on RHIC data, with the RHIC shape of the fragmentation region imposed, rather than by a boundary condition.

What is not shown. Broad depressions in $\mathrm{d}N_{\mathrm{ch}}/\mathrm{d}y$ are only weakly constrained, because the null families can absorb them, and the test is one-sided. A $\mathrm{d}N_{\mathrm{ch}}/\mathrm{d}y$ with a distinct forward component (HIJING) or a kernel that describes the forward region poorly (EPOS4) can produce a false dip, as both generators do at 200 GeV (Sect. 5). The split of the PHOBOS systematic uncertainties is not known and is scanned; the dip-test conclusion holds over the whole scan.

Assumptions of the limiting-fragmentation test. The test determines the normalization of the LHC fragmentation region relative to RHIC; it assumes that the shape is the same as at RHIC. A change of shape at fixed normalization would not be detected, and the data do not show that the LHC follows the RHIC shape. Its result depends on the correlations of the PHOBOS uncertainties in the window, which are not published (Sect. 7). It further assumes that the bias of the energy estimator transfers from RHIC to the LHC; the two generators give factors of 1.06 and 1.17, and the spread between them is a model uncertainty that current data cannot reduce. It also assumes the net-baryon profile at constant rapidity loss and the kinematic-edge form of the fall-off. The absolute value of $f_E$ differs between the data ($\kappa$) and the generators ($R$ in Table 2). The difference is not understood; much of $f_E$ comes from beyond the measured BRAHMS range, where neither composition nor spectra are measured (Sect. 7.2), and only the transfer is taken from the generators.

Spectra at RHIC. At 62.4 and 200 GeV the STAR spectra cover a narrow $p_T$ range, and their Lévy–Tsallis extrapolation inflates $\langle p_T\rangle$. At 200 GeV it gives a proton $\langle p_T\rangle$ of 1.42 GeV, above the 1.35 GeV of the 2.76 TeV spectra, whereas the PHENIX spectra give 1.01 GeV; at 62.4 GeV the extended spectra give 0.95 GeV instead of 1.30 GeV. EPOS4, which has radial flow, shows a bias of the same sign (protons +25% at 200 GeV, against +40% in the data), and HIJING a much smaller one (Sect. 5). The nominal kernel therefore uses the PHENIX spectra at 200 GeV and the extended spectra at 62.4 GeV, with the STAR-only spectra as a variant. The dip result does not depend on this choice. The energy calibration does, with $\kappa_{\mathrm{STAR}} = 1.28 \pm0.15$ for the STAR-only spectra, although the limiting-fragmentation ratio changes little (Table 3).

Outlook. The precision of the limiting-fragmentation test is limited by the RHIC uncertainties in the test window and their unknown correlations, which a published breakdown of the PHOBOS systematic uncertainties would constrain, and, at 5.02 and 5.36 TeV, by the coverage of the LHC data. Measurements of $dN_{\mathrm{ch}}/d\eta$ closer to beam rapidity at the LHC would test the extrapolation directly, and identified spectra over a wider rapidity range at RHIC would constrain both the kernel and the energy calibration.

Summary

Without an event generator, we have tested whether the midrapidity dip of $dN_{\mathrm{ch}}/d\eta$ in central heavy-ion collisions requires a structure in $dN_{\mathrm{ch}}/dy$, using a kernel built from measured identified spectra, and have quantified the sensitivity of the test. At all seven energies from 19.6 GeV to 5.36 TeV the data do not require it: they are consistent with a dip-free $dN_{\mathrm{ch}}/dy$ folded with the $y \to\eta$ kinematics of the measured spectra. For dips with widths up to one unit of rapidity, the depth giving $p < 0.05$ in half of the pseudo-experiments is 2.9–8.2% at the LHC and 5.9–12.5% at 62.4–200 GeV; broad depressions are only weakly constrained. Closure tests in HIJING and EPOS4 probe the kernel and the procedure and show a failure mode, a false dip at top RHIC energy that disappears when the most forward points are excluded from the fit.

For limiting fragmentation at the LHC the data alone leave the fragmentation region largely unconstrained. An energy sum rule calibrated on RHIC data, with the RHIC shape of the fragmentation region imposed, gives best-fit LHC/RHIC ratios of 1.09 at 2.76 TeV, 1.00–1.03 at 5.02 TeV and 0.84–0.92 at 5.36 TeV, with 95% intervals that include unity at all three energies, with or without the calibration transfers indicated by HIJING ($\kappa\times1.06$) and EPOS4 ($\kappa\times1.17$). Under these assumptions $s = 1$ is thus allowed, although at 2.76 and 5.02 TeV the LHC fits prefer 2.5–2.9 $\sigma_{\kappa}$ more estimated energy than the calibration. The result depends on the unpublished correlations of the PHOBOS uncertainties in the test window and on the kernel’s rapidity dependence: over these variants the best-fit ratios span 1.00–1.23, 0.91–1.23 and 0.75–1.14, and for some of them $s = 1$ lies outside the 95% interval of at least one family at each energy, although no single setting excludes it

excludes it only at 2.76 TeV ($s > 1$) for strongly correlated PHOBOS uncertainties and at 5.36 TeV ($s < 1$) for uncorrelated PHOBOS uncertainties. Limiting fragmentation is therefore neither established nor excluded by these data. At 2.76 TeV the result is consistent with a direct comparison at the most forward measured point; at 5.02 and 5.36 TeV it rests on the sum-rule extrapolation, which measurements closer to beam rapidity would test.

Acknowledgements The authors gratefully acknowledge Poojan Angiras, Manisha Rana and Sachin Rana for valuable discussions and feedback. The computations used the CERN batch computing service.

Data Availability Statement This manuscript has associated data in a data repository. [Authors’ comment: All experimental data used in this study are public. The LHC pseudorapidity densities, the identified-hadron spectra and the BRAHMS data are HEPData records; the PHOBOS tables are from the PHOBOS publication data archive, https://www.sdcc.bnl.gov/phobos/publications/Physics/.]

Code Availability Statement Code/software will be made available on reasonable request. [Authors’ comment: The analysis code and the scripts that produce every number, table and figure of this paper are available from the corresponding author on reasonable request.]

Competing interests The authors declare no competing interests.

Documented LHC covariances

Each LHC point’s systematic uncertainty $\sigma_j^{\mathrm{sys}}$ is split into a fully correlated normalization part and a remainder. For ALICE at 2.76 TeV the normalization part is 1% of the value; of the remaining variance a share of 80% is fully correlated within $|\eta| < 2$ and of 77% beyond, and the rest is uncorrelated; the two $|\eta|$ regions are correlated with each other only through the normalization. At 5.02 TeV the normalization part is 0.6%, with a correlated share of 80% within $|\eta| < 2$ and none beyond. For CMS at 5.36 TeV the normalization part is 0.1%, and the remaining variance is split into uncorrelated, exponential-shape (correlation length 1 unit of $\eta$) and fully correlated parts in the proportions 0.37, 0.37 and 0.26. The variants change the correlated share within $|\eta| < 2$ to 50% and 95%, the normalization at 2.76 TeV to 2%, and the CMS correlation length to 0.5 and 2. The diagonal always equals the published total.

Changes after the decision rule was fixed

The method and decision rule were fixed in the authors’ analysis notes before any $p$-value was computed. They were not registered publicly, so this ordering rests on the authors’ own records. Each later change was fixed before the corresponding rerun: (1) the covariance model was replaced by the three-way split of Sect. 3, because the original form was nearly singular; (2) the RHIC energies were added, with their fit ranges and variants; (3) a fitter defect that could leave H0 in a local minimum was fixed; in the data it had inflated $q$ (one spurious $p = 0.04$ at 19.6 GeV), in the pseudo-experiments it had biased $p$ upward; (4) digitized PHOBOS points were replaced by the official tables; (5) the symmetrized LHC tables are used with each $|\eta|$ once, and the documented LHC covariances became nominal; (6) the PHENIX spectra at 200 GeV and the extended spectra at 62.4 GeV became nominal, after the extrapolation of the narrow STAR $p_T$ ranges was found to inflate $\langle p_T \rangle$ (Sect. 8). None of these changes altered a decision. Two further departures from the original plan were made after the results: the injection grid was extended from three to seven widths, and the planned 95% confidence-level upper limits on $d$ from the profile likelihood are not quoted, because they are not reliable at the boundary $d \ge0$; $\mathrm{CL}_s$ limits from the injection pseudo-experiments are quoted instead (Sect. 4). The mirror-symmetric covariance, the injection into a truth outside the families and the rejection rate at the smallest depths (Sects. 3.5 and 4) were also added after the results. The generalized-Gaussian family was added as a robustness check after the plan was fixed; it does not enter the decision. The edge-family check of Sect. 4 and the limiting-fragmentation analysis are not part of the original plan.