The Erlang distribution is a direct result of the multi-source thermal model [20–22]. One or two-component Erlang distribution can describe the narrow or wide p _T spectra of particles, where the narrow (wide) p _T spectra refers to range less than a few GeV/c (more than 10 GeV/c) [22]. The multi-source thermal model assumes that multiple sources are formed in high-energy collisions. These sources can be nucleons or partons if we study the formation of nucleon clusters (nuclear fragments) or particles.

In this work, we assume that a few (n _s ) partons (partons-like) contribute to p _T of a given particle [22]. The contribution of the jth parton is assumed to be an exponential function with variable p _tj which depends on j, and average value ⟨p _t ⟩ which is independent of j. We have the normalized exponential function f p t j = 1 〈 p t 〉 exp − p t j 〈 p t 〉 . (1)

Here, ⟨p _t ⟩ represents the average contribution of participant partons to ⟨p _T ⟩ of the considered particles.

The contribution sum ( p t 1 + p t 2 + ⋯ + p t n s ) of n _s partons is p _T of a given particle. The result convoluting the contributions of n _s partons is the Erlang distribution. We have the Erlang p _T distribution to be f p T = 1 N d N d p T = p T n s − 1 n s − 1 ! 〈 p t 〉 n s exp − p T 〈 p t 〉 . (2)

Here N is the number of particles, and the form of (1/N)dN/dp _T results in the normalization of f(p _T ) to 1. In fact, the normalization of the Erlang distribution is naturally 1.

We would like to emphasize here the difference between “n _s ”, the number of partons and “N”, the number of particles. In γ*p collisions, if three quarks in the proton contributed to p _T , we have n _s = 3. If another q q ̄ pair also contributed to p _T , we have n _s = 3 + 2 = 5. Even in nucleus-nucleus collisions, the value of n _s is not large due to it being determined by the number of contributor partons in a nucleon-nucleon pair, but not collision system itself. This makes sense, in the Fock’s first two terms of the development of the wave function of the proton, as composed by 3 quarks and then 3 quarks plus a q q ̄ pair [41, 42]. As for N, its value may be small in small collision system or peripheral nucleus-nucleus collisions. The value of N may be very large in central nucleus-nucleus collisions at high energy.

As we know, both the average transverse momentum ⟨p _T ⟩ and initial-state temperature T _i [11–15] describe the excitation degree of the system. In particular, in the Erlang distribution, ⟨p _T ⟩ can be easily obtained by 〈 p T 〉 = ∫ 0 ∞ p T f p T d p T = n s 〈 p t 〉 , (3) where f (p _T ) is normalized to 1. Similarly, ⟨p _t ⟩ reflects the excitation degree of participant partons.

According to Refs. [43–45], with a color string percolation method [46], T _i can be regarded as T i = ⟨ p T 2 ⟩ 2 F ξ , (4) where ⟨ p T 2 ⟩ = ∫ 0 ∞ p T 2 f p T d p T (5) due to f (p _T ) is normalized to 1 and ⟨ p T 2 ⟩ is the root-mean-square of p _T . In Eq. 4, F(ξ) is the color suppression factor [46].

In the process of using color string method to obtain T _i in this work, only one string is used, i.e., F(ξ) = 1, in the formation of particle. Although there are probability to have any other strings, they do not affect noticeably T _i . If we consider other strings, according to Ref. [46], one has the minimum F(ξ) ≈ 0.6. This will cause the maximum increase of 29.1% in T _i . Considering the fraction of one string is very large, that of two strings is relative small, and that of multiple strings is very small, the increase in T _i will be much smaller than 29.1%.

In the center-of-mass reference frame, in two-body process 2 + 1 → 4 + 3 or two-body-like process of high-energy collisions, there are three Mandelstam variables defined based on the four-momenta of these particles. They have the forms to be s = − P 1 + P 2 2 = − P 3 + P 4 2 , (6) t = − P 1 − P 3 2 = − − P 2 + P 4 2 , (7) u = − P 1 − P 4 2 = − − P 2 + P 3 2 , (8) where P ₁, P ₂, P ₃, and P ₄ are four-momenta of particles 1 (target proton), 2 (incident γ*), 3 (emitted nucleon), and 4 (emitted meson), respectively. Here, we assume that particle 1 is incident along the Oz direction and particle 2 is incident along the opposite direction. After collisions, particle 3 is emitted with angle θ relative to the Oz direction and particle 4 is emitted along the opposite direction.

The three Mandelstam variables have different physical meaning. For instance, s refers to the center-of-mass energy, and −u is defined as the squared momentum transfer between particles 1 and 4. Here, selected variable −t (the squared momentum transfer between particles 1 and 3) is calculated to fit the experimental data. For convenience, we have | t | = | E 1 − E 3 2 − p ⃗ 1 − p ⃗ 3 2 | = | m 1 2 + m 3 2 − 2 E 1 p 3 T sin ⁡ θ 2 + m 3 2 + 2 E 1 2 − m 1 2 p 3 T tan ⁡ θ | , (9) where E ₁ and E ₃, p ⃗ 1 and p ⃗ 3 , as well as m ₁ and m ₃ are the energy, momentum, and rest mass of particles 1 and 3, respectively. In addition, p _3T referred to be perpendicular to the Oz direction component of the transverse momentum of particle 3, which obeys the Erlang distribution, that is Eq. 2 in which p _T = p _3T .

In this paper, the squared momentum transfer spectra of light meson at different squared photon virtuality Q ² and Bjorken variable x _B are fitted by calculated results with the Monte Carlo method. Here, Q ² is a reflection of hard scale of reaction [47–54]. The harder the reaction is, the higher the excitation degree is. In fact, Q ² is the absolute value of the squared mass of γ* (particle 2) that is exchanged between the scattered electron and the target proton (particle 1), and it effectively represents the transverse size of the probe [38]. In addition, −Q ² is also the squared momentum transfer to the target proton (particle 1) by the scattered electron [37].

As for the Bjorken variable x _B , it represents contrarily the momentum of particle 1. The lower the x _B is, the higher the momentum of particle 1 is. Generally, x B = Q 2 / ( 2 P 2 ⋅ − Q 2 ) ∝ Q [37]. In a symmetric frame, importing ξ′ as skewness, it is half of the longitudinal momentum fraction transferred to the struck parton. The skewness ξ′ can be used to express x _B approximately. That is x _B ≈ 2ξ′/(1 + ξ′) [37].

In the calculations of squared momentum transfer, the analytical expression of p _T distribution is difficult to be transformed to that of squared momentum transfer distribution directly by using Eq. 9. Alternatively, we may use the Monte Carlo method to transform p _T to squared momentum transfer. Let R _1,2 and r 1,2,3 , … , n s be random numbers distributed evenly in [0,1]. Then, many concrete transverse momentum p _3T satisfied with Eq. 2 and θ are produced. Other quantities such as E ₁, m ₁, and m ₃ in the equation are fixed, though E ₁ is treated as a parameter in the present work.

Generally, we may solve the equation ∫ 0 p T f p T ′ d p T ′ < R 1 < ∫ 0 p T + δ p T f p T ′ d p T ′ , (10) where δp _T is a small shift relative to p _T . Conveniently, there is a simpler expression due to Eqs 1, 2. In fact, solving the equation ∫ 0 p t j f p t j ′ d p t j ′ = r j j = 1,2,3 , … , n s , (11) we have p t j = − ⟨ p t ⟩ ln r j j = 1,2,3 , … , n s . (12)

The simpler expression is p T = ∑ j = 1 n s p t j = − 〈 p t 〉 ∑ j = 1 n s ln r j = − 〈 p t 〉 ln ∏ j = 1 n s r j . (13)

The distribution of θ satisfies with the half-sine function f θ θ = 1 2 sin ⁡ θ (14) which is obtained under the assumption of isotropic emission in the source’s rest frame. Solving the equation ∫ 0 θ f θ θ ′ d θ ′ = R 2 , (15) we have θ = 2 ⁡ arcsin R 2 (16) which is needed in the calculations.

We have check the consistency and correctness of the above expressions in the Monte Carlo method in terms of illustration which is not presented here. After obtaining concrete values of p _3T and θ, and using E ₁, m ₁, and m ₃, the value of |t| can be obtained from Eq. 9. Through repeating the calculations many times, the distribution of |t| is obtained statistically. Based on the method of least squares, the parameter ⟨p _t ⟩ and n _s are extracted naturally. Meanwhile, T _i can be obtained from Eq. 4 and ⟨ p T ⟩ ( ⟨ p T 2 ⟩ ) can be obtained from Eq. 3 [(5)] or from the statistics. The errors of parameters are obtained by the general method of statistical analysis.

It should be noted that the above Monte Carlo calculation is only performed in the transformation from transverse momentum to |t|, in which the physics process such as the radiative corrections for reactions induced by electrons has been taken into account naturally. In fact, the effects of the mentioned process and all other processes are included in the Erlang distribution which is a result of multi-factor interactions. In other words, the Monte Carlo calculation used here is not a simulation for the system evolution from initial to final stages, but the numerical transformation in the final stage.

Figure 1 shows the differential cross-section, dσ/d|t|, in squared momentum transfer |t| of γ*p → π ⁰ p process produced in 5.75 GeV electron beam induced collisions in a 2.5 cm long liquid-hydrogen target (ep collisions at beam energy of 5.75 GeV) in different ranges of squared photon virtuality, 1.0 < Q ² < 1.5, 1.5 < Q ² < 2.0, 2.0 < Q ² < 2.5, 2.5 < Q ² < 3.0, 3.0 < Q ² < 3.5, 3.5 < Q ² < 4.0, and 4.0 < Q ² < 4.6 GeV², from bottom to up sub-panels, as well as in different ranges of Bjorken variable, 0.10 < x _B < 0.15, 0.15 < x _B < 0.20, 0.20 < x _B < 0.25, 0.25 < x _B < 0.30, 0.30 < x _B < 0.38, 0.38 < x _B < 0.48, and 0.48 < x B < 0.58, from left to right sub-panels. The sample at the top-left sub-panel shows repeatedly the result in the range of squared photon virtuality, 1.5 < Q ² < 2.0 GeV², and the range of Bjorken variable, 0.20 < x _B < 0.25, as an example. The symbols represent the experimental data measured by the CLAS Collaboration [37] and the curves are the statistical results of squared momentum transfer |t|.

In Eq. 9, p _3T satisfies the Erlang distribution (Eq. 2) and we obtain it by the Monte Carlo method (Eq. 13). Then, the squared momentum transfer |t| is obtained statistically. In the fitting process, two main parameters, i.e., the average transverse momentum ⟨p _t ⟩ contributed by each participant parton and the number n _s of participant partons are extracted naturally. To obtain a better fit result, E ₁ is extracted as an insensitive parameter. In addition, a non-free parameter is the normalization constant σ ₀. The values of parameters with selection condition (Q ² and x _B ), χ ², and the number of degree of freedom (ndof) are listed in Table 1, where the number of parameters is always 4 which includes E ₁, ⟨p _t ⟩, n _s , and σ ₀. In the case of ndof being less than or equal to the number of parameters, we obtain the curve from a “prediction” or extrapolation based on other reasonable fits in which the tendency of parameters is available. Meanwhile, in these cases, the number of points (nop) is given in a bracket to replace ndof in the table. One can see that the values of χ ² are small in most cases, though the (necessary) dense log scale is not easy to judge. The model results are in agreement with the experimental data. From the values of parameters, the average transverse momentum ⟨p _T ⟩ and initial temperature T _i are obtained naturally.

Figure 2 presents the differential cross-section, dσ/d|t|, in |t| of γ ^∗ p → π ⁺ n process produced in ep collisions at beam energy of 6 GeV at different squared photon virtuality, Q ² = 1.75, 2.05, 2.35, 2.65, 2.95, 3.35, 3.85, and 4.35 GeV², as well as at different Bjorken variable, x _B = 0.25, 0.31, 0.37, 0.43, 0.49, and 0.55. The symbols represent the experimental data measured by the CLAS Collaboration [38]. As those in Figure 1, the curves in Figure 2 are also the statistical results of |t| in which p _3T satisfies the Erlang distribution and is obtained by the Monte Carlo method. The values of parameters with selection condition (Q ² and x _B ), χ ², and ndof are listed in Table 1. One can see that the model results are in agreement with the experimental data.

Figure 3 displays the differential cross-section, dσ/d|t|, in |t| of γ ^∗ p → ηp process produced in ep collisions at beam energy of 5.75 GeV in different Q ² and x _B ranges shown in the panel. As an example, the sample at the top-left sub-panel shows repeatedly the result in 1.5 < Q ² < 2.0 GeV² and 0.20 < x _B < 0.25. The symbols represent the experimental data measured by the CLAS Collaboration [39]. The curves are the statistical results of |t| in which p _3T satisfies the Erlang distribution and is obtained by the Monte Carlo method. The values of parameters with selection condition (Q ² and x _B ), χ ², and ndof are listed in Table 2. One can see that the model results are in agreement with the experimental data.

Similar to Figures 1–3, Figure 4 presents the differential cross-section, dσ/d|t|, in |t| of γ ^∗ p → ρ ⁰ p process produced in ep collisions at beam energy of 5.754 GeV in different Q ² and x _B ranges shown in the panel. As an example, the sample at the top-left sub-panel shows repeatedly the result in 2.8 < Q ² < 3.1 GeV² and 0.40 < x _B < 0.46 range. The symbols represent the experimental data measured by the CLAS Collaboration [40]. The curves are the statistical results of |t| in which p _3T satisfies the Erlang distribution and is obtained by the Monte Carlo method. The values of parameters with selection condition (Q ² and x _B ), χ ², and ndof are listed in Table 2. One can see that the model results are in agreement with the experimental data.

In Figures 1–4, the cross-sections for π ⁰ p, π ⁺ n, ηp, and ρ ⁰ p are fitted to show some differences in concrete values and parameters, and some common features among them in the tendency of curves also appear. This is caused by the fact that different channels have different fraction ratios, and all of them are from the same ep collisions, though the collision energies are slightly different.

The dependences of ⟨p _T ⟩ (A, C, E, G) and T _i (B, D, F, H) on Q ² in γ*p collisions with different emitted channels [π ⁰ p (A, B), π ⁺ n (C, D), ηp (E, F), and ρ ⁰ p (G, H)] are shown in Figure 5, where ⟨p _T ⟩ = n _s ⟨p _t ⟩ due to Tables 1, 2 and the values of T _i are from Tables 1, 2. Different symbols represent the results for different x _B . One can see that ⟨p _T ⟩ and T _i increase generally with an increase in Q ². Because Q ² represents the hard scale (violent degree) of collisions and a harder scale results in a higher excitation degree, it is natural that larger ⟨p _T ⟩ and T _i appear at higher Q ².

Figure 6 is similar to Figure 5, but it shows the dependences of ⟨p _T ⟩ (A, C, E, G) and T _i (B, D, F, H) on x _B in γ*p collisions with emitted channels π ⁰ p (A, B), π ⁺ n (C, D), ηp (E, F), and ρ ⁰ p (G, H). Different symbols represent the results for different Q ². One can see that ⟨p _T ⟩ and T _i increase generally with an increase in x _B . Because x _B ∝ Q, we may think that x _B also represents the hard scale of collisions and a harder scale results in a higher excitation degree. It is understandable that larger ⟨p _T ⟩ and T _i appear at higher x _B .

In addition, x _B also represents the longitudinal momentum fraction transferred to the struck parton. In the considered γ*p → meson + nucleon process in ep collisions at given energy, the larger x _B means the larger longitudinal momentum transfer to the struck parton or the system, and hence the more energy deposited to the system. The system naturally stays at higher excitation degree. As a result, larger ⟨p _T ⟩ and T _i are observed.

Generally, ⟨p _T ⟩ > T _i ≥ T _ch ≥ T ₀. If the evolution time of the system is 0, that is if the initial-state, chemical freeze-out, and kinetic freeze-out happen simultaneously, we have T _i = T _ch = T ₀. If the evolution time is not negligible, we have T _i > T _ch > T ₀. The difference between ⟨p _T ⟩ and temperature is explained as the contribution of flow effect. According to Ref. [55], in the final-state, the expected real T ₀ ≈ ⟨p _T ⟩/3.07. Then, we have the contribution of flow effect to be ⟨p _T ⟩ −T ₀ ≈ 2.07⟨p _T ⟩/3.07. One can see that the flow effect contributes largely to ⟨p _T ⟩. It is expected that the contribution of flow effect increases with the increase of evolution time, if ⟨p _T ⟩ is fixed from the initial- to final-states.

From Tables 1, 2, we note that the values of n _s are 3–5 for different channels. As the number of participant partons, n _s is constrained to be integer with uncertainty of 0. For a given channel, n _s is independent of Q ² and x _B in most cases. The channel independent n _s renders that the number of participant partons is not too small or big. The number of struck parton(s) is usually regarded as 1 or 2, which is very small. The struck parton(s) and the partons around the struck parton(s) are participant partons. The partons far away from the struck parton(s) are remainder or spectator partons.

Before summary and conclusion, we would like to point out that the discussion about the temperature and flow in this paper is applicable. Although the multiplicity in ep collisions at a few GeV is very limited and the final particles are in a state far from thermal equilibrium, we may use the grand canonical ensemble for lots of events in which the number of total particles is very large and the whole system is in a homogeneous and equilibrium state. Therefore, the temperature used in this paper is comparable to the freeze-out temperatures used in nucleus-nucleus collisions. Of course, we may also regarded the temperature used here as a fitting parameter if necessary.

The initial-temperature T _i is extracted from the root-mean-square of p _T , which is independent of model, though the relation between T _i and ⟨ p T 2 ⟩ is from the color string percolation method [43–46]. As deep inelastic scattering, ep collisions are head-on collisions, and may be harder than nucleus-nucleus collisions at similar energy per nucleon due to the fact that some non-head-on nucleon-nucleon collisions exist in the later. As a hybrid state of head-on and non-head-on nucleon-nucleon collisions, nucleus-nucleus collisions may be weaker than head-on ep collisions. In addition, cold spectator nuclear effect also causes the temperature in nucleus-nucleus collisions to reduce. This renders that T _i obtained in this paper is higher than that in nucleus-nucleus collisions.

It should be emphasized that the parameter T _i reflects the violent degree of collisions. To our knowledge, other groups and other studies where T _i is extracted for hadronic collisions is not available at present, though T _i for nucleus-nucleus collisions is available. In terms of T _i , Erlang distribution, and Monte Carlo calculation, the present work has proposed an alternative method to describe light meson electroproduction data obtained with the JLab-CLAS facility. Typically those data are interpreted in terms of handbag diagram within the formalism of generalized parton distributions, whereas here statistical methods, that were developed for high-energy nucleus-nucleus collisions, are applied. At least, the present work has significance in the application of statistical methods.