Both hypothesis are rather vague as far as numerical predictions of hadronic yields are concerned. They will be confronted in the following sub-sections with the complete set of $\pi^-$ data treated in the preceding sections.

§.§ The Hypothesis of Limiting Fragmentation

The approach starts from the excitation of the target and projectile nucleons into (unspecified) masses $M > m_p$ and their subsequent decay

\begin{equation}
\label{eq:mdec}
M \rightarrow n + 1
\end{equation}

where $n$ concerns secondaries and "1" specifies a final state nucleon thus ensuring baryon number conservation. More specifically the following decays of $M$ are considered:

\begin{equation}
\label{eq:mdecmult}
M \rightarrow p+\pi,\;p+2\pi,\;p+3\pi,\;...,\;p+n\pi
\end{equation}

In this sense the hypothesis is similar to preceding ideas like the "fireball" [82], "isobar" [83] or "diffraction dissociation" [84] models. It is specific in the sense that it makes a precise prediction as to the behaviour of particle densities in the target and projectile laboratory systems in the high energy limit. For this purpose a connection to the scaling limit of inelastic electron-proton scattering is established where it is shown that the excitation of target or projectile into a mass $M$ approaches a limit at high energy. This connection assumes of course the independence of target nucleon excitation on the type of projectile particle ("factorization"), one of the many conjectures contained in [81].

In the target rest frame with the coordinates $p_T$ and $p_L^{\textrm{lab}}$, the kinematic limit is described by elastic scattering defining a parabola in $p_T$ and $p_L^{\textrm{lab}}$ centred at $p_L^{\textrm{lab}}$ = 0 GeV/c for protons and negative $p_L^{\textrm{lab}}$ for lighter secondaries. For pions the region of interest is defined  by $p_L^{\textrm{lab}} <$ 0 GeV/c [81]. At the time of publication, only very scarce experimental data in this area (which corresponds to the very forward region in the experiments) were available, essentially for the two beam momenta of 19 and 30 GeV/c. This led to the statement that the expected limit was already reached at the higher beam momentum, a rather bold claim also in view of the normalization problems of the data (see Sect. <ref>).

Given the global data interpolation elaborated above over the wide range of interaction energies from 3 $< \sqrt{s} <$ 63 GeV or 1 $< \log{s} <$ 3.6, a new sensible test of the hypothesis of limiting fragmentation may be performed. As at the time of publication [81] the interaction cross sections were assumed to be constant, the known energy dependence of the inelastic cross section represents an additional challenge to any scaling approach, see Sect. <ref>. In the following Figs. <ref>–<ref> therefore both the quantities $f/\sigma_{\textrm{inel}}$ and $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$ are presented although the variation of $\sigma_{\textrm{inel}}$ over the given $s$ range is only about 20%.

$\pi^-$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.1 GeV/c and $\log(s)$ from 1 to 3.6, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

$\pi^-$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.2 GeV/c and $\log(s)$ from 1 to 3.6, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

$\pi^-$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.3 GeV/c and $\log(s)$ from 1 to 3.6, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

$\pi^-$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.4 GeV/c and $\log(s)$ from 1 to 3.6, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

$\pi^-$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.5 GeV/c and $\log(s)$ from 1 to 3.6, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

The closeness of the $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$ distributions over the full available range of interaction energies is striking (panels a). This provides strong evidence for the hypothesis of Limiting Fragmentation. The normalized cross sections $f/\sigma_{\textrm{inel}}$ (panels b) on the other hand spread considerably more with $\log(s)$. This indicates a preference for constant invariant cross sections in this specific area of phase space rather than constant hadronization over the full increasing interaction area.

A further, more stringent test may be provided by the $\pi^0$ data of LHCf (Sect. <ref>). Here the predicted $\pi^0$ cross sections (<ref>), (<ref>) at the NA49 and ISR energies can be compared to the LHCf data. This corresponds to a range of more than 2 orders of magnitude in $\sqrt{s}$. The results are presented in Figs. <ref>–<ref>.

$\pi^0$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.1 GeV/c for NA49 and ISR (predictions) and LHC energies, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

$\pi^0$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.2 GeV/c for NA49 and ISR (predictions) and LHC energies, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

$\pi^0$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.3 GeV/c for NA49 and ISR (predictions) and LHC energies, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

$\pi^0$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.4 GeV/c for NA49 and ISR (predictions) and LHC energies, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

$\pi^0$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.5 GeV/c for NA49 and ISR (predictions) and LHC energies, a) $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$, b) $f/\sigma_{\textrm{inel}}$

As expected from the discussion of the $p_T$ distributions for different large rapidities in Sect. <ref>, the $d^2\sigma/(dp_L^{\textrm{lab}}dp_T)$ distributions are surprisingly $s$-independent for $p_T \gtrsim$ 0.2 GeV/c up to 7 TeV and for all $p_T$ at 2.76 TeV. The convergence with increasing $\sqrt{s}$ improves with decreasing $p_L^{\textrm{lab}}$ and $p_T$, see Fig. <ref>. On the other hand, the $f/\sigma_{\textrm{inel}}$ distributions split very considerably in these regions.

At this point some general remarks are in place. Pions are in fact rather special items to be regarded at negative $p_L^{\textrm{lab}}$ as it is rather difficult to transport them there. Their yields are orders of magnitude below the proton ones at low $p_L^{\textrm{lab}}$. Indeed a widespread belief would have it that pions from resonance decay ((<ref>) and (<ref>)) are all centred at low $x_F$ and low $p_T$. This is true for $\Lambda$ decay, Figs. <ref> and <ref> where the pion yields are effectively cut off at $p_T \sim$ 0.3 GeV/c and $x_F \sim$ 0.3. The $\Lambda$ is a low-$Q$ resonance in its p+$\pi$ decay with a mass distribution of zero width. In contrast, already the K$^0_S$ decay with its modest $Q$ of 200 MeV and its symmetric 2$\pi$ decay yields pions at large $x_F$ and considerable $p_T$, see Fig. <ref>. For strong decays, the wide mass distributions with the very long Breit-Wigner tails deconfine the decay pions to cover the full phase space. This is demonstrated in the Sect. <ref> dealing with resonance decay contributions. Another point is the decay chain p+n$\pi$ (<ref>). In a two body decay the center of the pion decay ellipse is placed with respect to the protons at

\begin{equation}
\label{eq:pi2p}
\langle x_F^\pi \rangle \sim \frac{m_\pi}{m_p} x_F(p)
\end{equation}

such that a diffractively produced baryon resonance yields pions with an average $x_F \sim$ 0.15 only. This value moves to smaller $x_F$ for multi-pion decays.

As a consequence several general remarks may be made as to the production of pions at negative $p_L^{\textrm{lab}}$ concerning:

*  Two body decays of high-$Q$ resonances.

*  In the decay p+$\pi$ backward decay of the proton in the overall cms.

* Strong resonance decays with important relativistic Breit-Wigner tails.

* Decays which are symmetric in the decay particle masses.

* Contribution of both baryonic and mesonic resonances.

* With increasing interaction energy the rate of low multiplicity decays – especially two-body decays – is strongly reduced thus de-populating the region of negative $p_L^{\textrm{lab}}$.

The authors of [81] have clearly pointed out points <ref> and <ref> claiming that at higher $n$ (<ref>) the pions will move towards the central rapidity area. As far as resonance mass distributions are concerned they see only a modest "widening of the mesa-like structures" in phase space corresponding to the decay process. Mesonic resonances are not foreseen, nor are decays with symmetric decay particle masses.

§.§ The hypothesis of Feynman scaling

In his original suggestion [1] Feynman took as well as [81] reference to deep inelastic lepton-proton scattering, however in a completely different context. In fact the existence of scaling longitudinal "structure functions" of proton constituents in this process had just been established [85] with $x_{Bj}$ (Bjorken $x$) as characteristic variable giving the momentum fraction of a parton relative to the proton momentum.

The conjecture of Feynman suggested a similar behaviour of secondary hadrons in p+p interactions if plotted in the fractional longitudinal momentum variable $x_F = 2p_L/\sqrt{s}$ (<ref>) thus creating an analogy with the partonic structure of the proton. In 1972 Berman, Bjorken and Kogut [86] predicted in fact partonic scattering in large momentum transfer hadronic interactions, first with electromagnetic coupling but already foreseeing strong coupling by "gluons" with correspondingly much higher cross sections.

The rapid evolution of these ideas into a full-blown theory of the strong interaction (Quantum Chromodynamics, QCD) cannot be followed here. Unfortunately the application of QCD is limited, due to the variable strong coupling constant $\alpha_s$, to the perturbative sector of QCD. This was followed up by Field and Feynman in 1977 [87] for e+e collisions and by Feynman, Field and Fox in 1978 [88] for the production of hadrons at high $p_T$ in hadronic interactions. A crucial problem in all applications of parton dynamics is the fact that for the prediction of final state hadronic yields one is invariably confronted with non-perturbative QCD in the final stage of parton hadronization.

For soft hadronic interactions as they are studied in this paper there are basically two approaches. Either higher-order QCD is pushed to the limit ("next-to-next-to leading order") in order to make predictions at lower and lower $p_T$. Today it is believed that this is possible down to about 1 GeV/c. Or one assumes colour exchange as the source of hadronic interactions which leads to the breakup of both participating hadrons into quarks and, in the case of baryons, diquarks thus connecting the target and projectile systems by "strings" fragmenting into the final state secondaries.

This paper does not attempt to refer to any of the practically unlimited number of "string" models which are today available. These are rather characterized by a large number of adjustable parameters and very limited predictive power. If for the quark end of the strings experimental data from leptonic interactions may be invoked, the fictitious diquark systems which have by the way to ensure baryon number conservation have no referable analogue in different sectors of particle physics. It is a question why the baryon number should be contained in a partonic subsystem of the nucleon.

In addition, the point-like interactions assumed to describe the hadronic collisions do not carry any connection to the fact that hadrons are extended objects where the impact parameter is a decisive variable for any collision. This is apparent in the sector of peripheral interactions [6] ("diffraction") which is not really describable by string models.

In the following the energy dependence of longitudinal $\pi^-$ distributions will be shown.

§.§.§ $y_{\textrm{lab}}$ distributions

It should be remembered that $y_{\textrm{lab}}$ is equivalent to $x_F$ (Fig. <ref>) in the fragmentation zone for $\sqrt{s} \gtrsim$ 10 GeV and $x_F \gtrsim$ 0.2. The corresponding cross sections are readily available from the global interpolation, Sect. <ref>. Both the invariant cross sections $f(y_{\textrm{lab}},p_T)$ and the yields per inelastic event, $f(y_{\textrm{lab}},p_T)/\sigma_{\textrm{inel}}$ will be shown.

$\pi^-$ cross sections as a function of $y_{\textrm{lab}}$ for $p_T$ = 0.05, 0.1, 0.2 GeV/c and different $\log(s)$ values in the range 1 $< \log{s} <$ 3.6; a) $f(y_{\textrm{lab}},p_T)$, b) $f(y_{\textrm{lab}},p_T)/\sigma_{\textrm{inel}}$. The distributions for different $p_T$ values are successively scaled down by factor of 1/3 for better separation

$\pi^-$ cross sections as a function of $y_{\textrm{lab}}$ for $p_T$ = 0.4, 0.6, 0.8 GeV/c and different $\log(s)$ values in the range 1 $< \log{s} <$ 3.6; a) $f(y_{\textrm{lab}},p_T)$, b) $f(y_{\textrm{lab}},p_T)/\sigma_{\textrm{inel}}$. The distributions for different $p_T$ values are successively scaled down by factor of 1/3 for better separation

$\pi^-$ cross sections as a function of $y_{\textrm{lab}}$ for $p_T$ = 1.0, 1.2 GeV/c and different $\log(s)$ values in the range 1 $< \log{s} <$ 3.6; a) $f(y_{\textrm{lab}},p_T)$, b) $f(y_{\textrm{lab}},p_T)/\sigma_{\textrm{inel}}$. The distributions for different $p_T$ values are successively scaled down by factor of 1/3 for better separation

$\pi^0$ cross sections as a function of $p_L^{\textrm{lab}}$ for $p_T$ = 0.1, 0.2, 0.3, 0.4, 0.5 GeV/c for NA49 and ISR (predictions) and LHC energies; a) $f(y_{\textrm{lab}},p_T)$, b) $f(y_{\textrm{lab}},p_T)/\sigma_{\textrm{inel}}$. The distributions for different $p_T$ values are successively scaled down by factor of 1/3 for better separation

There is rather precise scaling over a wide range of $y_{\textrm{lab}}$ and $p_T$ values. In the higher $p_T$ range there are of course deviations which are just reflecting energy-momentum conservation at the lower interaction energies. And it is again plotting the invariant cross sections $f(y_{\textrm{lab}},p_T)$ rather than the yield per inelastic event which gives narrower $s$-dependences. Extending the $s$ range up to LHC energies the $\pi^0$ distributions (Sect. <ref>) show corresponding performance for the invariant cross sections as presented in Figs. <ref> – <ref> as a function of $y_{\textrm{lab}}$ both for $f(y_{\textrm{lab}},p_T)$ and for $f(y_{\textrm{lab}},p_T)/\sigma_{\textrm{inel}}$. The lowest $\log(s)$ value in each plot is given on the lowest-lying curve. Several points are noteworthy in this context:

* Again as in Sect. <ref> the invariant cross sections $f(y_{\textrm{lab}},p_T)$ show a more precise scaling behaviour than the yields per inelastic event $f(y_{\textrm{lab}},p_T)/\sigma_{\textrm{inel}}$.

* For medium $p_T$ values the spread is only on the few percent level.

* The distributions split up, depending on $p_T$, at different $y_{\textrm{lab}}$ values in the approach to the non-scaling central region (Sect. <ref>).

* With increasing $p_T$ the distributions at the lowest $log{s}$ values indicated in the Figures move down. This is an effect of energy-momentum conservation.

The corresponding distributions for the $\pi^0$ data from LHCf [80] are shown in comparison with the $\pi^0$ predictions from NA49 and ISR(63 GeV) in Fig. <ref> for the invariant cross sections $f(y_{\textrm{lab}},p_T)$ and $f(y_{\textrm{lab}},p_T)/\sigma_{\textrm{inel}}$. The approach to scaling over this very extended $\sqrt{s}$ range with increasing $p_T$ and decreasing $y_{\textrm{lab}}$ (Figs. <ref> and <ref>) is clearly evident.

§.§.§ $x_F$ and $x_F'$ distributions

The global interpolation in ($y_{\textrm{lab}}$, $p_T$) (Sect. <ref>) has been interpolated to constant $x_F$ values for a grid of 24 $x_F$ and the standard 26 $p_T$ values. As the distributions $f(x_F,p_T)/\sigma_{\textrm{inel}}$ are rather steep at low $x_F$ the following bins have been used:

$x_F$ =    0,       0.01,    0.02,    0.03,    0.04,    0.05,    0.075,    0.1,    0.125,    0.15,    0.2,    0.25,

0.3,     0.35,    0.4,     0.45,    0.5,     0.55,     0.6,    0.65,    0.7,    0.75,    0.8,    0.85

An overall impression of the $x_F$ distributions is given by Fig. <ref> for four $\log(s)$ values where the $x_F$ definition, see (<ref>):

\begin{equation}
\label{eq:xf1}
x_F = \frac{2p_L}{\sqrt{s}}
\end{equation}

is used. As $x_F$ is defined in the region 0 $< x_F <$ 1, this allows a unified overview over the complete phase space at each energy.

$f(x_F,p_T)$ as a function of $x_F$ for the four interaction energies $\log{s)}$ = 1, 1.8, 2.6, 3.4 or $\sqrt{s}$ = 3.2, 7.9, 20.0 and 50.1 GeV. The full lines give the global interpolation for $p_T$ values between 0.05 and 1.3 GeV/c

Three distinct zones are distinguishable:

* The central area at $x_F$ = 0 where an increase of the cross section by one order of magnitude over the range 3.6 $< \sqrt{s} <$ 50 GeV is visible.

* The forward or "fragmentation" area at $x_F \gtrsim$ 0.2 where the cross section is about $s$-invariant.

* An intermediate area at 0.05 $< x_F <$ 0.2.

The central and intermediate zones will be treated in the subsequent Sect. <ref>. Several features of the distributions shown in Fig. <ref> are noteworthy:

* The cross sections at $p_T <$ 0.2 GeV/c are hardly distinguishable in the intermediate and forward regions.

* At the lowest value at $\log(s)$ = 1.0 there are no data for $p_T >$ 0.75 GeV/c.

* At this energy the distributions are shrinking in their $x_F$ range with increasing $p_T$.

* For $p_T \gtrsim$ 1 GeV/c there is a clear change in the shape of the $p_T$ distributions.

A more detailed view of this situation is presented in Figs. <ref> to <ref> where the invariant cross sections $f(x_F,p_T)$ and $f(x_F,p_T)/\sigma_{\textrm{inel}}$ are shown as a function of $x_F$ for seven values of $\log(s)$ between 1 and 3.4. For the higher range of $\log(s)$ there is again a smaller spread of the yields for $f(x_F,p_T)$ than for $f(x_F,p_T)/\sigma_{\textrm{inel}}$. However, for $\sqrt{s}$ below about 12 GeV ($\log(s)\lesssim$ 1.8) there is an explicit non-scaling effect that reduces the yields especially at higher $p_T$ values. In addition, the $x_F$ range in the low-$s$ range is confined to successively smaller $x_F$ ranges with increasing $p_T$.

a) $f(x_F,p_T)$ as a function of $x_F$ for $p_T$ = 0.05, 0.1 and 0.2 GeV/c for the global interpolation at $\log(s)$ values from 1.0 to 3.6. The lines at subsequent $p_T$ values are scaled down by 1/3 for better separation. b) $f(x_F,p_T)/\sigma_{\textrm{inel}}$

a) $f(x_F,p_T)$ as a function of $x_F$ for $p_T$ = 0.4, 0.6 and 0.8 GeV/c for the global interpolation at $\log(s)$ values from 1.0 to 3.6. The lines at subsequent $p_T$ values are scaled down by 1/3 for better separation. b) $f(x_F,p_T)/\sigma_{\textrm{inel}}$

a) $f(x_F,p_T)$ as a function of $x_F$ for $p_T$ = 1.0 and 1.2 GeV/c for the global interpolation at $\log(s)$ values from 1.0 to 3.6. The lines at subsequent $p_T$ values are scaled down by 1/3 for better separation. b) $f(x_F,p_T)/\sigma_{\textrm{inel}}$

This has to do with the definition of $x_F$ (<ref>) being in principle only valid in an infinite momentum frame – as is the definition of $x_{Bj}$ for partons. For finite energies, the maximum final state momentum $p_L{\textrm{max}}$ is $s$-dependent and unequal to $\sqrt{s}$/2 as energy-momentum, charge and baryon number conservation have to be taken into account. This leads, for final state protons, to the value

\begin{equation}
p_L^{\textrm{max}} (p_T=0) = \frac{1}{2} \sqrt{s-4m_p^2} \: \neq \: \frac{1}{2} \sqrt{s}
\end{equation}

For $\pi^-$ production, there have to be, in addition to the $\pi^-$ mass, 2 nucleons for baryon number conservation and at least one $\pi^+$ for charge conservation in the final state,

\begin{equation}
p_L^{\textrm{max}} (p_T=0) = \frac{1}{2\sqrt{s}} \sqrt{s-(2m_{\pi}+2m_p)^2} \;\sqrt{s-4m_p^2}
\end{equation}

see [30] for further explanations. The definition of $x_F$ should  therefore be replaced by

\begin{equation}
x_F' = \frac{p_L}{p_L^{\textrm{max}}}
\end{equation}

see also (<ref>). This leads to sizeable deviations from definition (<ref>) at $\log(s) <$ 3 and subsequently to a dependence on $p_T$ as shown in Fig. <ref>.

The ratio $x_F'$/$x_F$ as a function of $\log(s)$ for different $p_T$ values

This traces the corresponding deviations from the correlation between $x_F$ and $y_{\textrm{lab}}$ shown in Fig. <ref>. A tolerable correctable deviation is only reached for the SPS energy range $\sqrt{s} \gtrsim$ 15 GeV. In fact the values $p_L^{\textrm{max}}(s,p_T)$ are corresponding to a scale $y_{\textrm{lab}}^{\textrm{min}}(s,p_T)$,

\begin{equation}
y_{\textrm{lab}}^{\textrm{min}} = y_{\textrm{beam}} - \frac{1}{2}ln\frac{ E + p_L^{\textrm{max}} }{ E - p_L^{\textrm{max}} }
\end{equation}

\begin{equation}
E = \sqrt{ (p_L^{\textrm{max}})^2 + p_T^2 + m_\pi^2}
\end{equation}

This is quantified in Fig. <ref> where the invariant cross section $f(y_{\textrm{lab}},p_T)$ is presented as a function of $y_{\textrm{lab}}$ for the relatively low interaction energy $\sqrt{s}$ = 4.0 GeV or $\log(s)$ = 1.2. For each $p_T$ the corresponding $y_{\textrm{lab}}^{\textrm{min}}$ value is indicated as a vertical line on the abscissa.

$f(y_{\textrm{lab}},p_T)$ as a function of $y_{\textrm{lab}}$ for $p_T$ values between 0.1 and 1.1 GeV/c in steps of 0.1 GeV/c and $\log(s)$ = 1.2. For each $p_T$ the corresponding $y_{\textrm{lab}}^{\textrm{min}}$ is indicated as a vertical line

The approach of the measured cross sections to the limiting $y_{\textrm{lab}}$ values as imposed by the conservation laws is clearly visible.

The invariant cross sections $f(x_F',p_T)$ are shown in Figs. <ref> to <ref> as a function of $x_F'$ for 8 values of $p_T$.

$f(x_F',p_T)$ as a function of $x_F'$ for $p_T$ = 0.05, 0.1 and 0.2 GeV/c for the global interpolation at $\log{s}$ values from 1.0 to 3.6. The plots at subsequent $p_T$ values are scaled down by 1/20 for better separation. The NA49 interpolation at $\log(s)$ = 2.48 is shown as a green line

$f(x_F',p_T)$ as a function of $x_F'$ for $p_T$ = 0.4, 0.6 and 0.8 GeV/c for the global interpolation at $\log(s)$ values from 1.0 to 3.6. The plots at subsequent $p_T$ values are scaled down by 1/5 for better separation. The NA49 interpolation at $\log(s)$ = 2.48 is shown as a green line

$f(x_F',p_T)$ as a function of $x_F'$ for $p_T$ = 1.0 and 1.2 GeV/c for the global interpolation at $\log(s)$ values from 1.0 to 3.6. The plots at subsequent $p_T$ values are scaled down by 1/5 for better separation. The NA49 interpolation at $\log(s)$ = 2.48 is shown as a green line

The $\pi^0$ data of LHCf [80] are presented in Fig. <ref> as a function of $x_F'$ for different $p_T$ values. As already evident from the $y_{\textrm{lab}}$ distributions (Fig. <ref>), scaling in $x_F'$ is approached for the higher $x_F'$ and $p_T$ regions if plotting the invariant cross sections $f(x_F',p_T)$.

Invariant $\pi^0$ cross sections $f(x_F',p_T)$ as a function of $x_F'$ for $p_T$ values between 0.1 and 0.5 GeV/c

For the lower $x_F'$ and $p_T$ regions the data exceed the NA49 and ISR results by factors of up to 1.7 corresponding to about 73% of the ratio of the inelastic cross sections. This indicates that the scaling data come from a constant "core" area of the increasing proton size with a transition region approaching a major fraction of the total interaction zone at lower $x_F'$ and $p_T$.

§.§ The non-scaling central region

The $x_F'$ distributions of Figs. <ref> to <ref> which are extended to the complete kinematical range clearly demonstrate a non-scaling region at low $x_F'$ with a sharply defined transition to approximate scaling located between $x_F \sim$ 0.16 at low $p_T$ approaching $x_F' \sim$ 0.3 in the higher $p_T$ range (always with the exception of the lowest $\log(s)$ values). This corresponds to a "rising rapidity plateau" [2, 3, 4] and has led to the assumption of a special source of central particle production ("pionization") and in consequence to an undue concentration of interest in the central rapidity area in contrast to the forward/backward production regions. This is specially flagrant for Heavy Ion collisions where a central "hot" zone is postulated as a "Quark-Gluon Plasma".

A detailed inspection of the low-$x_F'$ region allows however for a more differentiated picture especially if the $p_T$ dependence is taken into account. This is demonstrated in Figs. <ref> and <ref> where the evolution of the $s$-dependence with respect to the NA49 data is shown using the cross section ratio at $y$ = 0:

\begin{equation}
R(y=0,s,p_T) = \left(\frac{f}{\sigma_{\textrm{inel}}}\right)(s,p_T)\Bigg/\left(\frac{f}{\sigma_{\textrm{inel}}}\right)(NA49,p_T)
\label{eq:rat_y0}
\end{equation}

$R(y=0,s,p_T)$ is plotted as a function of $p_T$ for the $\log(s)$ range of the reference data in Fig. <ref>.

$R(y=0,s,p_T)$ as a function of $p_T$ for the reference data, $\log(s)$ from 1.0 to 3.6

Two main contributions to the evolution with increasing $s$ are clearly evident;

* An enhancement at low $p_T$ in the range $p_T <$ 0.3 GeV/c.

* A steady increase of yields up to the limit of $p_T$ at 1.3 GeV/c.

These trends continue up to LHC energies, Fig. <ref>, with a strong evolution of the higher $p_T$ region which exceeds the low $p_T$ enhancement above RHIC energies (notice the log scale).

$R(y=0,s,p_T)$ as a function of $p_T$ for $\sqrt{s}$ values from 3 GeV to 7 TeV

The approximately exponential increase of the cross section ratio above $p_T \sim$ 0.6 GeV/c may be used to separate the two components by extrapolating to low $p_T$ as shown in Figs. <ref> to <ref>.

$R(y=0,s,p_T)$ as a function of $p_T$ for $\log(s)$ = 7.7 with exponential fit to $p_T >$ 0.6 GeV/c, extrapolation to low $p_T$ (broken line) and subtracted ratio (dotted line)

$R(y=0,s,p_T)$ as a function of $p_T$ for $\log(s)$ = 6.9 with exponential fit to $p_T >$ 0.6 GeV/c, extrapolation to low $p_T$ (broken line) and subtracted ratio (dotted line)

$R(y=0,s,p_T)$ as a function of $p_T$ for $\log(s)$ = 3.6 with exponential fit to $p_T >$ 0.6 GeV/c, extrapolation to low $p_T$ (broken line) and subtracted ratio (dotted line)

The subtraction of the fit to higher $p_T$ yields exponential distributions with inverse slope parameters as shown in Fig. <ref>

Exponential fits to the cross section ratio (<ref>) as given in Figs. <ref> to <ref> a) inverse slopes for the low $p_T$ region $p_T <$ 0.5 GeV/c and b) slopes for the higher $p_T$ region at $p_T >$ 0.6 GeV/c

The extracted exponential behaviour of $R(y=0)$ in the low-$p_T$ region may be converted into invariant cross sections using the $p_T$ distribution from NA49 at $y$ = 0 with its inverse slope of 0.17 GeV/c, Fig. <ref>a, and finally into the $m_T-m$ distribution of Fig. <ref>b.

a) $f/\sigma_{\textrm{inel}}(y=0,s,p_T)$ as a function of $p_T$ for different values of $\log(s)$ compared to the invariant cross section of NA49; b) corresponding distribution as a function of $m_T-m$

The inverse slope or "temperature" of the $m_T$ distribution in the low-$p_T$ region is with 0.08 GeV/c about half the value of the inclusive data. In fact it is located about half-way between the "temperature" of the feed-down pions from K$^0_S$ and $\Lambda$ decay, Fig. <ref>. This gives a strong indication concerning the origin of pions in this part of the central rapidity region. In fact it is the decay of excited baryons [81] which gives preferentially low-$x_F$ pions due to the small mass ratio $m_\pi/m_p$ (<ref>). There are many channels open with small $Q$-values to pions like

\begin{equation}
\begin{split}
\Sigma(1385) &\rightarrow \Lambda + \pi \\
&\rightarrow \Sigma  + \pi
\end{split}
\end{equation}

(in a consistent picture these strange resonances ought to be subtracted from the inclusive pion sample like $\Lambda$, $\Sigma$ and K$^0$ decays. They are, however, as strong decays on-vertex).

\begin{equation}
\Xi(1320) \rightarrow \Lambda + \pi
\end{equation}

as well as heavy flavour baryons in the RHIC/LHC energy region. Concerning "normal" baryonic resonances there are many channels open with small $Q$-values like

\begin{equation}
\begin{split}
N^*(1440) &\rightarrow \Delta + \pi \\
N^*(1520) &\rightarrow p + \pi + \pi \\
&\rightarrow \Delta + \pi  \\
&\rightarrow p + \rho
\end{split}
\end{equation}

and of course higher resonances.

In this context it is interesting to look at comparable results for charged kaons [52] as shown in Fig. <ref> for the ratio $R(y=0,s,p_T)$ which are directly comparable to the $\pi^-$ data (Fig. <ref>).

$R(y=0,s,p_T)$ for charged kaons [52] as a function of $p_T$ for $\sqrt{s}$ values between 6.8 and 63 GeV

Also in this case the low-$p_T$ enhancement has been connected to low-$Q$ resonance decay, namely

\begin{equation}
\begin{split}
\Phi(1020)    &\rightarrow K^+ + K^- \\
\Lambda(1520) &\rightarrow p + K^-
\end{split}
\end{equation}

The quantitative results depend on the kaon charge due to the relatively large difference between the inclusive K$^+$ and K$^-$ yields in face of the charge symmetric $\Phi$ decay and the fact that $\Lambda(1520)$ exclusively feeds into K$^-$.

The energy dependence of the invariant cross section at low $p_T$ is shown in Fig. <ref> for $p_T$ = 0.1 GeV/c as a function of $\log(s)$.

$f/\sigma_{\textrm{inel}}$ at $p_T$ = 0.1 GeV/c as a function of $\log(s)$

The yield rises sharply from the threshold at $\log(s)$ = 0.67 or $\sqrt{s}$ = 2.2 GeV. It reaches an almost linear $\log(s)$ dependence above $\log(s)$ = 1.3 of the form

\begin{equation}
f/\sigma_{\textrm{inel}}(\log{s}) \sim -0.307 + 0.802\log{s}
\end{equation}

with small ondulating deviations between $\log{s} \sim$ 1.25 and 3.2. The yields at RHIC and LHC energy have negligible statistical errors but sizeable systematic uncertainties on the level of 5%.

For the higher $p_T$ region a more complex energy dependence is evident. This is presented for $p_T$ = 1.3 GeV/c in Fig. <ref>.

$f/\sigma_{\textrm{inel}}$ at $p_T$ = 1.3 GeV/c as a function of $\log(s)$

The invariant cross section again rises sharply from the threshold at $\log(s)$ = 1.2 or $\sqrt{s}$ = 4.0 GeV to a first plateau at $\sqrt{s} \sim$ 40 GeV with a following non-linear increase up to $\sqrt{s}$ = 13 TeV. It will be shown in Sect. <ref> that resonance decay governs the $\pi^-$ yields up to a transverse momentum of at least 3 GeV/c. In fact the total yields of hadronic resonances show a strongly non-linear energy dependence (see for instance Fig. <ref>) with a fast increase from threshold followed by a flattening off towards higher energies.

§.§ Pionization and factorization

The increase of central particle density in the ISR energy range came as a surprise [2, 3, 4], as cosmic ray data had indicated, on the contrary, a depletion. In consequence, a specific central production process ("pionization") was postulated as opposed to the independent hadronization of the target and projectile hemispheres with a central overlap. This independent fragmentation could also be described as "factorization" [81] in the sense that target (or projectile) particle densities would be independent of the type of projectile (or target) present in the collision. This would mean that experimentally for the inclusive reactions

\begin{align}
p + p                   &\rightarrow \pi^- + X \\
\langle \pi \rangle + p &\rightarrow \pi^- + X \\
\gamma + p              &\rightarrow \pi^- + X
\end{align}

the cross sections would be equal in the target hemisphere. The shape and extent of the central overlap region would be accessible via different experimental signatures such as long-range two-particle and multiplicity correlations, charge ratios or the extent of net baryon production. Also hadron-nucleus interactions may be used,

\begin{equation}
p + A \rightarrow \pi^- + X
\end{equation}

for the same purpose although the hadronization itself is of course different for both the target and projectile in this case.

A very detailed study of factorisation and overlap functions has been conducted in [5, 89]. Some examples are shown in Figs. <ref> to <ref> concerning the feedover-distributions for pions and for net protons. extracted from p+p, $\langle \pi \rangle$ + p and p + C collisions.

a) $\pi^+$/$\pi^-$ ratio in p+p and $\langle \pi \rangle$+p collisions; b) target and projectile components of charge-averaged pion production and c) target component with respect to the total yield of $\pi^+$ and $\pi^-$ in p+p interactions

$p_T$ integrated $\pi^+$/$\pi^-$ ratio as a function of $x_F$ in p+p and p+C: a) without $p_T$ cut off and b) with $p_T$ cut off included as the data points in p+p collisions are plotted only up to $x_F$ = 0.02. The data points concerned by an extrapolation of the low $p_T$ region are given as open symbols

The fact that the $\pi^+/\pi^-$ ratio has to be unity for the projectile contribution in $\langle \pi \rangle$ + p interactions is experimentally  verified in Fig. <ref>a such that only the target component will contribute to the measured $\pi^+/\pi^-$ ratio. The symmetric target and projectile distributions in p + p reactions, Fig. <ref>b, allow for the direct extraction of the target feed-over, Fig. <ref>c both in shape and range in $x_F$. A similar argument may be used in p + C interactions where the isoscalar Carbon nucleus has a $\pi^+/\pi^-$ ratio equal to unity (now in the target hemisphere) as presented in Fig. <ref>.
Note that this argument holds since there is no pion exchange at SPS energies as demonstrated in [90].

Concerning protons, $\langle \pi \rangle$ + p interactions should not have any net protons in the projectile hemisphere. In p + p collisions, a complementary method consists of fixing net baryon number by selecting protons either forward or backward at $x_F$ values outside the range of pair produced protons [91]. These methods allow for a completely model-independent way of measuring the feed-over distributions for net protons as shown in Fig. <ref>.

a) Target component for net protons measured with different methods and b) target components of $\langle \pi \rangle$ and net protons

This study shows that there is an important difference between the feed-over distributions of pions and net protons, Fig. <ref>b. In both cases the $x_F$ range of feed-over is well outside the accessible $x_F$ range of high energy colliders.

Note that also a test of factorization becomes possible between $\pi$ + p and p + p reactions as shown in Fig. <ref> in backward direction below $x_F \sim$ -0.5 where proton identification is possible in bubble chambers ($p_{\textrm{lab}}^{\textrm{p}} <$ 1.3 GeV/c). In addition, [53] has shown factorization for baryon production between p + p and leptoproduction at HERA [92, 93].

In the context of the present argumentation it is important to realize that it has been experimentally verified that no central "pionization" region is present in addition to the (factorizing) target and projectile contributions as already hypothesized in [81]. This is also true for p + C interactions.

§.§ Summary on scaling

The critical comparison of the complete set of existing data and the resulting establishment of a global interpolation scheme with a systematic uncertainty below the 5% level allow for a confrontation of certain physics hypotheses with the experimental reality over the full energy range available today. This has been attempted here with special emphasis on energy dependence or, rather, the absence thereof which is usually labelled as "scaling".

Several aspects are of importance here:

* The complete coverage of not only the energy range but also of a maximum of the available phase space is absolutely mandatory.

* Even on the single-particle inclusive level the complexity of the experimental results precludes any simplistic description of a physics situation which is not amenable to theoretical predictions.

* The increase of the total inelastic cross section by a factor of three over the available energy range imposes a very careful scrutiny of the geometrical aspects in terms of dependencies on the impact parameter and its consequences on different sectors of phase space, especially at LHC energies.

* In view of this fact any description in terms of point-like interactions as they are for instance imposed by parton dynamics is problematic.

Two hypotheses which were phrased rather early on in the development of particle physics (1969) have been confronted with the data in detail.

The first hypothesis concerns "limiting fragmentation" [81].

* This hypothesis was for the first time tested over the full energy range (Sect. <ref>).

* The predicted energy invariance is verified within tight limits if plotting the invariant cross section as a function of $p_L^{\textrm{lab}}$.

* On the other hand if using the cross section per inelastic event ($f/\sigma_{\textrm{inel}}$) sizeable deviations from $p_L^{\textrm{lab}}$ scaling are observed.

* This is a first indication of a yield dependence on the production point in the overlapping hadronic disks.

Ref. [81] is an astounding document as it contains a number of conjectures that were not testable at the time of publication. To be mentioned here is hadronic factorization and the absence of central pionization (Sect. <ref>). The inclusion of baryon  number conservation as a basic ingredient as well as resonance production and decay are further assets.

* The study of $y_{\textrm{lab}}$ distributions (Sect. <ref>) is complementary to $p_L^{\textrm{lab}}$ scaling and shows similar results, again in preference of invariant cross sections as compared to yields per inelastic event.

The second hypothesis, "Feynman scaling", is inspired by the partonic structure of the colliding hadrons as it is measured in deep inelastic lepton scattering.

* In consequence the conjecture postulates scaling in longitudinal momentum if referred to the maximum available momentum (Sect. <ref>).

* In the very definition of this relative variable it is necessary to take full account of energy and baryon number conservation.

* Nevertheless there are strong deviations from scaling at all energies below the ISR range.

* An inherent weakness of this hypothesis is the baryonic sector. In fact there are no informations to be gained on baryon distributions and baryon number conservation from deep inelastic lepton scattering. Hypothetical "di-quark" states have therefore to be introduced which carry unknown fragmentation functions.

* The sector of diffraction which is characterized by large impact parameters [6] is in no way describable in this approach.

In this context of impact-parameter dependencies it is worth while looking as a side remark at diffractive proton production which would stem [6, 7] from the expanding rim of the interaction region. One would therefore, contrary to the $\pi^-$ production, expect cross sections to scale in the yield per inelastic collision rather than the invariant cross section itself which should increase with the interaction area. This has been shown in [94] reproduced here as Fig. <ref>.

Ratio of proton cross sections $R$ relative to NA49 data as a function of $\sqrt{s}$ in the diffractive sector at $x_F$ = 0.92 and 0.95 including data from UA4 [94]. The full lines correspond to the invariant cross sections, the broken lines to the cross sections per inelastic collision

The UA4 experiment at the CERN p$\overline{\textrm{p}}$ collider published data [94] only at $\sqrt{s}$ = 530 GeV. It would be of considerable interest to have data also up to the highest LHC energies.

§ TRANSVERSE MOMENTUM

The available data concerning double differential inclusive $\pi^-$ cross sections present – given sufficient precision on the experimental
level – an overall phenomenology of great complexity defying attempts
at simple interpretations or algebraic approximations. It is however
the largely asymmetric occupation of the phase space in longitudinal
and transverse direction, the "longitudinal phase space" (<ref>), which seems
to offer a natural approach to some underlying physics. The longitudinal
part has been analysed in the preceding Sect. <ref> in terms of energy dependence and "scaling" with only marginal success. The transverse
part shows a steep $p_T$ dependence that has been from early on tentatively approximated by an exponential form

\begin{equation}
\frac{f}{\sigma( x_F', p_T, s )} = A( x_F', s )e^{-B( x_F', s ) p_T}
\end{equation}

Of course this simplistic form is only defendable as a zero-order approximation to the rather complex overall $p_T$ dependence that follows from the global interpolation (Sect. <ref>). This is shown in
Fig. <ref> where the invariant cross section is presented as a function of $p_T$ and several values of $x_F'$ and $\log(s)$.

Invariant $\pi^-$ cross section $f/\sigma_{\textrm{inel}}$ as a function of $p_T$ for several $x_F'$ and $\log(s)$ values

In fact there is only a very small region at $\sqrt{s} \sim$ 11 GeV where
an exponential shape is realized at $x_F' \sim$ 0, in this area even up to large $p_T$ values of $\sim$ 4 GeV/c [33]. Always of course exception made for small $p_T$ where the invariant cross section has to approach $p_T$ = 0 GeV/c with
tangent zero. Towards higher interaction energies, where only data at central rapidity are available, the again very complex $p_T$ and $\log(s)$
dependence has already been shown in Fig. <ref>.

Notwithstanding this complex experimental situation, rather general
attempts at understanding the transverse phenomenology by way of statistical or thermal "models" have been and still are enjoying
widespread interest. These attempts are characterized by only one
or a few parameters besides the particle mass. They will be confronted
with the results of resonance decay in the following Sections.

§.§ General considerations concerning the choice of co-ordinates

The in-depth analysis of transverse momentum phenomena needs
first of all some clarification, in particular concerning the
coordinate systems to be used. Unlike the longitudinal momentum
distributions discussed in the preceding Section in different
reference systems, the use of orthogonal coordinates is mandatory. This is evident in Fig. <ref> where the mean transverse momentum distributions are shown as functions of $x_F$ and rapidity for the
NA49 data at $\sqrt{s}$ = 17.2 GeV.

Mean transverse momentum $\langle p_T \rangle$ as a function of
a) $x_F$ and b) rapidity integrated in the limit 0 $< p_T <$ 1.3 GeV/c

The completely different functional behaviour is explained by the fact that the rapidity variable is not orthogonal in $p_L$ and $p_T$ but links both variables such that, for pions, the cross sections are closely following, for constant rapidity, fixed cms angles thus folding the $p_T$ distributions with the strongly decreasing longitudinal yields as presented in Fig. <ref>.

Kinematical relation between the orthogonal longitudinal and transverse momentum variables $p_L$ and $p_T$ and lines of constant rapidity (full lines). For different interaction energies $\log(s)$ the kinematic limits are indicated (dashed lines)

For the NA49 energy the kinematical limit is reached for the limiting $p_T$ value of this study of 1.3 GeV/c at $y$ = 2.56 at which point the pion yield has to vanish. Nothing could be more erroneous than to interpret the decreasing $\langle p_T \rangle$ as a function of $y$ (Fig. <ref>b) as a reduction of "hadronic temperature" in forward direction. On the contrary, there is a strong increase of $\langle p_T \rangle$ with increasing $p_L$ or, as shown in Fig. <ref>a), with $x_F$ – an effect known since decades as "seagull effect" [95, 96].

In connection with the determination of the mean transverse momentum the pion densities as functions of $x_F$ and rapidity are used:

\begin{align}
\frac{d^2n}{dx_Fdp_T} &= \frac{p_T}{E}f(x_F,p_T) \\
\frac{d^2n}{dydp_T} &= p_T f(y,p_T)
\end{align}

These two pion densities have a different shape as shown in Fig. <ref>.

a) Pion densities $d^2n/dx_Fdp_T$ for $x_F$ = 0 and b) $d^2n/dydp_T$ for $y$ = 0 as a function of $p_T$ at $\sqrt{s}$ = 17.2 GeV/c

This results in $\langle p_T \rangle$ values as

\begin{align}
\label{eq:mpt_xf}
\langle p_T \rangle_{x_F} &= \frac{\displaystyle \int p_T \,\frac{dn}{dx_Fdp_T} \,dp_T}
{\displaystyle \vphantom{\sum_{a}}\int \frac{dn}{dx_Fdp_T}\, dp_T} = \frac{\displaystyle \int \frac{p_T^2}{E} \,f\, dp_T}
{\displaystyle \int \frac{p_T}{E} \,f \,dp_T} \\
\label{eq:mpt_y}
\langle p_T \rangle_{y}   &= \frac{\displaystyle \vphantom{\sum^{a}}\int p_T \,\frac{dn}{dydp_T}\, dp_T}
{\displaystyle \vphantom{\sum_{a}}\int \frac{dn}{dydp_T}\, dp_T} = \frac{\displaystyle \int p_T^2 \,f \,dp_T}
{\displaystyle \int p_T\, f\, dp_T}
\end{align}

The resulting $\langle p_T \rangle$ values differ appreciably by 0.1 GeV/c for $\pi$, 0.06 GeV/c for K and 0.03 GeV/c for p($\overline{\textrm{p}}$). This difference depends on the particle mass as presented in Tab. <ref>.

$\langle p_T \rangle_{x_F}$ [GeV/c]      $\langle p_T \rangle_{y}$ [GeV/c]
$\pi^+$                       0.258        0.352

$\pi^-$                       0.256        0.351

K$^+$                         0.406        0.466

K$^-$                         0.395        0.451

K$^0_S$                       0.402        0.460

p                             0.495        0.528

$\overline{\textrm{p}}$       0.477        0.507

$\langle p_T \rangle_{x_F}$ for $x_F$ = 0 and $\langle p_T \rangle_{y}$ for $y$ = 0 at $\sqrt{s}$ = 17.2 GeV of feed-down subtracted data for different particle species

A last general remark concerns the influence of the feed-down correction on the determination of mean $p_T$. As shown in Fig. <ref> at $x_F$ = $y$ = 0 the $p_T$ distributions are different with and without feed-down correction, see Sect. <ref>.

$f/\sigma_{\textrm{inel}}$ for $\pi^-$ as a function of $p_T$ at $x_F$ = $y$ = 0 with and without feed-down subtraction

The corresponding $\langle p_T \rangle_{x_F}$ and $\langle p_T \rangle_{y}$ values are different as given in Tab. <ref>.

$\langle p_T \rangle_{x_F}$ [GeV/c]      $\langle p_T \rangle_{y}$ [GeV/c]
feed-down subtracted        0.256        0.351

feed-down included          0.249        0.343

$\langle p_T \rangle_{x_F}$ for $x_F$ = 0 and $\langle p_T \rangle_{y}$ for $y$ = 0 at $\sqrt{s}$ = 17.2 GeV of feed-down subtracted and feed-down included data

§.§ Statistical and thermodynamical models

Transverse momentum distributions, if plotted for constant $s$ and different $y_{\textrm{lab}}$ (Fig. <ref>) or different $x_F'$ (Fig. <ref>) or at central rapidity for different s (Fig. <ref>), exhibit an extraordinary variety of shapes that seems to defy any attempt at a simple algebraic description covering the full phase space. In the preceding Sections, mostly concentrated on longitudinal momentum dependencies, the transverse component was rather characterized by its rapid, quasi-exponential and quasi invariant, decrease ("longitudinal phase space") without attempts to describe the detailed shape which is, of course, not calculable in non-perturbative QCD.

On the other hand, early on in the 1950's, rather general statistical models [97, 98] based on thermodynamical analogous were developed and applied to hadronic interactions and particle production.

This approach has been extended to predictions concerning momentum distributions, principally also the transverse component, in the mid-1960's [99] and quantified in subsequent publications [100], [74].

The evolution of interaction energies through the ISR and RHIC regions up to the LHC necessitated, however, important modifications and extensions concerning the shapes of the transverse distributions. This was attempted by adding at least one further parameter to the original thermodynamic fits, or taking reference to another statistical mechanical formulation [101, 102].

To date the original claim of thermodynamical behaviour of hadronic production, in particular concerning the notion of "hadronic temperature" and, in connection with a "limiting temperature", the purported existence of a hadronic phase transition, is still widely proclaimed. In the following sections, a critical view at this problematics will be attempted for the sector of elementary hadronic interactions.

§.§.§ Hagedorn's Statistical Bootstrap Model (SBM)

The "bootstrap" yields an exponentially increasing mass spectrum of "fireballs" or "resonances" whose inverse slope defines a limiting "temperature" – hence the inherent possibility of a hadronic phase transition should this limiting temperature ever be exceeded. The momentum distributions of fireballs, resonances and their decay products are postulated to obey the characteristics of an ideal gas in equilibrium with only longitudinal motion. This leads, taking account of the different particle masses via the transverse mass $m_T$ (<ref>) to the prediction of exponential $m_T$ distributions with a unique inverse slope defined by the hadronic temperature $T$. This means that the SBM model should describe the transverse momentum distributions by a single parameter, independent of interaction energy and particle type.

A first confrontation with the experimental reality is given in Fig. <ref> where $m_T-m$ distributions for pions, kaons and protons are shown for the NA49 data at $\sqrt{s}$ = 17.2 GeV and central rapidity.

Inverse slopes of the transverse mass distributions as a function of $m_T-m$ for a) pions [17], b) kaons [52] and c) protons and anti-protons [53] at $\sqrt{s}$ = 17.2 GeV

Although these distributions have comparable shapes they are far from showing a constant inverse slope. In fact minima are reached at 160 MeV for $m_T-m$ between about 150 and 350 MeV (pions), at 150 MeV for 50 MeV (kaons), at 140 MeV for 150 MeV (protons) and at 130 MeV for 100 MeV (anti-protons). The strong increase below these minima is due to the flat approach to $p_T$ = 0 GeV/c (Fig. <ref>). All plots increase to temperatures of about 200 MeV above the minima at the $p_T$ limit of this study. This value is far above the "limiting temperature" given by Hagedorn [103] as 160$\pm$10 MeV.

The situation gets more involved when using the global interpolation and LHC results over the complete energy region as shown in Fig. <ref>.

Inverse slopes of the transverse mass distributions as a function of $m_T-m$ for $\pi^-$ from $\log(s)$ = 1.0 to 8.0

Only the region of $\sqrt{s} \sim$ 10 GeV (Serpukhov) exhibits a rather constant inverse slope which corresponds to the exponential $p_T$ distribution in this energy range (Figs. <ref> and <ref>). Further information is given in Fig. <ref> where $T$ is presented for constant values of $m_T-m$ as a function of $\log(s)$.

Inverse slopes $T$ for constant $m_T-m$ as a function of $\log(s)$

Here again only at $\log(s)$ = 2 there is a reasonably small range of inverse slopes around 160 MeV. Below and above this energy there is a wide spread of temperatures from 100 to 300 MeV. A similar result has been obtained for kaons [52] as shown in Fig. <ref>.

Inverse slopes $T$ for constant $m_T-m$ as a function of $\sqrt{s}$ for a) K$^+$ and b) K$^-$

The increase of inverse slopes at low $m_T-m$ (Fig. <ref>) may be attributed to the fact that a purely exponential approach becomes unphysical at $m_T-m$ = $p_T$ = 0 since the invariant cross sections have to flatten out such that the axis at $p_T$ = 0 GeV/c is crossed horizontally. The increase beyond $m_T-m \sim$ 0.4 GeV/c$^2$ poses a more principle problem as it continues to $p_T$ values much beyond the limit imposed in this paper. Originally this was attributed to occasional thermodynamic turbulence, to local temperature fluctuations or even to the "escape" of hadrons from a plasma state. An explanation seemed to be offered by QCD in the perturbative sector which allowed the description of the "high $p_T$" results at ISR and the SPS $\overline{\textrm{p}}$p collider as basically a power-law behaviour. Hence a second parameter beyond the limiting temperature was introduced [100] in this sense. A problem is however posed by the fact that perturbative QCD is certainly not applicable in a $p_T$ range of 0.4 to 1.3 GeV/c. In addition it will be shown in Sect. <ref> that particle yields are saturated at least up to $p_T$ = 3 GeV/c by resonance decay. For a practical application of statistical thermodynamics in an adapted two-parameter form see Sect. <ref>.

There are other inconsistencies inherent to the SBM approach. The "fireballs" are supposed to have a sequence cascading from higher to lower mass states. In the thermodynamic sense each fireball should follow a Boltzmann-type thermal transverse distribution such that already in a second generation chain m$^*$ $\rightarrow$ m$^{**}$ $\rightarrow$ final state hadrons these hadrons would have a non-thermal $m_T$ distribution. This has lead to the introduction of an "effective temperature" that is far above the thermal one up to cascading masses m$^*$ of 5 GeV and more [100]. Finally, only the very last step to final hadrons is taken into account with m$^*$ having no transverse momentum. Furthermore, even in this approximation, two-body decays pose a grave problem as they would yield specific lines in $p_T$. In order to describe these cases which are common in hadron spectroscopy, $m^*$ has to have at least a thermal $m_T$ distribution. Even then the resulting final state hadron distributions are non-thermal. It is supposed that two-body decays only happen at the very end of the cascading chain and would not have incidence on the overall $m_T$ distributions. However for most of the heavier hadronic decays the multibody final states regroup into effective two-body decays. This is apparent both for baryonic and mesonic resonances, for instance

\begin{equation}
\begin{split}
N^*(2050) \rightarrow& N+3\pi       \\
&\rightarrow N+\eta       \\
&\rightarrow N+\omega     \\
\rightarrow& N+2\pi       \\
&\rightarrow N+\rho       \\
&\rightarrow \Delta+\pi    \\
&\rightarrow N+2\pi       \\
\textrm{or}         &                         \\
f_2(2300) 	\rightarrow& \phi+\phi     \\
&\rightarrow K+\overline{K} + K+\overline{K}
\end{split}
\end{equation}

which contain effective two-body decays with large branching fractions.

It is also claimed [100] that already for three-body decays the secondary distributions would be thermal which is not the case, see Sects. <ref> and <ref>.

§.§.§ Application of Tsallis statistics

The use of the Tsallis form of statistical thermodynamics [101] to hadronic production [102] is not uncontested [104].

Nevertheless it offers a two-parameter fit similar to extension of the purely thermal fits given in [100] of the form

\begin{equation}
\label{eq:tsallis}
f = C m_T\left( 1 + (q-1)\frac{m_T}{T} \right)^{-\frac{1}{q-1}}
\end{equation}

Fits of this type are widely used in the energy region from RHIC to LHC. In the lower energy range, $\log(s)$ = 1.0–3.5, using the general interpolation of $\pi^-$ results, the resulting fits with the parameters $T$ and $q$ are given in Fig. <ref>.

Tsallis fit parameters $T$ and $q$ (<ref>) for $\log(s)$ from 1 to 3.5 as a function of rapidity

A wide spread of the fit parameters is obtained which are not described by a common behaviour. In particular at the lower interaction energies the approach to the kinematic limit has to be taken into account leading to a progressive shrinking of the $p_T$ range (see also Fig. <ref>). A line corresponding to a $p_T$ range of 0.2 GeV/c is indicated in Fig. <ref>.

Further indications of the quality of the fits are given in Fig. <ref> for $\log(s)$ = 2.5 (NA49 energy range). The cross sections $d^2n/p_Tdp_Tdy$ are given for rapitities from 0 to 3.6 together with the residual distributions as a function of both $m_T*\cosh(y)$ and $p_T$ indicating the strong decrease of the $p_T$ range at forward rapidities.

Tsallis fits of $\pi^-$ cross sections for rapidities between 0 and 3.6 at $\log{s})$ = 2.5. The abscissa shows both the quantity $m_T*\cosh(y)$ and $p_T$. For each rapidity the fit parameters $T$ and $q$ as well as the residual distributions of the fits are shown

§ $P_T$ INTEGRATED DISTRIBUTIONS AND TOTAL $\PI^-$ YIELDS

Following the discussion of the double differential results in the preceding chapters integrated distributions up to the determination of the resulting total $\pi^-$ yields will be presented here. There are several constraints and boundary conditions to be considered:

* The global interpolation is limited in $p_T$ to an upper limit of 1.3 GeV/c. The influence of this limit on the integration will be quantified by extending the range to 1.9 GeV/c using an exponential extrapolation.

* There is an upper limit also in $x_F$ which depends on $\sqrt{s}$ and is imposed by the limit of published cross sections at about 1 $\mu$b. Also here an extrapolation to higher $x_F$ will be used to quantify the influence on the integral.

* Results both with and without feed-down subtraction will be presented in order to quantify this important contribution.

* Due to the limited phase space coverage of the high-energy colliders the results only reach up to ISR energy at $\log(s)$ = 3.6.

§.§ $dn/dx_F$ distributions

$dn/dx_F$ distributions are shown in Fig. <ref> as a function of $x_F$ for several values of $\log(s)$ from the lowest available value to the highest ISR energy for feed-down subtracted data.

a) $dn/dx_F$ distributions as function of $x_F$ for several $\log(s)$ values from 1 to 3.6 using the $p_T$ integration up to 1.3 GeV/c, b) percent increase of $dn/dx_F$ extending the $p_T$ range from 1.3 to 1.9 GeV/c, c) percent increase of $dn/dx_F$ by adding the feed-down contribution, d) variation of $dn/dx_F$ in the high-$x_F$ region relative to $\log(s)$ = 2.5 normalizing all distributions to unity at $x_F$ = 0.4

The distributions span five orders of magnitude in $dn/dx_F$ and feature a strong $\log(s)$ dependence for $x_F \lesssim$ 0.2 followed by a region of shape similarity up to the highest $x_F$ at around 0.75. The $p_T$ integration is performed up to the limit of the global interpolation at $p_T$ = 1.3 GeV/c. The increase of $dn/dx_F$ when extending the upper $p_T$ limit to 1.9 GeV/c is shown in panel b). The increase is limited to 2% in the upper range of $\log(s)$. Panel c) presents the % increase of $dn/dx_F$ when adding the feed-down contribution. The percentage increases steadily with $\log(s)$ up to a maximum
of about 10% at low $x_F$ followed by a strong decrease towards high $x_F$. Panel d) addresses the shape similarity of the $dn/dx_F$ distributions in the high-$x_F$ region. The distributions are normalized to unity at $x_F$ = 0.4 and their ratio to a reference at $\log(s)$ = 2.5 is plotted.

This last panel allows a rather precise check of longitudinal scaling. This concept has been treated in Sect. <ref> on the double differential level of cross sections. With respect to the reference line at $\log(s)$ = 2.5 or $\sqrt{s}$ = 17.8 GeV a steady evolution of the $\pi^-$ density is evident. The line is approached from below with increasing energy and a short plateau between 0.5 and 0.7 in $x_F$ and between 2.4 and 2.8 in $\log(s)$ is followed by a further increase at ISR energy. This phenomenology does not correspond to a manifestation of the down-quark structure function as postulated  by some models [122, 123]. In this context it is regrettable that the far forward region is not attainable at the high energy colliders above the ISR energy range.

§.§ $y$ and $y_{\textrm{lab}}$ distributions

$p_T$ integrated rapidity distributions are shown in Figs. <ref> and <ref> for a subset of interaction energies ranging from $\sqrt{s}$ = 3 GeV to the highest ISR energy at 63 GeV. The linear plot of Fig. <ref> brings out the well-known features of such distributions: a strong increase of the central $\pi^-$ density with interaction energy developing into a "rapidity plateau" followed by a steep decrease into the forward direction.

$dn/dy$ distributions as functions of $y$ for several $\log(s)$ values from 1 to 3.6 using the $p_T$ integration up to 1.3 GeV/c a) with a linear scale on the ordinate  b) with a logarithmic scale for a dense coverage in $\log(s)$ from 1 to 3.6 in steps of 0.1

The presence of an increasing rapidity plateau fostered ideas of a "central" production mechanism as opposed to a "fragmentation region" in forward direction with little if any energy dependence.

The study of the central rapidity region using double differential cross sections in Sect. <ref> shows however a more complex situation. In fact there are two components contributing to the increasing particle yield: a first part at low $p_T$ typical of low-$Q$ resonance decay as it is also seen in the feed-down from weak decays, Sect. <ref>, followed by a steady increase towards higher $p_T$ again due to the decay of resonances with increasing mass hence high $Q$ decay.

The steady broadening of the $y$ distributions with $\log(s)$ is a trivial consequence of the logarithmic increase of the extent of the longitudinal phase space with energy (<ref>). The decrease of the cross section towards higher rapidity is closely conformal as shown in Fig. <ref>b) for a dense sequence of $\log(s)$ values. Hence the ideas of "energy scaling" in this region.

The scaling hypothesis is further detailed by replacing the rapidity variable by $y_{\textrm{lab}} = y_{\textrm{beam}} - y$ (<ref>) thus taking account of the increasing phase space and eliminating one trivial component. This is shown in Fig. <ref> which in addition presents two normalizations of the yield: In panel a) the density per inelastic event,

\begin{equation}
\frac{dn}{dy} = \frac{\pi}{\sigma_{\textrm{inel}}} \int f dp_T^2
\end{equation}

and panel b) the integrated invariant cross section

\begin{equation}
\frac{d\sigma}{dy} = \pi \int f dp_T^2
\end{equation}

$y_{\textrm{lab}}$ distributions of a) the rapidity density per inelastic event and b) the integrated invariant cross section as a function of $y_{\textrm{lab}}$ for seven values of $\log(s)$

As already evident in the double differential case, Sect. <ref>, the integrated invariant cross section shows less $s$-dependence than the rapidity density. This is another evidence for the necessity to properly address the variation of the inelastic cross section as a function of interaction energy and, connected to this, the question of the impact parameter dependence of a given inclusive phenomenology.

§.§ Total $\pi^-$ yields

The $p_T$ integrated distributions discussed above may now be integrated over $x_F$ or $y$ in order to determine the total $\pi^-$ yield as a function of $\log(s)$. The resulting numerical values are given in Table <ref> both with and without feed-down contribution. The extrapolation from the $p_T$ limit at 1.3 GeV/c to 1.9 GeV/c as well as the extrapolation from $x_F$ = 0.75 to $x_F$ = 0.95 contribute each less than 0.4% to the yield.

$\log(s)$     1.0       1.1       1.2       1.3       1.4       1.5       1.6       1.7       1.8       1.9
with FD		   0.2649    0.3574    0.4736    0.6077    0.7473    0.8844    1.0214    1.1659    1.3229    1.4920

without FD	   0.2581    0.3446    0.4527    0.5767    0.7049    0.8305    0.9557    1.0874    1.2306    1.3847
$\log(s)$     2.0       2.1       2.2       2.3       2.4       2.5       2.6       2.7       2.8       2.9
with FD		   1.6710    1.8572    2.0479    2.2422    2.4395    2.6417    2.8496    3.0685    3.2986    3.5403

without FD	   1.5474    1.7161    1.8883    2.0631    2.2406    2.4238    2.6125    2.8113    3.0209    3.2419
$\log(s)$       3.0       3.1       3.2       3.3       3.4       3.5       3.6
with FD		   3.7895    4.0504    4.3200    4.5963    4.8778    5.1648    5.4547

without FD	   3.4700    3.7094    3.9578    4.2137    4.4757    4.7440    5.0158

Total $\pi^-$ yield as a function of $\log(s)$ with and without feed-down contribution

These values are shown in Fig. <ref> as a function of $\log(s)$ together with the contribution of the feed-down component in percent.

Total $\pi^-$ yields as a function of $\log(s)$ with and without feed-down contribution (left scale) together with the feed-down component in percent (right scale)

§ RESONANCES AND THEIR DECAY

It has been known since decades that final state hadrons are decay products of resonances [75, 76, 105] to a fraction which ranged, at the time, from 80 to 90%. To date one may assume that all of these hadrons stem from resonances, in other words there is no "direct" production. Resonances form a cloud of states above the ground state of stable hadrons which is of extraordinary diversity and extent concerning their masses, quantum numbers and decay branching fractions. Their study is thus indispensable for any progress of understanding that goes beyond the realm of purely inclusive single particle physics.

In order to determine the contribution of a given resonance to the final state inclusive particle cross section several quantities have to be known:

* Central mass and width

The resonance width may change from zero for weak decays (Sect. <ref>) to several hundred MeV for high-mass, strongly decaying states.

* Mass distribution

Resonances are in general characterized by asymmetric Breit-Wigner type mass distributions which may extend up to masses far above the nominal value. Their important influence on tail effects in the inclusive cross sections will be evoked below.

* Decay branching fractions

For non-exotic decays as they are treated in this paper, branching fractions span a range from a few to almost 100%. They are generally well measured for mesonic and heavy flavour resonances as they are readily accessible in e+e colliders. For baryon resonances the situation is less favourable. Here in general only a few channels have been isolated with rather large error margins.

* Resonance production cross sections over full phase space

Here the production yields have in principle to be known, on a double-differential level, over the full phase space. Given the problems of coverage for the final state hadrons evoked in this paper, it is not surprising that this condition is hardly fulfilled for most resonances. Here single-differential distributions integrated over $p_T$ or $x_F$ will have to do in most cases.

In the following sections, a particular resonance, namely the $\Delta^{++}$(1232), will be used to clarify some details concerning the problems mentioned above. This resonance has been studied by about 30 experiments over a wide range of interaction energies and for different projectiles on proton targets. The $\Delta^{++}$ has an almost 100% branching fraction into the two-body channel p + $\pi^+$. The asymmetry in the decay particle masses allows, in addition, to elaborate the mass effect on the final state particle distributions.

Some references for p+p interactions are given in [106, 107, 108, 109, 110, 111, 112, 113, 114, 115].

These earlier publications concern mostly bubble chamber work and are concentrated in the 1970's, see the discussion in Sect. <ref>, and Fig. <ref>. Refs.[112, 113] come from spectrometer experiments at the ISR and SPS.

Experimental information on the $\Delta$(1232) quadruplet at $\sqrt{s}$ = 17.2 GeV is available from two Theses [114, 115] in the context of the NA49 experiment.

§.§ $\Delta^{++}$ decay: General considerations

The phase space distribution of hadrons from resonance decay is characterized by a two-step process.

The first step is given by the disintegration process in the resonance cms which is defined by the resonance and the decay particle masses. For a two body decay, the secondary particles have equal momentum independent of their respective masses. This disintegration is modulated by a spin dependent decay angular distribution which is in general not uniform.

In a second step a Lorentz transformation from the resonance cms to the overall cms as defined by the colliding initial state hadrons has to be performed. As the resonance has both longitudinal and transverse momentum in this system the transformation is three-dimensional.

There are two complications involved with this straight-forward decay process. Firstly, the resonance has in general not a fixed mass but features a complex asymmetric mass distribution with a long tail towards high masses. Secondly, a three dimensional boost is not easily amenable to simple algebraic formulations. Here it is specifically the resonance transverse momentum which creates uncomfortable complications. It is therefore rather customary to neglect the one or the other complication (or both) in hadron production models. Although it has generally to be admitted that a hadronic collision results in a spectrum of massive objects (often named "fireballs") their decay into the final state hadrons is rarely specified in detail. As an example  the "NOVA" model [116] specifies a rather narrow Gaussian-type "fireball" mass spectrum centred at relatively low masses, but it does not allow these objects to have transverse momentum. On the other hand Hagedorn's "Thermal" model [99] allows for a fireball spectrum but again neglects their transverse momentum. In addition final state particle production is assumed to be described by Boltzmann radiation from an equilibrium thermodynamic state where mass dependence is introduced at this level only.

In view of this situation $\Delta^{++}$ decay will be treated in the following Sections in some detail by starting with the fixed nominal mass 1.232 GeV/c$^2$ and an $x_F$ distribution without allowing for transverse momentum. The transverse momentum and the Breit-Wigner mass distributions will then be successively introduced.

§.§ $\Delta^{++}$ resonance: $x_F$ and $p_T^2$ distributions

From the rich sample of experimental results mentioned above the phase space distribution of the Delta resonance may be extracted, in the present argumentation at the NA49 energy of 17.2 GeV. These distributions are shown in Fig. <ref>.

Single differential $\Delta^{++}$ distributions at $\sqrt{s}$ = 17.2 GeV: a) $dn/dx_F$ as a function of $x_F$, b) $dn/dp_T^2$ as a function of $p_T^2$

The $dn/dx_F$ distribution integrates to 0.22 $\Delta^{++}$ per inelastic event at $\sqrt{s}$ = 17.2 GeV. This contributes, given the 100% branching fraction, about 7% to the total $\pi^+$ yield and 20% to the total p yield at this energy [17, 53].

§.§ $\Delta^{++}$ decay: Mass as a delta function at 1.232 GeV and no $p_T$ allowed

The Lorentz transformation between the resonance cms and the p+p cms is defined by the momentum vector of the resonance and decay particle masses via the energy and momentum values in the resonance cms. For the nominal $\Delta^{++}$ mass the decay momentum is $q$ = 0.227 GeV/c independent of the particle mass $m_{\textrm{dec}}$ which enters only via the energy factor

\begin{equation}
E_{\textrm{dec}} = \sqrt{ q^2 + m_{\textrm{dec}}^2}
\label{eq:edec}
\end{equation}

Excluding the resonance transverse momentum, the $p_T$ of the decay particles is limited by the decay momentum $q$ to 0.227 GeV/c, independent on particle mass. This results in the mean $p_T$ as a function of $x_F$ as shown in Fig. <ref> for the decay protons and pions.

Mean $p_T$ for decay protons and pions as a function of $x_F$

The Lorentz boost is purely longitudinal and yields, for a given resonance momentum $p_{l,\textrm{res}}$, the Lorentz factors

\begin{align}
\beta =& \frac{p_{l,\textrm{res}}}{\sqrt{ p_{l,\textrm{res}}^2 + m_{\Delta}^2 }} \label{eq:beta}\\
\textrm{and} & \nonumber \\
\gamma =& \frac{1}{\sqrt{ 1 - \beta^2 }} \label{eq:gamma}
\end{align}

and hence the decay particle longitudinal momenta between

\begin{align}
p_{l,\textrm{min}} &= \gamma( -q + \beta E_{\textrm{dec}} )  \label{eq:pmin} \\
\textrm{and} & \nonumber \\
p_{l,\textrm{max}} &= \gamma( +q + \beta E_{\textrm{dec}} )  \label{eq:pmax}
\end{align}

This leads to an explicit mass dependence in the $x_F$ dependence as shown in Fig. <ref> for the $p_T$ integrated quantity $dn/dx_F$.

$dn/dx_F$ for decay protons and pions as a function of $x_F$

The neglecting of the resonance transverse momentum and mass distribution is rather common in model calculations, albeit for practical considerations concerning the possibility of algebraic solutions [99, 116]. It leads to the wide-spread belief that low-mass secondaries are all centred at low $x_F$ and low $p_T$ [105].

§.§ $\Delta^{++}$ decay: adding resonance transverse momentum

In taking account of the measured $p_T$ distribution of the $\Delta^{++}$ the Lorentz transformation gets a transverse component. At $x_F$ = 0. this transverse boost replaces $p_L$ with $p_T$ in Eqs. <ref> to <ref> thereby creating a mass dependence also in transverse direction. This is shown in Figs. <ref> and <ref> for the mean $p_T$ and $dn/dx_F$ as a function of $x_F$ both for the decay pion and proton.

Mean $p_T$ for decay protons and pions as a function of $x_F$

$dn/dx_F$ for decay protons and pions as a function of $x_F$

It is to be noted that the increase of $\langle p_T \rangle$ for pions varies from 16% at $x_F$ = 0 to 300% at $x_F$ = 0.4 whereas for protons it is much larger at $x_F$ = 0 with 300% decreasing slightly to 250% at $x_F$ = 0.7. This originates from the mass dependence of the energy factor also in the transverse Lorentz boost (<ref>) and explains the observation that heavy particles have in general considerably higher $\langle p_T \rangle$ than lighter ones at low $x_F$.

For $dn/dx_F$ there is no noticeable variation both for protons and for pions as compared to Fig. <ref> without resonance $p_T$.

§.§ Delta decay: adding the resonance mass distribution

The Breit-Wigner (BW) mass distribution of resonances has a decisive influence both on mean $p_T$ and $dn/dx_F$ for pions whereas for protons there are only relatively small modifications. This is again due to the fact that tails in the $Q$ value of the Lorentz boosts (<ref>) are most effective for small decay mass values $m_{\textrm{dec}}$.

Limitations to the shape and mass range of the BW distribution will be discussed in the following Sect. <ref>. In the present explanatory section a BW distribution with linear mass damping up to 3 GeV for the $\Delta^{++}$ will be used.

The most important effect of the resonance mass distribution consists in the fact that the decay pions now cover the complete $x_F$ range and also reach the $\langle p_T \rangle$ values of the protons in the higher $x_F$ region. This is demonstrated in Fig. <ref> for $dn/dx_F$ and in Fig. <ref> for $\langle p_T \rangle$.

$dn/dx_F$ for decay protons and pions as a function of $x_F$

Mean $p_T$ for decay protons and pions as a function of $x_F$

The $dn/dx_F$ distribution for protons increases by 20% at $x_F$ = 0 and decreases by about the same amount at $x_F$ = 0.55 with respect to Fig. <ref>. For pions there is a decrease of about 10% up to $x_F \sim$ 0.3. Beyond this limit there is a dramatic increase of the yield such that the cross section of the decay pions reaches a level of about 20% of the total inclusive $\pi^+$ yield (also shown in Fig.<ref> up to the highest $x_F$ range.

For the mean $p_T$ a similar picture emerges. There is a strong increase of 30% already at $x_F$ = 0 compared to Fig. <ref> which continues rising until it crosses and exceeds the proton value at $x_F$ above 0.45. The relative factor between $\langle p_T \rangle$ for protons and pions,

\begin{equation}
R_{\langle p_T \rangle}(x_F) = \frac{\langle p_T \rangle_{\textrm{prot}}(x_F)}{\langle p_T \rangle_{\textrm{pion}}(x_F)}
\end{equation}

is shown in Fig. <ref>.

$R_{\langle p_T \rangle}$ as a function of $x_F$ for the different decay configurations Sect. <ref>–<ref>

An overview of the double-differential pion cross section $f(x_F,p_T)$ is presented in Fig. <ref> showing the total inclusive $\pi^+$ yield, the decay pion distribution at fixed resonance mass and its evolution allowing for a Breit-Wigner mass distribution.

Invariant pion cross section as a function of $p_T$ and $x_F$ for data [17] and decay of the baryon resonance $\Delta(1232) \rightarrow N + \pi$ with zero width and Breit-Wigner mass distribution

From the above discussion it should be clear that a detailed treatment of resonances and of their decay is absolutely mandatory in order to come to an understanding of their contribution to the measured inclusive data. In this sense resonances offer a model-independent step towards a better understanding of non-perturbative hadronic phenomena.

§ THE RESONANCE MASS SPECTRUM

The discussion of $\Delta^{++}$ decay in the preceding Sections has demonstrated the necessity of properly taking into account both the transverse momentum and the mass dependencies. Whereas the transverse momentum distribution is experimentally accessible this is in general not true for the mass distribution due to problems in extracting its tails from an important background. Here experimental constraints as well as energy-momentum conservation and cascading decays play a role. This will be discussed in the following Sections.

§.§ The unconstrained mass distribution

Performing an energy scan over the mass range of a resonance in a "formation experiment" combined with a partial wave analysis the unbiased Breit-Wigner mass distribution may be obtained:

\begin{equation}
BW(m) = \frac{m m_0\Gamma(m)}{(m_0^2 - m^2)^2+m_0^2\Gamma^2(m)}
\label{eq:bw_mass}
\end{equation}

where $m_0$ is the resonance mass and $\Gamma(m)$ is the mass and spin dependent width

\begin{equation}
\Gamma(m) = \Gamma_0 \left(\frac{Q}{Q_0}\right)^{2L+1} \frac{2Q_0^2}{Q^2+Q_0^2}
\end{equation}

$\Gamma_0$ is the central width value, $Q$ and $Q_0$ the momenta of the decay products in the resonance rest frame and $L$ the spin of the resonance. This mass distribution is shown in Fig. <ref> for the decay

\begin{equation}
\Delta(1232) \rightarrow p + \pi
\end{equation}

BW mass spectrum of the $\Delta(1232)$ resonance with full width 118 MeV

It is limited at low mass by the threshold $m_{\textrm{thr}}$ = $m_p$ + $m_\pi$ and features a long, unconstrained tail. This tail is subject, for production experiments as they are discussed here, to several constraints presented in the following sections.

§.§ Mass limitations due to the available interaction energy

As production experiments are performed at fixed interaction energy there is, by energy-momentum conservation a cut-off mass smoothly approached from below according to the $x_F$ distribution of the resonance. This is demonstrated in the original Jackson paper [117] for the interaction K$^+$ + p $\rightarrow$ K$^0$ + $\pi$ + p at 1.14 and 3 Gev/c incident K$^+$ momentum. As shown in Fig. <ref> both the resonance shape and the position of the central mass are depending on the interaction energy.

Relativistic Breit-Wigner mass distributions of the final state $\Delta^{++}$ in the reaction K$^+$ + p $\rightarrow$ K$^0$ + $\pi$ + p at $\sqrt{s}$ = 1.86 and 2.61 GeV, [117]

At higher interaction energy the central mass distribution stays largely unchanged but the high mass tails are curtailed by energy-momentum conservation as as shown in Fig. <ref> for p+p interactions in the range of beam momenta between 5 and 158 GeV/c. At the lowest beam momentum the mass distribution is conformal up to about 1.4 GeV and vanishes at about 1.7 GeV whereas at 158 GeV/c the change is smaller than 20% even at a mass of 5 GeV.

Relativistic Breit-Wigner mass distributions for final state Delta baryons in p+p interactions at beam momenta between 5 and 158 GeV/c

Correspondingly the $x_F'$ distributions (Sect <ref>) are progressively suppressed above $x_F'\sim$ 0.4 as shown in Fig. <ref>.

Suppression of $x_F'$ distributions with respect to the unconstrained case for beam momenta between 5 and 158 GeV/c

§.§ Constraints due to proton momentum cuts

In bubble chamber experiments the positive identification of protons via bubble density is only feasible for lab momenta below 1.0–1.3 GeV/c. In consequence both the high mass tails of the relativistic Breit-Wigner distribution are curtailed and the decay angular distributions become suppressed in the backward direction. This is demonstrated in Figs. <ref> for the mass distributions and Fig. <ref> for the Gottfried-Jackson angle $\Theta_{GJ}$.

Measured relativistic Breit-Wigner mass distributions of a $\Delta$(1232) resonance imposing a lab momentum cut on the decay protons at $p_{\textrm{lab}}$ = 1.0 GeV/c for beam momenta between 10 and 158 GeV/c

Suppression factors of the Gottfried-Jackson decay angular distribution in backward direction as a consequence of a proton $p_{\textrm{lab}}$ cut at 1.0 GeV/c as a function of cos($\Theta_{GJ}$) and beam momenta between 10 and 158 GeV/c

The resulting modification of the relativistic Breit-Wigner shape has an opposite energy dependence as the one imposed by energy-momentum conservation, Fig. <ref>.

The lab momentum cut also influences the acceptance of the decay angle distribution by suppressing backwards decays. This is shown in Fig. <ref> where the resulting suppression factor is plotted as a function of the Gottfried-Jackson angle.

Evidently the resulting distortions have to be corrected for in the analysis of the experimental data.

§.§ Cascading

An important aspect of resonance production and decay is cascading. Given sufficient cms energy, the lower-mass resonances are decay products of higher mass states. This means that they are not directly produced but rather cascade down from heavier objects. It is interesting to estimate the amount of cascading steps prior to the appearance of the final state particles. In the baryonic sector this may be done for final state "net" baryons invoking baryon number conservation. Given the total sum of net baryonic resonance cross sections $\sum \sigma(\textrm{res})$, the number of cascading steps $n_{\textrm{casc}}$ is given by

\begin{equation}
n_{\textrm{casc}} = \frac{\sum \sigma(res)}{2\sigma_{\textrm{inel}}}
\end{equation}

At SPS energy this may be estimated to $n_{\textrm{casc}} >$ 2 [118] where this value is to be regarded as a lower limit because the number of measured resonances gets rather limited towards high masses.

§.§.§ Cascading: Production mechanism of resonances - charge exchange versus Pomeron exchange

The interplay of charge and Pomeron exchange has been studied for baryons and pions in connection to Regge theory [90]. For baryons the $s$-dependence of the diffractive production of two-body final states of the type

\begin{align}
p + p &\rightarrow n + (p + \pi^+) \\
&  \nonumber \textrm{or} \\
p + p &\rightarrow (p + \pi^-) + (p + \pi^+)
\end{align}

approaches a slope compatible with Pomeron exchange at SPS energy [6] whereas at lower energies Reggeon exchange prevails as shown in Fig. <ref>.

Invariant cross sections of charge exchange and single and double dissociation in nucleon-nucleon interactions as a function of $s$ at a momentum transfer $t$ = 0.032 GeV$^2$. The	full line represents an interpolation of the data points. The insert gives the local slope $\beta$ in the parametrization $f \sim s^{-\beta}$ as a function of $s$

For pions the $\pi^+/\pi^-$ ratio in the target hemisphere of the isoscalar Carbon nucleus in the reaction

\begin{equation}
p + C \rightarrow \pi^+, \pi^-
\end{equation}

may be used for a similar separation of charge and Pomeron exchange [90]. As shown in Fig. <ref> the $\pi^+/\pi^-$ ratio approaches unity again in the SPS energy region whereas at lower energy it increases to values typical of charge exchange.

$\langle \pi^+/\pi^- \rangle$ as a function of $1/\sqrt{s}$ for five values of $p_{\textrm{lab}}$ between 0.2 and 1.0 GeV/c. The full lines are hand interpolations through the data points for p+C interactions in the target frame [90]

It may therefore be concluded that Pomeron exchange prevails over quantum number exchange both for diffractive and central interactions at SPS energies and above. This has the important consequence that N$^*$ production dominates there over $\Delta$ production, see also [6]. It is in sharp contrast to most "microscopic" production models where N$^*$ production tends to be neglected with respect to $\Delta$ by isospin counting arguments. In fact $\Delta$ resonances turn up in the final state as decay products of higher mass N$^*$ states with their mass distributions being inscribed in the respective N$^*$ decay distributions. Indeed three body decays of the type

\begin{equation}
N^* \rightarrow N + \pi + \pi
\end{equation}

are characterized by

\begin{equation}
N + \pi + \pi \rightarrow \Delta + \pi
\end{equation}

with large branching fractions [119].

§.§.§ Cascading: Consequences for the Breit-Wigner mass distribution

In a situation where resonances proper are decay products of higher mass states, a modification of their mass distributions with respect to the unconstrained Breit-Wigner shape, Sect. <ref>, must be expected. Due to the absence of theoretical predictivity in the soft sector of QCD it is difficult to estimate the effect of cascading in the Breit-Wigner mass tails. A look at particle production with different assumptions concerning the extent of the mass tails may help to at least establish some limits. Taking up the example of $\Delta^{++}$ decay into p and $\pi^+$ the yield of the decay products may be studied for different upper limits of the Breit-Wigner tails using a simple, linear tail suppression of the form

\begin{equation}
f_{\textrm{supp}}(m) = 1 - \frac{ m - m_{\textrm{0}} }{ m_{\textrm{up}} - m_{\textrm{0}}}
\end{equation}

where $f_{\textrm{supp}}(m)$ is the suppression factor applied to the form (<ref>) and $m_{\textrm{up}}$ an upper mass limit at which the Breit-Wigner distribution vanishes.

The resulting $dn/dx_F$ distributions for the decay pion and protons are shown in Fig. <ref> for different values of the upper mass limit $m_{\textrm{up}}$ at SPS energy ($\sqrt{s}$ = 17.2 GeV). For comparison the total inclusive yields are also indicated.

$dn/dx_F$ distributions from $\Delta^{++}$ decay as a function of $x_F$ for different damping factors $f_{\textrm{supp}}(m)$ at $\sqrt{s}$ = 17.2 GeV, a) for decay pions, b) for decay protons. The respective measured inclusive cross sections are indicated as full lines

The effect of the resonance mass tail on the decay particle yields is strong for $\pi^+$ above $x_F$ = 0.3 whereas for protons there is a smaller dependence spread over the complete $x_F$ range. This is quantified in Fig. <ref> where the yield ratios relative to the total inclusive cross section are given.

Decay particle yields relative to the total inclusive cross sections as a function of $x_F$ a) for the decay pions and b) for the decay protons at $\sqrt{s}$ = 17.2 GeV and different damping factors

It is interesting to note that the yield ratios reach more than 20% for the decay protons over most of the $x_F$ range. For $\pi^+$ the ratios increase from a few percent at low $x_F$ to about 20% at $x_F$ = 0.3 and then fan out for higher $x_F$ up to unphysical values at high $m_{\textrm{up}}$.

Further constraints on the mass cut-off may be obtained by inspecting the double-differential invariant cross sections as functions of $p_T$ for different $x_F$ values shown in Fig. <ref> for the decay pions.

Invariant decay pion cross sections as a function of $p_T$ for different $x_F$ values and mass cuts at $\sqrt{s}$ = 17.2 GeV. The measured total inclusive cross sections are shown as full lines

The corresponding yield ratios relative to the measured inclusive cross section are presented in Fig. <ref> as a function of $p_T$ for different mass cuts

Yield ratio $R_{\pi^+} = f_{\pi^+}(\Delta^{++})/f_{\textrm{incl}}(\pi^+)$ as a function of $p_T$ for $x_F$ between 0.0 and 0.7 and different mass cuts

The yield ratio exceeds the total inclusive level at $x_F >$ 0.4 for $m_{\textrm{up}}$ = 5.0 GeV/c$^2$ and reaches 1.0 for $m_{\textrm{up}}$ = 4.0 GeV/c$^2$ at $x_F >$ 0.6. At low $p_T$ the ratio stays consistently on a 10–20% level, decreases to a flat minimum at $p_T \sim$ 0.5 GeV/c and stabilizes again at 10–20% for higher $x_F$ and $m_{\textrm{up}}$ = 2.5–3 GeV.

For the decay protons the influence of the mass cut-off is less important as shown in Fig. <ref> which gives directly the yield ratio $R_{\textrm{p}} = f_p(\Delta^{++})/f_{\textrm{incl}}(p)$ as a function of $p_T$ for different values of $x_F$ and three different mass cuts.

Yield ratio $R_{\textrm{p}} = f_p(\Delta^{++})/f_{\textrm{incl}}(p)$ as a function of $p_T$ for $x_F$ from 0.1 to 0.9, a) mass delta function,  b) $m_{\textrm{up}}$ = 3.0 GeV/c$^2$ and c) $m_{\textrm{up}}$ = 5.0 GeV/c$^2$

The ratio $R_{\textrm{p}}$ is in general on the level of 20–30% of the inclusive proton yield increasing in the higher $p_T$ range from 30 to 50% at high $x_F$ and $m_{\textrm{up}}$ up to 3.0 GeV/c$^2$. At $m_{\textrm{up}}$ = 5.0 GeV/c$^2$ this increase exceeds the limit of 100% in the medium $x_F$ range. There is a general decrease of $R_{\textrm{p}}$ with values below 10% towards $x_F$ = 0.9 which is due to the transition to single proton production for highly peripheral collisions.

§.§.§ Mean $p_T$ and inverse $m_T$ slopes

The influence of the resonance mass distribution pervades, in addition to the examples discussed above, all secondary particle distributions. As two further examples the mean $p_T$ and the inverse slopes of the transverse mass distributions ("hadronic temperature") are presented below for $\Delta^{++}$ decay.

The mean transverse momentum of the decay pions and protons is shown in Fig. <ref> as a function of $x_F$ for different mass cuts $m_{\textrm{up}}$, a) for the decay pions and b) for the decay protons.

Mean transverse momentum as a function of $x_F$, a) for the decay pions and b) for the decay protons from the $\Delta^{++}(1232)$ resonance

While there is a strong dependence of $\langle p_T \rangle$ on $m_{\textrm{up}}$ for pions in particular for $x_F >$ 0.3, the mean $p_T$ of protons is less affected over the full $x_F$ range. As already discussed in Sect. <ref> (Fig. <ref>) the sizeable difference of the $\langle p_T \rangle$ values for pions and protons at small $x_F$ vanishes with increasing $x_F$ for $m_{\textrm{up}}$ between 2.0 and 3.0 GeV/c$^2$, while again the extreme assumptions of a delta function in mass or of $m_{\textrm{up}} >$ 3 GeV/c$^2$ allow for the limitation of the range of possible resonance mass distributions.

A similar conclusion may be drawn from the inverse $m_T$ slopes shown in Fig. <ref>.

Inverse slopes $T$ of the $m_T$ distributions a) for decay pions and b) for decay protons from the $\Delta^{++}$ resonance as a function of $x_F$ for different values of $m_{\textrm{up}}$

While the $T$ values of the decay protons tend to level out for $m_T - m >$ 0.3 GeV/c$^2$ after a strong increase for all $m_{\textrm{up}}$, the inverse slopes of the decay pions show, in contrast, a  strong $m_{\textrm{up}}$ dependence in the same $m_T - m$ region contradicting any "thermal" behaviour.

§.§ Summary of Sects. <ref> and <ref> concerning the decay of the $\Delta^{++}$ resonance

The in-depth discussion of several features of $\Delta^{++}$ decay as presented in Sects. <ref> and <ref> points out the importance of several resonance parameters beyond the mass, half width and branching fraction. This concerns in particular two important ingredients:

* In order to predict the inclusive distribution of the decay particles the knowledge of the resonance cross section over the full phase space is mandatory

* In strong decays a long upper mass tail (Breit-Wigner distribution) exists which depends on several parameters like energy-momentum conservation, experimental cuts and in particular cascading decay from higher mass resonances. The effect of cascading on the mass distribution is as yet not calculable as it makes part of non-perturbative QCD. It may at least be approximated by a detailed study of the phase space distribution of the decay particles with different assumptions on the extent of the upper mass tail. This way a range of possible decay parameters may be established

In the present study of $\Delta^{++}$ decay a damping of the unconstrained Breit-Wigner distribution up to mass values between 2 and 3 GeV/c$^2$ has been shown to be compatible with a wide range of secondary particle distributions.

§ INVERSE SLOPES FROM THE DECAY OF DIFFERENT RESONANCES ("TEMPERATURE")

Following the preceding discussion of the influence of the resonance parameters on a number of inclusive quantities of the final state hadrons, the inverse slope parameters of decay pions from a set of different baryonic and mesonic resonances will be discussed. In addition to the $\Delta$ baryon decay (Sect. <ref> and <ref>) the inverse slopes of $\pi^-$ from the weak decays of K$^0$, $\Lambda$ and $\Sigma^-$ (Sect. <ref>) and from $\rho^0$(770) and f$_2^0$(1270) are presented in Fig. <ref>. Here the weakly decaying strange hadrons have zero width whereas for the strong decays a linear damping of the respective Breit-Wigner distributions up to 3 GeV/c$^2$ mass is introduced. The inverse slopes $T$ are given at 158 GeV/c beam momentum and central rapidity.

Inverse slopes $T$ (GeV) as a function of transverse mass at central rapidity and $p_{\textrm{beam}}$ =158 GeV/c for decay pions from different baryonic and mesonic resonances. The measured inclusive value for pions is indicated by a thick dashed line

Several conclusions may be drawn from Fig. <ref>:

* Each resonance decay creates its proper inverse slope distribution.

* These distributions are far from flat at a given "temperature" as it would be imposed by a Boltzman-type emission from an equilibrium thermal ensemble.

* The corresponding inverse slope values are imposed by a set of resonance parameters like the $Q$ value, the central width and the mass extent of the Breit-Wigner distribution as discussed in Sect. <ref> and <ref> above.

* The actual $T$ values show a very large spread from a few tens of MeV up to several hundred MeV.

* The inclusively measured $T$ distributions of a particle must be seen as an incidental result of the overlap of many distinct individual contributions.

Given the fact that most if not all final state hadrons stem from resonance decay the very notion of a thermodynamic origin is not tenable.

§ MEAN TRANSVERSE MOMENTUM

§.§ $\langle p_T \rangle$ for pions, kaons and baryons from NA49

The NA49 experiment has provided precision data of mean transverse momenta for $\pi^+$, $\pi^-$, K$^+$, K$^-$, protons and anti-protons at $\sqrt{s}$ = 17.2 GeV [17, 52, 53]. These data are presented in Fig. <ref> as a function of $x_F$.

Mean $p_T$ as a function of $x_F$ a) $\pi^+$, K$^+$, protons and b) $\pi^-$, K$^-$, anti-protons

Several features are noteworthy:

* Opposite charges give very similar $\langle p_T \rangle$ to within about 20 MeV/c.

*  At low $x_F$ there are substantial differences in $\langle p_T \rangle$ between particle species of about 150 MeV/c between pions and kaons and about 250 MeV/c between pions and baryons.

*  These differences tend to vanish with increasing $x_F$ such that at $x_F \sim$ 0.4 all particle types yield $\langle p_T \rangle$ values which are equal to within about 50 MeV/c.

*  There are two exceptions to this effect: $\langle p_T \rangle$ for anti-protons continues to increase with $x_F$ whereas for protons there is a pronounced decrease for $x_F >$ 0.4.

In the framework of resonance decays treated in Sects. <ref> and <ref> above, the effects <ref> and <ref> are explained by the resonance transverse momentum and the resonance mass distribution, respectively, as shown in Figs. <ref> and <ref>. In this sense the data of Fig. <ref> represent a model independent check of the importance of resonance decay for the final state inclusive data. The deviations <ref> are a consequence of the detailed production mechanism. For protons, there is the approach to diffraction at large $x_F$ with sharper $p_T$ distributions [53] whereas anti-proton production is characteristic of high-mass isospin triplets of mesonic origin [120] with a steeply declining $x_F$ distribution. In this case it is the high mass tail of the respective resonances which is exploited.

In this context it should be stressed that many salient features of inclusive final state hadron distributions follow naturally from resonance decay in contrast to statistical and thermal models which have to take reference to an initial equilibrium state and neither predict the experimental inverse slopes (Sect. <ref>) nor the dependence of mean $p_T$ on $x_F$.

§.§ $\langle p_T \rangle$ for $\pi^-$ as a result of the global interpolation for different values of $\log(s)$

The interpolation scheme presented in Sect. <ref> allows for the calculation of mean transverse momentum in a range of $\log(s)$ from 1 to 3.6 or $\sqrt{s}$ from 3 to 63 GeV over the full phase space. The result is shown in Fig. <ref> as a function of $x_F$ for $\log(s)$ = 1 to 3.6
in steps of 0.4 for feed-down corrected data.

$\langle p_T \rangle$ of $\pi^-$ as a function of $x_F$ for 7 values of $\log(s)$ from 1 to 3.6

Fig. <ref> shows a definite progression of the mean $p_T$ by factors which vary from 10% at $x_F$ = 0 to 70% at $x_F$ = 0.4 over the range of $\sqrt{s}$ from 3 to 63 GeV.

Several remarks are to be made concerning the strong dependence of $\langle p_T \rangle$ on both the longitudinal momentum and the interaction energy as is evident from Fig. <ref>:

* $\langle p_T \rangle$ increases with $x_F$ from $x_F$ = 0 to $x_F \sim$ 0.4 where the increase reaches a plateau at all energies.

* There is a strong dependence on $\sqrt{s}$ which saturates at $\log{s} >$ 3.0.

* There is as well a strong increase with $x_F$ which ranges from about 20% at $\log(s)$ = 1 to 90% at $\log(s)$ = 3.6.

This situation is visualized in Fig. <ref> which shows the mean $p_T$ as a function of $\log(s)$ for several values of $x_F$.

$\langle p_T \rangle$ of $\pi^-$ as a function of $\log(s)$ for several values of $x_F$

The saturation in $\log(s)$ is due to two underlying effects. Firstly the invariant cross sections show only a small $s$-dependence above $\log(s)$ = 2.5, see Fig. <ref> for $p_T$ = 1.3 GeV/c. Secondly the limit at $p_T <$ 1.3 GeV/c leaves the tail of the $p_T$ distributions unaccounted for. The effect of this $p_T$ cut is shown in Fig. <ref> where the percentage increase of $\langle p_T \rangle$ from an upper limit at 1.3 GeV/c to an upper limit at 1.9 GeV/c is presented.

The increase of $\langle p_T \rangle$ between an upper integration limit at 1.3 GeV/c and a limit at 1.9 GeV/c in % as a function of $\log(s)$

There is a steady increase of $\langle p_T \rangle$ between the two integration limits as a function of $\log(s)$ which reaches more than 3.5% or about 0.017 GeV/c at $\log(s)$ = 3.6.

As far as the $x_F$ dependence is concerned ("seagull effect") two underlying components may be isolated using the detailed $p_T$ dependence of the cross section ratios (<ref>) shown in Figs. <ref>–<ref> at various $x_F$ values, Fig. <ref>.

The ratio $R$ of the invariant $\pi^-$ cross section to the NA49 data (<ref>) as a function of $p_T$ for different values of $\log(s)$ for a) $x_F$ = 0, b) $x_F$ = 0.05, c) $x_F$ =0.1 and d) $x_F$ = 0.3

At central rapidity, a first component at low $p_T$ exhibits an exponential shape with inverse $p_T$ slopes of about 0.18 GeV/c which is independent of the interaction energy, Figs. <ref>–<ref>. This corresponds, if converted into invariant cross sections and plotted against $m_T$ (Fig. <ref>), to inverse $m_T$ slopes of about 0.08 GeV/c. The (non-thermal) value is similar to the results for feed-down pions from weakly decaying strange hadrons (Fig. <ref>) and indicates the presence of further resonant states with low $Q$ values and small width.

A second component corresponds to the higher $p_T$ range and is clearly visible at $p_T >$ 0.5 GeV/c. It features an exponential increase with $p_T$ whose inverse slope increases with the interaction energy as shown in Fig. <ref>b. This slope becomes negative below $\sqrt{s} \sim$ 10 GeV. At the base of this behaviour is of course energy-momentum conservation at low interaction energy together with the production of resonances with increasing mass towards higher energies.

In the central region both the exponential decrease at low $p_T$ and the exponential increase at higher $p_T$ tend to compensate such that the mean $p_T$ remains small with little dependence on $\log(s)$.

In forward direction both components vanish as shown in Fig. <ref>b) to d) for $x_F$ = 0.3 up to $\sqrt{s}$ = 63 GeV, see also Figs. <ref>, <ref> and <ref>. The behaviour towards higher $\sqrt{s}$ will probably never be accessible at present-day colliders.

Here the prevailing effect of resonance decay (Sect. <ref> and the discussion in Sect. <ref>, Figs. <ref> and <ref>) leads to the strong increase of $\langle p_T \rangle$ with $x_F$ which is mostly driven by the resonance mass distribution via the Lorentz transformation from the resonance cms to the lab system.

A further remark concerns the approach to the plateau in $\langle p_T \rangle$ up from $x_F$ = 0, Fig. <ref>. With increasing interaction energy this approach becomes successively steeper indicating that the relative longitudinal momentum $x_F$ is not the correct variable to describe this phenomenon. In fact there is a factor of 20 in longitudinal momentum for constant $x_F$ between the lowest and the highest interaction energy. The mean transverse momentum is therefore plotted in Fig. <ref> as a function of cms energy $E^*$ rather than $x_F$ in order to cover the increase of $\langle p_T \rangle$ from the lowest longitudinal momentum, $E^*$ = $m_\pi$ up to $E^*$ = 3 GeV.

$\langle p_T \rangle$ of $\pi^-$ as a function of cms energy $E^*$ for seven interaction energies $\sqrt{s}$ from 3.2 to 63 GeV

The $E^*$ dependence shows a common, $s$-independent increase from $E^*_{\textrm{min}} = m_\pi$ up to $\sqrt{s}$ = 12.6 GeV. The deviations at 3.2, 5, 7.9 and 12.6 GeV above $E^*$ =0.4, 0.8, 1.4 and 2.4 GeV are imposed by energy-momentum conservation which cuts in at about 35% of the available cms energy $E^*_{max} = \sqrt{s}/2$. The saturation at $\log(s)$ above 2.6, Fig. <ref>, transforms into a subsequent lowering of the mean $p_T$ in the $E^*$ scale such that similar $\langle p_T \rangle$ values are reached multiplying $E^*$ by the $\sqrt{s}$ ratio ($x_F$ scaling).

§ PHYSICS BEYOND THE INCLUSIVE LEVEL: CONTRIBUTION FROM KNOWN RESONANCES TO THE $\PI^-$ YIELD

In the preceding sections some aspects of resonance decay and its influence on different inclusive quantities of the final state particles have been discussed. If a complete set of hadronic resonances decaying into given final state hadrons would be known including detailed yields and phase space distributions, quantum numbers, branching fractions, decay widths and mass distributions, a new level of inclusive physics beyond the single particle inclusive state might be defined. This could ultimately yield a deeper understanding of non-perturbative QCD on a purely experimental level.

Unfortunately the present state of knowledge concerning resonance production and decay does not allow such a complete study. There are however quite a number of experimental results available which may be used to obtain at least a preliminary picture of the salient features of the resulting inclusive quantities.

The following argumentation is based on a set of eight mesonic and five hadronic resonances which has been discussed some years ago [121] concerning final state $\pi^-$ production in connection with the NA49 experiment at $\sqrt{s}$ = 17.2 GeV. Although the corresponding results are necessarily incomplete some important conclusions may nevertheless be drawn.

§.§ The resonance sample

The following resonances decaying into $\pi^-$ have been included in the present study:

mesonic                   baryonic
1-1 3-3
$\eta$(548)               $\Delta^0$(1232)

$\rho^0$(770)             $\Delta^-$(1232)

$\rho^-$(770)             N$^*$(1440)

$\omega$(782)             N$^*$(1520)

f$_2$(1270)               N$^*$(1680)

In order to avoid double counting, only two-body decays into final state hadrons have been included (three-body decays for $\eta$ and $\omega$). This excludes cascading decays like

\begin{align}
\rho_3 &\rightarrow \omega+\pi \label{eq:rho3}\\
N^* &\rightarrow \Delta+\pi  \label{eq:nstar}
\end{align}

The resulting $\pi^-$ yields present therefore a lower limit, where cascading is expected to contribute preferentially to the low $p_T$ and $x_F$ regions.

§.§ Resonance yields as functions of $x_F$ and $p_T$

As there are practically no double-differential resonance cross sections available the measured $p_T$ integrated yields $dn/dx_F$ and $dn/dp_T^2$ are used for the generation of complete resonance spectra. $dn/dx_F$ distributions are presented in Fig. <ref> as a function of $x_F$ for the mesons and baryons defined above, inter- or extrapolated to the NA49 energy, $\sqrt{s}$ = 17.2 GeV.

$dn/dx_F$ distributions as a function of $x_F$, a) for a set of 8 mesonic resonances b) for five baryonic resonances

The $dn/dx_F$ distributions give, when integrated over $x_F$, the total resonance yields that serve as normalization for the generated particle yields.

$dn/dp_T^2$ distributions are given in Fig. <ref> as a function of $p_T^2$, here normalized at $p_T^2$ = 0.025 GeV/c$^2$.

Normalized $dn/dp_T^2$ distributions as a function of $p_T^2$, a) for mesonic b) for baryonic states

For the mesonic states there is an important evolution of the $p_T^2$ slopes with resonance mass. This evolution is smaller for the baryonic states as the baryonic mass differences are smaller. The distribution for protons is also given in Fig. <ref>. Here the yield difference is a function of $p_T^2$ similar to the same effect between $\Lambda$ and proton as indicated in Fig. <ref> above.

The integration of the $dn/dx_F$ distributions results in total resonance yields and, given the two-body branching fractions, in the $\pi^-$ cross sections given in Table <ref>. All yields are given per inelastic event.

Resonance             Total yield     Branching fraction       $\pi^-$ yield
$\Delta^0$(1232)        0.123                0.333               0.041

$\Delta^-$(1232)        0.0584               1.0                 0.0584

N$^*$(1440)             0.157                0.43                0.0675

N$^*$(1520)             0.0848               0.40                0.0339

N$^*$(1680)             0.0716               0.43                0.0308
$\eta$(548)             0.316                0.273               0.0862

$\rho^0$(770)           0.837                1.0                 0.837

$\rho^-$(770)           0.286                1.0                 0.286

$\omega$(782)           0.333                0.907               0.302

f$_2$(1270)             0.0709               0.62                0.0440

$\rho_3^0$(1690)        0.0566               0.236               0.0134

$\rho_3^-$(1690)        0.0397               0.236               0.0094

f$_4$(2050)             0.0159               0.113               0.00180
total baryons           0.495         $\pi^-$ from baryons       0.232

total mesons            1.955         $\pi^-$ from mesons        1.579
3ltotal $\pi^-$ from 2-body decays (3-body from $\eta$ and $\omega$)          1.811
3ltotal inclusive $\pi^-$ [17]       2.36
3lfraction from 2-body decays (3-body from $\eta$ and $\omega$)             76.7%

Resonance yields

Several consequences follow from the results of Table <ref>:

* About three quarters of all inclusive $\pi^-$ are coming from 2-body decays of measured baryonic and mesonic resonances.

* As stated above this result has to be regarded as a lower limit since higher multiplicity decays are disregarded and clearly not all possible resonant states have been included. This is especially true for strange particle decays which are not included here.

* As decay particles from two-body resonance decays are non-thermal in the sense of Hagedorn's statistical bootstrap model [100] their overwhelming contribution to the final-state inclusive yields puts grave doubts on the applicability of statistical or thermal models to hadronic interactions.

* The biggest contributions in Table <ref> come from the low-mass mesonic resonances $\eta$, $\rho$ and $\omega$. It is clear from the available branching fractions of higher mass states that these mesons are cascading down from both baryonic and mesonic resonances. Hence the necessity to restrict their contribution to their direct decays into final state hadrons. On the other hand decays of the type $\Delta \rightarrow N^*(1440) + \pi$ or $b_1(1235) \rightarrow \omega + \pi$ yield – again via two-body decays – pion contributions which have to be counted into the total inclusive cross sections.

* The resonance cross sections given in Table <ref> span three orders of magnitude and their decay pions another 4 to 5 orders of magnitude. It has therefore to be realized that the relative contribution of different resonances are strongly depending on their $Q$ values, mass and phase space distributions such that the tails of the inclusive cross sections both in $x_F$ and in $p_T$ have a different resonance heritage than their mean values.

§.§ Predictions for double-differential cross sections from resonance decay

Invariant cross sections of $\pi^-$ from the decay of the resonances given in Table <ref> are given in Figs. <ref> to <ref> as a function of $p_T$ for fixed $x_F$. They are directly comparable to the NA49 data [17].

Predicted invariant $\pi^-$ cross sections from resonance decay as a function of $p_T$ for fixed $x_F$ = 0, 0.05, 0.1 and 0.15. The predictions are given as full lines, the NA49 data as dots and the corresponding interpolations as red lines

Predicted invariant $\pi^-$ cross sections from resonance decay as a function of $p_T$ for fixed $x_F$ = 0.2, 0.25, 0.3 and 0.35. The predictions are given as full lines, the NA49 data as dots and the corresponding interpolations as red lines

Predicted invariant $\pi^-$ cross sections from resonance decay as a function of $p_T$ for fixed $x_F$ = 0.45, 0.55, 0.65 and 0.75. The predictions are given as full lines, the NA49 data as dots and the corresponding interpolations as red lines

Note that the plots are continued to $p_T$ values above the limit of 1.3 GeV/c imposed in this paper as the NA49 data extend partially up to 2 GeV/c. The simulation of the resonance decay results in $\pi^-$ cross sections which are rather exponential at low $x_F$ and $p_T >$ 1 GeV/c. The prediction has hence been extended to $p_T$ beyond 2 GeV/C.

Several features are apparent from these figures:

* On first sight there is a surprising reproduction of the main features of the inclusive distributions as far as the dependencies on $p_T$ and $x_F$ are concerned.

* At $x_F >$ 0.4 the data and their interpolation coincide with the prediction to a precision of a few percent for $p_T >$ 0.5 GeV/c.

* There is an offset centred around $p_T \sim$ 1.2$\pm$0.3 GeV/c which is $x_F$ dependent. Such offset would be expected if the resonance sum of Table <ref> would lack contributions from further resonances which are of course to be expected.

* There is an enhancement below $p_T \sim$ 0.8 GeV/c which is strongly $x_F$ dependent.

* There is another enhancement at $p_T \gtrsim$ 1.3 GeV/c and $x_F <$ 0.6.

In order to quantify these effects the offset shown in Fig. <ref> as a function of $x_F$ has been subtracted from the data interpolation.

Offset between data and $\pi^-$ from resonance decay at $p_T$ = 1.2$\pm$0.3 GeV/c in percent as a function of $x_F$

After this subtraction the percent difference between data and resonance prediction is presented in Figs. <ref>–<ref>.

Percent deviation between data and resonance sum as a function of $p_T$ for $x_F$ = 0, 0.05, 0.1 and 0.15

Percent deviation between data and resonance sum as a function of $p_T$ for $x_F$ = 0.2, 0.25, 0.3 and 0.35

Percent deviation between data and resonance sum as a function of $p_T$ for $x_F$ = 0.45, 0.55, 0.65 and 0.75

It is also useful to show the $x_F$ dependence for fixed $p_T$, Figs. <ref> and <ref>.

Percent difference $\Delta$ between data and resonance prediction as a function of $x_F$ for $p_T$ = 0.1, 0.5, 0.7, 0.9 and 1.1 GeV/c

Percent difference $\Delta$ between data and resonance prediction as a function of $x_F$ for $p_T$ = 1.3, 1.5, 1.7. 1.9 and 2.0 GeV/c

Some remarks are in place with respect to Figs. <ref>–<ref>.

* Both in $p_T$ and in $x_F$ a strong decrease of the excess yield is apparent up to about 1 GeV/c in $p_T$ and 0.2 in $x_F$.

* This is followed by a systematic increase towards higher $p_T$ and $x_F$.

* The increase with $p_T$ is continuing well beyond the actual $p_T$ limit at 2.0 GeV/c.

* The increase with $x_F$ shows a maximum at $x_F \sim$ 0.3 – 0.4 and decreases strongly towards higher $x_F$.

The low $p_T$/low $x_F$ enhancements are reminiscent of the decay of low-$Q$ resonances without Breit-Wigner tails, see Fig. <ref>. In this context the resonance decays into final-state hadrons

\begin{equation}
R \rightarrow h_1 + h_2 \qquad \textrm{($R$ - resonance, $h_1$ and $h_2$ - final state hadrons)}
\label{eq:res1}
\end{equation}

as discussed above have to be confronted with the cascading decay

\begin{equation}
R \rightarrow R' + h \qquad \textrm{($R$ and $R'$ - resonances, $h$ - final state hadron)}	\label{eq:res2}
\end{equation}

Some examples for (<ref>) have been evoked in Sect. <ref>, (<ref>) and (<ref>). Whereas in (<ref>) the Breit-Wigner tail of the decaying resonance comes fully into play, Sect. <ref>, the cascading decays like (<ref>) have been shown to damp the extent of the mass tail, Sect. <ref>. By energy-momentum conservation, the second resonance in (<ref>) has to be strongly correlated in its mass distribution with the first one: indeed the mass of the second resonance cannot be too different from the first one. This would explain the mass cut-offs necessitated by unitarity arguments shown in Sect. <ref>. The final state hadron will therefore be emitted at low $Q$ with a correlation that keeps the masses of both resonances close together with a result corresponding to Fig. <ref>. It would of course be highly desirable to look with high-statistics at cascading reactions like (<ref>) directly assessing mass correlations between $R$ and $R'$

The increase towards high $p_T$ and medium $x_F$ is, on the other hand, rather suggestive of higher mass resonances not contained in the list of Table <ref>, see for instance Fig. <ref>.

§.§ Integrated quantities: mean $p_T$

The first moments of the $p_T$ distributions of $\pi^-$ from resonance decay, Table <ref>, are shown in Figs. <ref> for mesonic and Fig. <ref> for baryonic initial states.

Mean $p_T$ for $\pi^-$ from the decay of mesonic resonances as a function of $x_F$

Mean $p_T$ for $\pi^-$ from the decay of baryonic resonances as a function of $x_F$

The $\langle p_T \rangle$ values show important differences for the individual parent resonances depending on resonance mass and two- or three-body decays. There is no trace of a common behaviour neither for the shape of the $\langle p_T \rangle$ distributions nor for their absolute values: this corresponds of course to the problems encountered with the Statistical Bootstrap Model, Sect. <ref>.

The superposition of $\langle p_T \rangle$ as a function of $x_F$ for the complete resonance sample is presented in Fig. <ref> and compared to the $\pi^-$ data from NA49 [17]. The rather precise reproduction of the main features of the experimental result proves once more the necessity to go beyond the most general inclusive level in order to develop an understanding of the underlying physics.

Mean $p_T$ for $\pi^-$ from the complete resonance sample of Table <ref> as a function of $x_F$ in comparison to the measured inclusive data from NA49

§.§ Integrated quantities: $dn/dx_F$ distributions

The integration over $p_T$ at fixed $x_F$ yields the $dn/dx_F$ distributions of the decay pions shown in Fig. <ref> for the individual resonances used in Table <ref>:

\begin{equation}
dn/dx_F = \pi/\sigma_\textrm{inel} \cdot \sqrt{s}/2 \cdot \int{f/E \cdot dp_T^2
\label{eq:dndxf}
\end{equation}

$dn/dx_F$ distributions for $\pi^-$ from the decay of the resonances used in Table <ref> as a function of $x_F$

The $dn/dx_F$ distributions show the expected steepening when compared to the parent particles, Fig. <ref>. And again, as for $\langle p_T \rangle$, there are marked differences between the different mesonic and baryonic decays. The sum over all individual contributions is shown in Fig. <ref>, and compared to the NA49 data.

$dn/dx_F$ of $\pi^-$ as a function of $x_F$ for the inclusive NA49 data and for the decay of the resonances contained in Table <ref>. A fit to the data suggesting the down-quark structure function is superimposed

Several remarks should be made in this context:

* The distribution from resonance decay tracks the $\pi^-$ data rather precisely for $x_F >$ 0.4 with a constant difference of about 25%.

* Towards lower $x_F$ this difference increases to about 60% at $x_F$ = 0. This is to be expected since the difference between data and resonance decay increases strongly for $x_F <$ 0.2, Fig. <ref>.

* Both distributions reproduce the shape of the down-quark structure function where it is thought to apply at $x_F >$ 0.35 as given by the function fitted to the data in Fig. <ref>.

* The $\pi^+$ data of NA49 [17] are presented in Fig. <ref> including a fit corresponding to the up-quark structure function. Again a rather precise correspondence is visible beyond $x_F$ = 0.35.

* The interpretation of this shape similarity invoking parton dynamics is rather daring when realizing that the final state pions come – after several cascading steps – from parent resonances where each parent creates a  proper pionic $dn/dx_F$ distribution of different shape, Fig. <ref>.

$dn/dx_F$ of $\pi^+$ as a function of $x_F$ for the inclusive NA49 data. A fit to the data suggesting the up-quark structure function is superimposed

These examples illustrate the necessity to follow inclusive phenomena beyond the lowest level of simplification thus opening up a new way towards the understanding of the underlying physics processes. Instead, the shape similarities have given rise to "recombination" models [122, 123] and to "quark counting rules" which were widely proposed and used in the 1970´s.

§.§ The approach to "high" $p_T$

Another example for the use of inclusive data for the introduction of parton dynamics in order to explain an observed phenomenology is the so-called high-$p_T$ sector. About 50 years ago it was observed at the CERN ISR collider that the $p_T$ distributions at central rapidity developed non-exponential tails above about 1 GeV/c which were in contradiction to the scaling hypotheses of the time. This has been discussed in Sect. <ref>, see for instance Fig. <ref>.

At the same time results from deep-inelastic lepton-proton scattering were used to predict the appearance of just this phenomenon at high transverse momentum by hard parton-parton scattering inside the colliding hadrons [86]. Partonic scattering amplitudes were established using the strong coupling constant to first order and the partonic structure functions available from lepton-proton experiment.

The term "high $p_T$" is rather ill-defined. As the early data were very limited in their $p_T$ range the tendency was to push the application of parton-parton scattering down to the very limit of momentum transfer compatible with perturbative QCD. This led to some problems with the reproduction of the observed yields at $p_T \lesssim$ 2 GeV/c which were partially solved by the introduction of a parton transverse momentum $k_T$ [88, 124]. Transverse momentum is however not a priory contained in the definition of partonic structure functions.

In a later stage higher-order QCD graphs including gluon radiation were introduced thus opening up a source of transverse momentum. This approach was used to compare to the data even down to $p_T \sim$ 1.5 GeV/c.

In connection with the present discussion it is however absolutely mandatory to regard resonance decay as an eventual source of "high $p_T$" hadrons. The predicted $\pi^-$ yields from the two-body decay of the resonance sample given in Table <ref> have therefore been extended up to $p_T$ values of 3 GeV/c which should be safely contained in the regime of perturbative QCD. The corresponding cross sections (at $\sqrt{s}$ = 17.2 GeV) are visible in Figs. <ref> to <ref>. As the data from NA49 at this energy are limited to $p_T <$ 2 GeV/c, available data from the ISR and Serpukhov up to $p_T$ = 3 GeV/c have been combined in Fig. <ref>.

Invariant $\pi^-$ cross sections at $x_F$ = 0.0 as a function of $\sqrt{s}$ indicated as bands with a width of $\pm$15% at $p_T$ = 2.0, 2.5 and 3.0 GeV/c. The predicted yields from resonance decay are given as points at $\sqrt{s}$ = 17.2 GeV at the same $p_T$ values with estimated errors of 30%

It is apparent from Fig. <ref> that the predicted cross sections are well within an error margin of $\pm$15% from the measured yields for $p_T$ = 2 and 2.5 GeV/c and fall below by about 40% at $p_T$ = 3 GeV/c. Here it should be remembered that the resonance sample of Table <ref> is not complete and is lacking high-mass baryons and mesons as well as heavy flavours which would contribute substantially, notwithstanding their low cross sections and small branching fractions, see the study of K$^-$ from D mesons in [52].

In general the application of perturbative QCD for the prediction of inclusive processes has to be seen very critically. This concerns also one of the generally undoubted "successes" of QCD, the production of high mass lepton pairs via the Drell-Yan effect. It has been shown 30 years ago [125] that the production of dilepton states may be well described by semi-leptonic decays of charm-anticharm and beauty-antibeauty mesons, including characteristic inclusive quantities like mass and $s$-dependence, $x_F$ and helicity angle distributions as well as average transverse momentum.

§ CONCLUSIONS

The conclusions will follow the three main sections of this paper as specified in the introduction:

* Critical review of all available data and establishment of a general interpolation scheme including results from the high-energy proton colliders.

* Confronting the interpolation with various physics hypotheses both in longitudinal and transverse direction concerned with energy dependences and eventual manifestations of partonic effects.

* Extension of the study beyond the purely inclusive level by regarding resonance production and decay in its decisive influence on virtually all inclusive phenomena.

§.§ Data evaluation

In a first step the inclusiveness of the existing data sets with respect to the treatment of $\pi^-$ feed-down from weak decays is established by isolating two classes of experiments with and without feed-down correction. Corresponding corrections over the full phase space are worked out. In a second step a subgroup of experiments yielding internally consistent cross sections is defined as "reference data". In a third step a three-dimensional interpolation scheme in the co-ordinates interaction energy, longitudinal and transverse momentum is developed  using these data. In a forth step this interpolation is confronted with the remaining "spectrometer" data which are shown to essentially necessitate corrections connected to normalization problems. In this context data sets which fall far out of all other results in terms of systematical deviations are eliminated. A detailed statistical analysis shows that the interpolation has an unprecedented systematic precision on the 5% level from $\sqrt{s} \sim$ 3 GeV up to the highest ISR energy. The interpolation is made available both with and without feed-down correction in the coordinate triplets $y_{\textrm{lab}}$, $p_T$, $\log(s)$ and $x_F$, $p_T$ and $\log(s)$.

The situation at the high-energy colliders RHIC and LHC is less favourable as the phase-space coverage of the data is reduced to the very central or very forward directions and characterized by rather substantial discrepancies between experiments. Nevertheless an attempt is made to use these data up to $\sqrt{s}$ = 13 TeV at the LHC.

§.§ Confronting the data with different physics hypotheses

Considering the wide range in interaction energy a first rather fundamental question is concerned with the data re-normalization. The inelastic cross section grows by about a factor of three from the lowest to the highest $\sqrt{s}$ available. Recent conjectures seem to indicate that this increase involves mostly the peripheral regions of the hadrons as opposed to a constant, central component. Therefore both the invariant cross section proper and the cross section divided by the inelastic cross sections are shown: due to the high precision of the interpolation, differences between the two normalizations become visible already at ISR energy. Some interesting effects up to the LHC range are exploited for $\pi^0$ production.

As far as longitudinal momentum dependencies are concerned the rather vague connection to parton dynamics by introducing the "scaling" variable $x_F$ does not yield satisfactory results when regarded over the complete $\sqrt{s}$ range and confronted with data of sufficient precision.

On the other hand the hypothesis of "limiting fragmentation" proposed at a time where very few and uncertain data at low energy were available, shows surprising predictivity even up to LHC energy. This approach has the advantage to contain built-in baryon number conservation and the definition of a cascading chain of high-mass excited states including proper treatment of the involved Lorentz-transformations.

On the level of transverse momentum there used to be a general idea of exponential damping corresponding to the early establishment of the "longitudinal phase space". A strong predictivity has been introduced via the notion of "statistical thermodynamics". This model postulates unified exponential distributions in "transverse mass" where the inverse slope defines a "temperature" that should be independent of interaction energy and particle mass and smaller than a limiting value. In connection with this a hadronic phase transition makes part of the thermodynamic treatment.

The in-depth study of the available data over the full energy range shows that this simplistic picture is not tenable. Being applicable essentially at central rapidity, the model misses as an essential part the strong increase of the mean transverse momentum as a function of $x_F$ known as "seagull effect" which is followed here with precision over the full range of interaction energies.

§.§ Resonance decay as a first step beyond the single-particle inclusive level