Odderon and proton substructure from a model-independent Levy imaging of elastic $pp$ and $p\bar p$ collisions

We describe a new and model-independent L\'evy imaging method of quality fits to the published datasets and reconstruct the amplitude of high-energy $pp$ and $p\bar p$ elastic scattering processes. This method allows us to determine the excitation function of the shadow profile $P(b)$, the elastic slope $B(t)$ and the nuclear phase $\phi(t)$ functions of $pp$ and $p\bar p$ collisions directly from the data. Surprisingly, notable qualitative differences in $B(t)$ for $pp$ and for $p\bar p$ collisions point towards an Odderon effect. As a by-product, we clearly identify the proton substructure with two different sizes at the ISR and LHC energies, that has striking similarity to a dressed quark (at the ISR) and a dressed diquark (at the LHC). We present model-independent results for the corresponding sizes and cross-sections for such a substructure for the existing data at different energies.

"What we pay attention to is largely determined by our expectations of what should be present." Christopher Chabris

I. INTRODUCTION
The TOTEM Collaboration at the Large Hadron Collider (LHC) at CERN has released recently two new data sets from the first measurements of the total, elastic and differential cross sections, as well as the ρ-parameter, of elastic pp collisions at the currently highest available energy of √ s = 13 TeV [1,2]. Taken together, these papers indicate the influence of an odd-under-crossing (or C-odd) contribution to the elastic scattering amplitude, traditionally called the Odderon [3], or in more modern language of Quantum Chromo Dynamics (QCD), an odd-gluon (predominantly, three-gluon) bound state, a quarkless excitation sometimes also referred to as a vector glueball. As far as we know, the properties like the √ s and t dependences of the differential cross-section of an Odderon contribution at the LHC energies were determined from pp and pp collisions first in Ref. [4].
The new TOTEM preliminary results [1,2] generated a burst of high-level and vigorous theoretical debate about the correct interpretation of these data, see Refs. [5]- [18]. All possible extreme views were present among these papers, including claims for a maximal Odderon effect [14] and claims of lack of any significant Odderon effects, see Refs. [10,18].
In this paper, we investigate earlier, published data and a recently released, new TOTEM preliminary data set [19], to look for the Odderon effects in elastic pp collisions, and for the interpretation of the data, using a new and model-independent imaging technique, the Lévy series. Our analysis is based on the TOTEM preliminary data as presented by F. Nemes for the TOTEM Collaboration in Ref. [19].
We find that clear, but indirect signals of Odderon effects are present in the differential cross-section of elastic pp and pp scattering, as indicated by the difference of the fourmomentum transfer dependence of the elastic slope B(t) for pp and for pp collisions. A less evident but clear difference is also identified between the nuclear phase φ(t) of pp and pp collisions in the TeV energy range.
Although our analysis was motivated by the search for Odderon effects, our most surprising result is that we find a clear-cut evidence for a proton substructure having two distinct sizes in the GeV and TeV energy ranges, respectively. We estimate these sizes and the corresponding contributions to the total pp cross-section at the ISR and Tevatron/LHC energies.
The structure of the presentation is as follows. In Section II, we present the modelindependent imaging approach of Lévy expansions, in Section III, we apply this method in a comprehensive analysis of elastic pp and pp collisions, and in Section IV, we discuss and interpret the results, as well as present the indirect signals of Odderon effects in elastic scattering in the TeV energy range. Finally, we summarize and conclude in Section V. This manuscript is closed by four Appendices. Appendix A details the results of the Lévy expansion fits to elastic pp scattering data for the whole available region of t. Appendix B shows the individual fits to elastic pp scattering data. Appendix C describes fits to the tails (i.e. in the large-|t| region just after the dip and the bump structure) of the differential cross-section of elastic pp scattering for all the experimentally accessed energies. These fits indicate an evidence for a proton substructure. Appendix D details Lévy fits to the cone (small-|t|) regions of elastic pp scattering, indicating that the proton size grows self-similarly in the GeV energy range, but at the LHC energies of √ s = 7 and 13 TeV, something drastically changes, such that the protons keep on growing but their shape also becomes significantly different from their shape in the ISR energy range of √ s = 23.5 -62.5 GeV.

II. MODEL-INDEPENDENT LÉVY ANALYSIS OF ELASTIC pp AND pp SCAT-TERING
We describe a model-independent Lévy series, that is a generalization of the Laguerre, Edgeworth [20] and Lévy expansion [21] methods proposed to analyze nearly Lévy stable source distributions in the field of particle correlations and femtoscopy. The key point of this method is to have a look at the data, guess their approximate shape (for example, Gaussian, exponential or Lévy-stable shape) and then develop a systematic method to characterize the deviations from the approximate shape with the help of a dimensionless variable denoted in this paper by z, and a complete orthonormal set of polynomials that are orthogonal with respect to the assumed zeroth order shape or weight function w(z). We recommend Ref. [20] for detailed discussions and for the convergence criteria of such expansions given in general terms. This way one may find the minimal number of necessary expansion coefficients.
For example, if the data follow the guessed approximate shape precisely, that can be clearly demonstrated by fitting the series to the data and finding that all the expansion coefficients that measure deviations from the zeroth-order shape are within errors consistent with zero. Indeed, the PHENIX Collaboration analyzed recently Bose-Einstein correlations of √ s N N = 200 GeV Au+Au collisions [22], and found that they are well described by the Lévy stable source distributions. The accuracy of the Lévy description was tested by a Lévy expansion of the Bose-Einstein correlation functions, as proposed in Ref. [23], to find that the first-order deviations from the Lévy stable source distributions within errors are consistent with zero [22].
The data analysis method of Ref. [23] was developed first for the functions that may oscillate between positive and negative values. However, the differential cross-section of elastic scattering is measured as an angular-dependent hit distribution, so it is a positive definite function, related to the modulus square of the elastic scattering amplitude. Hence, we describe this expansion at the amplitude level, with complex expansion coefficients, and then take the modulus square of such an amplitude to get a positive-definite form.
The differential cross-section of elastic scattering at high energies is measured as a function of the four-momentum transfer squared, the Mandelstam variable t = (p 1 − p 3 ) 2 < 0. A differential cross-section of elastic scattering starts at the optical point at t = 0, decreases quickly and nearly exponentially in |t|, as typically characterized by an exponential slope parameter B, as follows: dσ dt = A exp(−B|t|). The region, where such a behaviour is approximately valid, is called a diffractive cone, and such a featureless exponential decay corresponds to a nearly Gaussian distribution of the centers of elastic scattering [24]. This region is followed by (one or more) alternating diffractive minima and maxima, finally at high four-momentum transfers, a diffractive tail might be seen as well. For more details and for an introduction and review of diffraction before the start of the LHC measurements, see Ref. [24].
We know from the TOTEM analysis of √ s = 8 TeV elastic pp scattering data, that the differential cross-section in the diffraction cone, at |t| values below the diffraction minimum, deviates significantly from an exponential shape [25]. This deviation is a subtle, but a more than 7σ effect [25]. Using this knowledge, we assume that the elastic scattering amplitude is nearly (but not completely) exponential in |t|, i.e. we assume that it is approximately described by a (Fourier-transformed) Lévy stable source distribution: The conventional exponential behaviour corresponds to the α = 1 special case. This way the deviation from the exponential behaviour can be quantified by a single parameter. If the value of the exponent α is significantly different from unity, it evidences a non-exponential behaviour of the differential cross-section of elastic scattering. Later on, we shall see that this is a very fortunate approach, and it connects the imaging of the differential cross-sections of high-energy pp and pp collisions with the Lévy stable source distributions that are ubiquitous in Nature [26].
We also know that the differential cross-section of high-energy pp and pp collisions has a diffractive minimum, followed by a second maximum, that is followed by an extended tail. Thus, the behaviour of the differential cross-sections at large |t| has structures that deviate from a simple Lévy stable source. In this paper, we attempt to characterize these structures with an orthonormalized Lévy expansion. This way we obtain a new imaging method, that we describe in detail below, and apply it to reconstruct the shadow (inelasticity) profile functions, the t-dependent slope parameters and the t-dependent nuclear phases in highenergy pp as well as pp collisions.
These physical and mathematical assumptions result in the following formulae for the differential cross-section of elastic pp and pp collisions: 3) c j = a j + ib j , (2.4) where w(z|α) is the Lévy weight function, and a dimensionless variable z is introduced as the magnitude of the four-momentum transfer squared |t| multiplied by a Lévy scale parameter R squared, where R is measured in units of fm (natural units c = = 1 are adopted here and below). Note that the w(z|α) is also called the stretched exponential distribution, which, for α = 1 limiting case, corresponds to the exponential distribution. This shape actually corresponds to a Fourier-transformed and modulus squared symmetric Lévy-stable source distribution centered on zero. The orthonormal set of Lévy polynomials, defined below, are denoted as l j (z|α). The complex expansion coefficients are c j , with a j and b j standing for the real and the imaginary parts of c j , respectively. In the forthcoming, we shall refer to the zeroth-order (c i = 0) Lévy expansion simply as a Lévy fit. Let us clarify here that these so-called Lévy fits actually correspond to the Fouriertransformed and modulus squared, symmetric Lévy stable source distributions t el (b), as approximations to the shape of the amplitude of elastic pp scattering. This elastic amplitude will be introduced and detailed in Subsection II A. For more details on the application of Lévy stable source distributions [27,28] in particle correlations and femtoscopic measurements, we recommend Refs. [21-23, 29, 30].
The expansion (2.1) is expected to converge to nearly Lévy shaped data, if the order of the series n is chosen to be sufficiently large, i.e. n → ∞. In practice, however, third-order (n = 3) Lévy series already converged to the data measured at √ s < 1 TeV, with confidence levels (i.e. the probability that the Lévy series or expansion of the elastic scattering amplitude converges to the differential cross-section under investigation) corresponding to a statistically acceptable description. In order to gain a statistically marginal or acceptable description of the high precision TOTEM data at 7 TeV and preliminary data at 13 TeV, we had to go to the fourth-order Lévy series, n = 4, in these two cases. For reasons of consistency, and in order to eliminate fitting artefacts that may show up if one compares different orders of the Lévy expansions with one another, we decided to re-fit all the pp elastic scattering data with the fourth-order Lévy expansion and to show only these results. However, when applying a similar procedure to pp elastic scattering, it turned out that the range of the data around the dip position was too limited in this case, and the fourth-order Lévy expansion terms could not be determined in a reliable and reasonable manner. So we decided to show only the fit results of the second-and third-order Lévy expansions for the pp elastic scattering data.
As we explicitly demonstrate in Appendices C and D, in certain limited intervals of |t|, Lévy fits without correction terms (i.e. for c i = 0, i ≥ 1) provided statistically not unacceptable, but in contrast, rather good quality fits and the corresponding confidence levels. These results suggest that the Lévy series is a reasonable representation of the scattering amplitude of elastic pp and pp collisions.
The first four orthogonal (but not yet normalized) Lévy polynomials denoted as L i (z | α) are found as follows and Euler's gamma function is defined as The normalization of these Lévy polynomials is straightforwardly expressed as follows: where D 0 = 1 and, in general, D j ≡ D j (α) stands for the Gram-determinant of order j, 14) These normalized Lévy polynomials l j (z|α) are, as far as we know, newly introduced in this work, while the unnormalized Lévy polynomials L j (z|α) were introduced earlier in Ref. [21].
The orthonormality of {l j (z|α)} ∞ j=0 with respect to a Lévy or stretched exponential weight is expressed by the following relation: (2.16) The first few of these orthonormal Lévy polynomials are illustrated in Fig. 1 for a specific value of α = 0.9. More details of them, in particular the α = 2 Gaussian, the α = 1 Laguerre special cases and their explicit forms are being described in Refs. [29,30]. Once we have a statistically acceptable description of the differential cross-section of elastic pp and pp collisions, we can build the elastic scattering amplitude with the help of the Lévy imaging method, and we can then compare the resulting shadow profiles of pp and pp collisions without any model-dependent assumptions. Similarly, we can extract the t-dependence of the nuclear slope B(t) directly from the data, as well as compare its behaviour for pp collisions with that of pp collisions. The same is true for the nuclear phase and for the t-dependent ρ-parameter as well.
In fact, in our analysis we rely only on the convergence of the Lévy series (or Lévy expansion), that we have tested by the usual χ 2 -optimization methods with the CERN Minuit package and by evaluating the confidence level.
Based on our experience with extracting B(t), ρ(t) and the shadow profile P (b) from pp and pp elastic collision data, we can definitely state that the precise reproduction of the measured data points, with a statistically acceptable confidence level of CL > 0.1% is a necessary condition for interpreting our fit results. We have achieved such good quality fits in many cases of the analysis of the published, final data, except for the 7 TeV pp elastic scattering data, where we reached a marginal confidence level of CL ≈ 0.02%, as indicated in Fig. 2. After scrutinizing this fit, presented in Fig. 2, we decided to interpret the parameters of this fit as well. But in principle, in order to get the final errors of our parameters, we may need to repeat the analysis by taking into account the full covariance matrix. Although the fit quality is marginal, CL ≈ 0.02%, the fitted curve follows the data so closely that we decided to interpret the fit parameters, noting that the errors on the best values of the parameters are likely underestimated.
In case of the TOTEM preliminary 13 TeV data set, we also warn the reader that these data points and their errors are still in a preliminary phase, so we have determined the best preliminary values of the parameters of the Lévy series from the minimum of the χ 2distribution. The preliminary value of the confidence level of fits to the TOTEM preliminary data at √ s = 13 TeV CL= 2%, as indicated in Fig. 3, satisfies the criteria for good quality fits.   4 summarizes the fits with a fourth-order Lévy expansion to all the differential crosssection measurements of elastic proton-proton collisions from √ s = 23.5 GeV up to 13 TeV. These fits are detailed in Appendix A, where the fits to each dataset are shown in detail, with the fit parameters and confidence levels (or p-values) are listed on the corresponding plots, together with the total and elastic cross-sections as well as the ρ(t = 0) values that are calculated from the fit parameters.
Similarly, Fig. 5 summarizes the fits with a third-order Lévy expansion to all the differential cross-section measurements of elastic proton-antiproton collisions from √ s = 53 GeV up to 1.96 TeV. The fits converged, error matrix was accurate and CL ≥ 0.1 % for these fits, that were obtained with fixed α = 0.9. These fits are described in greater details in Appendix B, where the fit parameters and confidence levels (or p-values) are listed on the corresponding plots, together with the total and elastic cross-sections as well as the ρ(t = 0) values that are calculated from the fit parameters. Fig. 6 represents the summary plot of the Lévy fits, dσ/dt = A exp(−(R 2 |t|) α ) , that correspond to the zeroth-order of the Lévy expansion detailed in this manuscript, to the tails of the elastic pp scattering data at ISR and LHC energies from 13 TeV, with α = 0.9 fixed. These results are detailed in Appendix C, and are explained in terms of the newly identified proton substructure in Subsection III D. Namely, rather elegantly and clearly, a smaller substructure is seen in the ISR energy range, that is invariant for the (relatively small) change of √ s, as evidenced by the dashed lines that (except an overall normalization factor) follow the same curves. It is apparent from the visualization of Fig. 6, that at √ s = 7 and 13 TeV, the slope of these dashed lines changes dramatically: a proton substructure of a different size is found that is apparently (within the errors) the same both at 7 and at 13 TeV. Fig. 7 represents a summary plot of the Lévy fits, dσ/dt = A exp(−(R 2 |t|) α ) to the cone or low-|t| part of the elastic pp scattering data at ISR and LHC energies from √ s = 23.5 GeV up to 13 TeV, with α = 0.9 fixed. All the low-t ISR fits are successful with the Levy leading order result where there are sufficient data sets at low-t, except those data at √ s = 44.7 and 52.8 GeV which exhibit a gap in the measured dataset at low values of |t|. These results are detailed in Appendix D. The gradual steepening of the slope of the fitted curves indicate, that the Lévy scale R, characterizing the overall size of the proton, The fits converged, error matrix was accurate and CL ≥ 0.1 % for these fits, that were obtained with fixed α = 0.9. These fits are detailed in Appendix B.
was increasing monotonically with increasing energy √ s. This behaviour can be explained in terms of the proton size growing self-similarly in the ISR energy range, as evident also from the excitation function of the shadow-profiles at 23.5 ≤ √ s ≤ 62.5 GeV. This effect is discussed in greater details in Subsection III B.
However, at √ s = 7 and 13 TeV, the leading-order Levy fits for fixed α = 0.9 have failed, indicating that not only the size of the proton increases with an increase of collision energies, but also the shape of the protons changes. A successful fit in this case is possible only if the first few data points are taken only into account. Note, that any data in the Coulomb-nuclear interference (CNI) region, i.e. at |t| < 0.01 GeV 2 , were not included in these fits, thus one may neglect any possible CNI effect when interpreting the results of the Levy expansion.
The observation that the impact parameter b dependent shadow profile function P (b) at 7 and 13 TeV deviates from a Levy shape even in the leading order confirms that the shape of protons changes at LHC energies in a way which is different from that at ISR energies. These results, detailed in Subsection III B, can be explained in terms of the evolving shape of the shadow profiles at 7 and 13 TeV which exhibits plateaux near b = 0 (which were not seen at lower energies). As shown in Subsection III B, one observes a saturation of the shadow profile functions P (b) ≈ 1 in the b ≤ 0.4 − 0.5 fm region at 7 and 13 TeV, respectively.
A. Differential, total and elastic cross-sections The conventional form of the elastic differential cross section  provides us with the key expression for the complex-valued elastic scattering amplitude T el (∆) in terms of a Lévy series Then, according to the optical theorem, the total cross section is found as 20) while the ratio of the real to imaginary parts of the elastic amplitude is known as the ρ-parameter, in consistency with the traditional form of the forward elastic differential cross section The t-dependent elastic slope B(t) is traditionally defined as There is a great current interest in the value of this function at t = 0 at LHC energies. Traditionally, the elastic slope is determined as B = B(t = 0). Let us note that this requires an extrapolation of the measured differential cross-sections to the t = 0 optical point. Frequently, an exponential approximation is applied, however, at the very low-|t| region the CNI complicates such an extrapolation as well [25]. At this point, it is important to emphasize, that the Lévy series utilized in this paper to represent the elastic scattering amplitude is not an analytic function at t = 0 if α < 1, hence formally our expressions for B may not exist, as B(t) is well-defined only for |t| > 0 in this case. However, if α = 1, the cone region decreases exponentially and the elastic scattering amplitude becomes an analytic function at t = 0. Hence, it is very important to determine the value of α precisely from the analysis of the elastic differential cross-section data.
Note that B is related to the root-mean-square (RMS) radius of the impact-parameter bdependent elastic amplitude t el (b). It is well known that for Lévy-stable source distributions, that are our zeroth-order choices for the impact parameter dependent elastic amplitudes, the RMS of the source is divergent, if the Lévy index of stability α L < 2 [27], with the exception of the Gaussian case, corresponding to α L = 2, when the RMS of the source is finite. Due to the importance of this point, we have dedicated Subsection II E to compare the differential cross-section of Lévy-stable sources with Gaussian sources and dedicated Appendices C and D to investigate if such a non-analytic model at t = 0 as the Lévy-stable source distribution t el (b) can describe reasonably well the pp elastic differential cross-section data in limited kinematic regions. Actually, we find that this is indeed a very good approximation to the data at low values of |t| in the ISR energy range of 23.5 ≤ √ s ≤ 62.5 GeV, as explicitly demonstrated with α = 0.9 fixed fits in Appendix D. Let us note here that the cone region of TOTEM data at √ s = 7 and 13 TeV can also be described in the same cone (or low-|t|) region only approximately, if the parameter α is released, corresponding to a change of the proton shape in the TeV energy range.
For more detailed examples and for the introduction of Lévy stable source distributions to femtoscopy in high-energy particle and nuclear physics, with an emphasis at their nonanalytic nature of their Fourier-transform, we recommend Ref. [31].
C. Four-momentum transfer dependent nuclear phase φ(t) The nuclear phase φ(t) can also be introduced at this point. A conventional definition of this phase φ(t) is referring to it as the phase of the elastic scattering amplitude in the complex plane, as follows: T el (t) = |T el (t)| exp(iφ(t)) .

(2.24)
An alternative definition was used recently by the TOTEM Collaboration, that related φ(t) to the t-dependent ρ(t) parameter as These definitions are equivalent if the nuclear phase satisfies 0 ≤ φ(t) ≤ π. The mathematical definition of the inverse tangent function arctan(x) is as follows. For a real number x, θ = arctan(x) represents the radian angle measure θ with − π 2 < θ < π 2 such that tan(θ) = x. Thus by definition, φ 2 (t) of eq. (2.25) satisfies 0 ≤ φ 2 (t) < π : it stands for the principal value of the nuclear phase φ(t). Given that for complex arguments, arctan(z) has branch cut discontinuities on the complex plane the principal value and the continuous definitions of the nuclear phase with the help of Eqs. (2.25) and (2.24) are, in general, inequivalent.
It is also interesting to evaluate the total elastic cross-section, σ el . In this calculation, we utilize the orthonormality of the Lévy polynomials l n (z|α) to obtain One may argue that the overall phase χ(z) of the elastic amplitude (2.18) is not constrained by the fits of the elastic cross section data. Indeed, an overall phase factor may in principle cause a redefinition of the complex expansion coefficients in the amplitude as such that the new coefficientsc i (z) = exp(iχ(z))c i could not be uniquely determined by the fits with the Lévy expansion. However, under certain additional assumptions, the Lévy expansion and the reconstructed phase is unique, as outlined below. By means of the optical theorem, the measurable total cross section uniquely determines the imaginary part of the elastic amplitude at t = 0, as indicated by Eq. (2.20). Measurements of the nuclear phase at t = 0 use the CNI region to determine ρ(t = 0) with typical values at LHC of the order of 0.1. This implies that not only the imaginary but also the real part of the forward scattering amplitude can also be uniquely determined at the optical point, t = 0, not allowing for an arbitrary phase at z = 0: the value of φ(t = 0) is fixed by measurements unambiguously. One can also determine the nuclear phase at the dip position t dip , where the differential cross-section of elastic proton-proton collisions has a diffractive minimum. At this point, the imaginary part of the elastic scattering amplitude approximately vanishes, hence the differential cross-section measures the real part of the forward scattering amplitude and the nuclear phase is an integer multiple of π. Under the additional assumption that the nuclear phase φ(t) is an analytic function of t, it can be continued with the help of our expansion method from its known value to arbitrary values of 0 ≤ −t < ∞.
In the subsequent applications, the nuclear phase φ(t) can thus be uniquely determined, if a well-defined diffractive minimum is included to the measured region, due to the following reasons: The arbitrary phase χ(z), based on the construction of Eq. (2.18), has to have a vanishing initial value both at small z, i.e. χ(z) → 0 for z → 0, as well as at the value of z dip that corresponds to t dip .
If a non-vanishing χ(z) emerges at z > 0, that would apparently destroy the orthonormality of the Levy terms (2.16) in the expansion for the amplitude, given that (2.27) which holds, by construction, for χ(z) = 0. Actually, the condition of the orthonormality may only be restored for functions that satisfy two conditions: cos(χ(z)) = 1, and χ(z = 0) = 0. There are infinitely many piecewise continuous functions that satisfy both requirements, by jumping between 0 and 2π at various values of z, however, the requirement that the arbitrary phase χ(z) is a continuous function of its argument z, uniquely determines that χ(z) ≡ 0. One may also observe that the overall normalization coefficient A and the zeroth order expansion coefficient a 0 + ib 0 would appear in the fits with Eq. (2.1) in a product form only. By absorbing these possible zeroth order constants into the overall normalization constant A, the zeroth order expansion simplifies as A exp(−z α ), that can be uniquely fitted to the data points. Furthermore all the fit parameters of Eq. (2.1) can be uniquely determined under these above listed conditions, that imply that χ(z = 0) ≡ 0. This proof indicates, that with the help of the Lévy expansion method, the nuclear phase of the elastic scattering amplitude can be unambiguously determined, if this nuclear phase is measurable at t = 0 (as evident from total cross-section and ρ 0 measurements at t = 0) and if it is assumed to be a continuous function of its argument t (apart from branch cuts, where uniqueness requires that we specify on which branch we define the value of the nuclear phase).
A posteriori, our method is validated by the excellent reproduction of the ρ 0 value measured in elastic pp collisions at √ s = 13 TeV by the TOTEM Collaboration [2] using the data in the CNI region. Fits detailed in Appendices A and B indicate that indeed the 4th and 3rd order Lévy expansions reproduce well measured values of ρ 0 not only at √ s = 13 TeV but at lower energies as well, if the fit quality is satisfactory. Results summarized in Appendices C and D also indicate that the zeroth order Lévy fits are not suitable to determine the nuclear phase and ρ 0 : one has to measure an interference to extract information about the phase, which interference is natural to find in the dip-bump region of elastic protonproton (or, proton-antiproton) reactions. Note that this proof however does not extend to the investigation of the convergence properties of the φ(t) measurements. We recommend to investigate the stability of the reconstructed φ(t) by fitting the data with higher and higher order Lévy expansions and to look for the numerical convergence of the reconstructed φ(t) functions and to determine the domain of convergence also numerically.

D. Shadow profile functions
Turning to the impact parameter space, we get This Fourier-transformed elastic amplitude t el (b) can be represented in the eikonal form where Ω(b) is the so-called opacity function, which is in general complex. Thus, a statistically acceptable description of the elastic scattering data provides us with a direct access to the opacity Ω(b) (known also as the eikonal function) and, in particular, to the shadow profile function defined as

E. A simple example -Gaussian versus Lévy stable source distributions
To gain intuition about the meaning of the characteristic Levy scale parameter of the proton R, before we go deeper to the data analysis, let us consider first the usual α = 1 case, neglecting all but the leading order (unity) term in the series that defines the Lévy expansion of the differential cross-section in Eq. (2.1).
In this Gaussian case, the differential cross section is apparently a structureless exponential in t, that corresponds to a Gaussian parametrization of the elastic scattering amplitude, based on Eq. (2.17): The Gaussian distributions correspond to central limit theorems, when several random elementary processes are convoluted to yield the final distribution. The Gaussian appears as a limiting distribution, if the elementary processes have finite means and variances, regardless of further details of the elementary probability distributions. Generalized central limit theorems describe limiting distributions for a large number of elementary processes, when the resulting elementary distributions have infinite means or variances. In these cases a limiting distribution exists, that remains stable for adding one more random elementary process. Due to this reason, such distributions are called stable, or, Lévy-stable distributions. They are denoted by S n (x|α L , β, γ, δ) where x is the variable of the distribution, α L stands for the Lévy index of stability, β is the asymmetry parameter, γ is the so called scale parameter, and δ is the location parameter of this distribution, while n determines the convention [27,28,31]. In this paper, we follow the conventions defined in Ref. [31], that correspond to n = 2. In this case, for α = 1, the Fourier-transformed Lévy stable source distribution reads asS By now, the Lévy stable distributions are implemented into commercially available software packages, for example, Mathematica [27]. In practice, we utilize the parameterization that is continuous in the Lévy index of stability α L [27]. Given that the Gaussians correspond to Lévy-stable source distributions with α L = 2 (the value of the exponent in the Fourier-transformed Gaussians) and taking into account, that in our analysis the Gaussian elastic amplitude t el (b) has the exponent α = 1, we conclude that the Lévy index of stability α L is simply twice the exponent of our Lévy series, i.e.
Recently, the TOTEM Collaboration published a high precision measurement of the low-|t| region of the differential cross-section of elastic pp scattering at √ s = 8 TeV [25]. This demonstrated a significant, more than 7σ deviation from a simple exponential cone behaviour, corresponding to a Gaussian representation of the elastic scattering amplitude. In our language, this means that α < 1, or, using the standard form of the Lévy index of stability, α L = 2α < 2 for this data set.
Subsequently, let us present the results of our data analysis and indicate the best Lévy expanison fits to elastic pp and pp differential cross-sections. Let us proceed to evaluate the shadow profiles P (b, s) and the slope parameters B(s, t) as well as the nuclear phases φ(s, t) for the available values of s and for both proton-proton and proton-antiproton elastic scattering reactions, to find their excitation functions, and to compare proton-proton and proton-antiproton results, as described in the next section.

III. DATA ANALYSIS
Let us test the power of our Lévy expansion method first on the already published differential cross-section data of elastic pp collisions at √ s = 7 TeV, utilizing the published TOTEM data set of Ref. [32]. Fig. 2 indicates that the 7 TeV TOTEM data set can be represented by a fourth-order Lévy expansion with a reasonable χ 2 /NDF = 224/154, that corresponds to a marginal confidence level of CL ≈ 0.02%. Inspecting Fig. 2 by eye suggests also that the parameters of the Lévy expansion in Eq. (2.1) can be interpreted as they closely represent the data. These parameters are printed on the right-hand side of the top panel of Let us also investigate in detail the recently released, new 13 TeV TOTEM preliminary data set [19], to look for crossing (C) odd effects in the comparison of elastic pp and pp collisions. In what follows, we consider four different aspects of the TOTEM data in comparison with elastic scattering data at lower energies, both for pp and pp collisions. Namely, we compare the shadow profile functions, the t-dependence of the elastic slope parameter B, the same for the ρ-parameter and the so-called nuclear phase φ(t), that measures the argument of the elastic scattering amplitude. Finally, we show in a simple and straightforward analysis of a large-|t| region beyond the diffractive minimum and maximum, that the differential cross-section of elastic pp scattering evidences a proton substructure of two distinct sizes for GeV and TeV energy ranges, respectively.

A. Looking for Odderon effects
As noted in Refs. [4,33], the only direct way to see the Odderon is by comparing the particle and antiparticle scattering at sufficiently high energies provided that the high-energy pp and pp elastic scattering amplitude is a difference or a sum of even and odd C-parity contributions. The even-under-crossing part consists of the Pomeron and the f Reggeon trajectory, while the odd-under-crossing part contains the Odderon and a contribution from the ω Reggeon, i.e.
It is clear from the above formulae that the odd component of the amplitude can be extracted from the difference of the pp and the pp scattering amplitudes. At sufficiently high energies, the relative contributions from secondary Regge trajectories is suppressed, as they decay as negative powers of the colliding energy √ s. The vanishing nature of these Reggeon contributions offers a direct way of extracting the Odderon as well as the Pomeron contributions, T O el (s, t) and T P el (s, t), respectively, from the elastic scattering data at sufficiently high colliding energies.
In Refs. [4], the authors argued that the LHC energy scale is already sufficiently large to suppress the Reggeon contributions, and they presented the (s, t)-dependent contributions of an Odderon exchange to the differential and total cross-sections at LHC energies. That analysis, however, relied on a model-dependent, phenomenological extension of the Phillips-Barger model [34] and focussed on fitting the dip region of elastic pp scattering, but it did not analyse in detail the tail and cone regions. In fact, that analysis relied heavily on the extrapolation of fitted model parameters of pp and pp reactions to exactly the same energies. Similarly, Ref. [13] also argued that the currently highest LHC energy of √ s = 13 TeV is sufficiently high to see the Odderon contribution, given that the Pomeron and the Odderon contributions can be extracted from the elastic scattering amplitudes at sufficiently high energies as One of the problems is that the elastic pp and pp scattering data have not been measured at the same energies in the TeV region so far. So, we strongly emphasize the need to run the LHC accelerator at the highest Tevatron energies of 1.96 TeV, in order to make such direct comparisons possible. Another problem is a lack of precision data at the low-and high-|t|, primarily, in pp collisions. Nevertheless, we show that robust features of the already performed measurements provide not only an Odderon signal, but they also indicate the existence of a proton substructure.
In this paper, we take the data as given and do not attempt to extrapolate the model parameters for their unmeasured values. Instead, we look for even-under-crossing and oddunder-crossing contributions by comparing pp and pp collisions at different energies, looking for robust features that can be extracted in a model-independent manner. In addition, we build upon Ref. [4] by assuming, as justified by that analysis, that the Reggeon contributions to the elastic scattering amplitudes are negligible if √ s ≥ 1.96 TeV. Let us first of all compare the behaviour of the shadow profile functions P (b) before investigating the four-momentum transfer dependent B(t) functions for both pp and pp reactions.

B. Excitation function of the shadow profiles
The excitation of the shadow profiles is obtained from the elastic scattering amplitude obtained by fits to pp elastic scattering cross-section data from √ s = 23.5 GeV to 13 TeV, as illustrated in Fig. 8. The excitation function of the shadow profile functions for pp reactions is indicated in Fig. 9. In pp collisions at the lower ISR energies, the shadow profile functions look nearly Gaussian, and their values at zero impact parameter are below unity, P (b = 0) < 1. The picture changes at the LHC energies of 7 and 13 TeV, where the shadow profile functions seem to saturate, with P (b = 0) > 99.9% in an extended range, for b < 0.4 fm at 7 TeV and b < 0.5 fm at 13 TeV. This indicates that the black disc limit is reached in the center of these collisions, corresponding to P (b) ≈ 1. However, outside the 0.4 or 0.5 fm saturated regions, the P (b) decreases nearly in the same manner, as at lower energies. One may conclude that a new, black region opens up in the TeV energy region, which increases with growing colliding energies, and it is surrounded by a gray hair or skin region, that has a "skin-width" that is approximately independent of the energy of the colliding protons. Thus, with an increase of colliding energies, the protons become blacker, they do not become edgier but become larger. This is the so called BnEL effect [35], which can be contrasted to the earlier expectations, the so-called BEL effect suggesting that with increasing energy of the collisions, the protons might become blacker, edgier and larger. A new trend opens up with the 7 TeV TOTEM data, that indicates a black region, P (b) 1 up to a radius of about 0.4 fm and the size of this black region is increasing with an increase of colliding energies. Note also that at the ISR energy range, √ s ≤ 62.5 GeV, the shadow profiles are very similar, however, at the TeV energy range, pp and pp collisions evolve somewhat differently. For example, in the shadow-profile of the elastic pp collisions at √ s = 1.96 TeV, the nearly flat region with P (b) ≈ 1 is not yet present, while this region is present and it is rather extended in the shadow profiles of elastic pp collisions at √ s = 7 and 13 TeV.
In both Figs. 8 and 9, one can observe that the proton becomes blacker and larger with increasing energies, however, its edge is apparently nearly constant, looks like a nuclear skin that has the same skin-width regardless of the energy of the collision. These results are similar to earlier observations, published in Refs. [35][36][37][38]. To highlight this point, we have plotted the shadow profile functions P (b|pp, 13 TeV), P (b|pp, 7 TeV) and P (b|pp, 1.96 TeV) together in Fig. 10 containing the effects that come from evolution of the structure of the proton with increasing √ s. We find that the energy evolution of the shadow profiles is similar for the C-even pp collisions and for the C-odd pp collisions.
We summarize that from the shadow profile functions it is very difficult to draw strong conclusions given that the model-independent method does not allow estimation of the collision energy dependence of the model parameters yet, so it is very hard to tell if the obvious difference between the pp and the pp collisions is due to the difference in the energy of the collisions or not.
This observation underlines the importance of data-taking at the LHC, a pp collider, with √ s decreased to as close as reasonably possible to 1.96 TeV or 1.8 TeV, the energy range of the pp collisions measured at Tevatron.

C. Results for the nuclear slope parameter B(t)
Let us move on to the analysis of another important characteristics, namely, the tdependent nuclear slope B(t).
In pp collisions, the analysis of the four-momentum transfer and center-of-mass energy dependent nuclear slope, B(t, s) is summarized in Fig. 11. Surprisingly, in the low-|t| region, where a diffraction cone is expected, we find that B(t) is actually not exactly constant, but a t-dependent function, so the exponential behaviour can only be considered as an approximation, as clearly shown in Fig. 11. In the ISR energy range 23.5 GeV ≤ √ s ≤ 62.5 GeV, the nuclear slope can be considered as roughly constant both in the |t| ≤ 1.0 GeV 2 (diffractive cone) and in the 2.0 ≤ |t| ≤ 3.0 GeV 2 (tail) region, with nearly ISR energy independent B cone (pp|ISR) ≈ 10 GeV −2 and B tail (pp|ISR) ≈ 2 GeV −2 , rather surprisingly. At the LHC energy scales √ s =7 and 13 TeV, the cone region shrinks, as expected, down to |t| ≤ 0.3 GeV 2 , however, rather unexpectedly and surprisingly, the tail region opens up with a featureless, nearly flat B(t) function that extends the tail region to 1.0 ≤ |t| ≤ 3.0 GeV 2 . An obvious and numerically very stable observation is that the low-|t| approximate value, B cone (pp|LHC) ≈ 20 GeV −2 is nearly a factor of two larger than the corresponding values at ISR, but also it is clear that these values are significantly t-dependent. On the other hand, in the large-|t| region, B tail (pp|LHC) ≈ 5 GeV −2 , valid in a broad, 2 GeV −2 wide range of |t|, with not larger than 20 % level variations over this range. This suggests the existence of some substructures inside the protons that we detail in the next subsubsection.

Substructures of protons from B(t) at large t
It is important to realize, that the asymptotic value B tail (pp|LHC) ≈ 5 GeV 2 is nearly independent of the 7 or 13 TeV colliding energies and it is apparently significantly larger than the same asymptotic values at the ISR energies. It is rather clear that these values of B(t), nearly constant over such rather extended 1 or 2 GeV 2 ranges (plateaux) of |t|, indicate two different and nearly Gaussian shaped proton substructures, a smaller substructure at ISR energies, that corresponds to B tail (pp|ISR) ≈ 2 GeV −2 , and a larger substructure at LHC energies, that corresponds to B tail (pp|LHC) ≈ 5 GeV −2 . It is remarkable that the size of these substructures is not changing when the center of mass energy of the collision is varied in the √ s = 23.5 to 62.5 GeV ISR range, or in the LHC energy range of √ s = 7 to 13 TeV. It is desirable to take more data and to investigate what happens in between these energy ranges. In particular, we suggest that the colliding energy of the LHC accelerator be varied in the broadest possible range, from √ s = 900 GeV to the designed top colliding energy of 14 TeV, to see particularly if a smooth or sudden transition is seen in B tail (pp) between 2.76 TeV and 7 TeV, or not.
In pp elastic collisions, a similar analysis of the four-momentum transfer squared and center-of-mass energy dependent nuclear slope, B(t, s) is summarized in Fig. 12. In this analysis, the data are less detailed and only one data set is analyzed in the ISR energy range, corresponding to √ s = 53 GeV. This data set has points both in the cone and in the tail regions, while the dataset at √ s = 546 GeV is detailed in the low-|t| region but lacks data in the tail. The data at √ s = 630 GeV lacks data in the cone region, but extends more to the tail range, finally the data that we analyze at the Tevatron energy scale, at √ s = 1.96 TeV, have limited |t|-range that only partially covers the cone and the tail regions.
For a kind of uniformity of the comparisons of pp elastic scattering data at various energies, we thus rely on fits and extrapolations with fixed α = 0.9, as detailed in Appendix B, and also summarized in Fig. 5. In the low-|t| region, where a diffractive cone is expected, we find that B(t) is actually not exactly constant, but a t-dependent function, so the exponential behaviour can only be considered as an approximation at all the considered energies, as clearly seen in Fig. 12. In pp collisions, the nuclear slope is approximately a constant both in the |t| ≤ 0.5 GeV 2 (cone), as well as in the 2.0 ≤ |t| ≤ 3.0 GeV 2 (tail) region. These data sets cover a broad energy range, and Fig. 12 clearly indicates, that the approximate average values of B cone (pp) increase monotonically with an increase of √ s. It is remarkable, that the slope parameter B(t) is, rather surprisingly, tending to an energy-scale independent asymptotic value of B tail (pp) ≈ 5 GeV −2 . The range, over which the asymptotic exponential region prevails, is apparently extending at least up to 3 GeV 2 or more. The larger the colliding energy, the broader this region, which starts to open at |t| 1.5 GeV 2 at √ s = 1.96 TeV. The asymptotic value B tail (pp|1.96 TeV) ≈ 5 GeV 2 is nearly the same as the value of B tail (pp|LHC) at √ s = 7 or 13 TeV colliding energies and is thus almost independent of the type of the collisions as well.
It is clear that these two different, but nearly constant asymptotic values of B(t), corresponding to B tail (pp|1.96 TeV) ≈ B tail (pp|LHC) ≈ 5 GeV −2 and B tail (pp|ISR) ≈ 2 GeV −2 over extended, 1 or 2 GeV 2 wide four-momentum transfer squared ranges exhibit a domain with a nearly exp(−B|t|) behaviour. Thus, this domain reveals the existence of a proton substructure with a nearly Gaussian elastic scattering amplitude distribution, t el (b) ∝ exp(−b 2 /(2R 2 )). As is well known, an approximate value of the slope parameter B tail is proportional to the squared Gaussian radius of such a substructure. The larger radius of this substructure is observed in the TeV energy range, both in pp and in pp collisions, while a smaller-size substructure is seen in pp collisions in the √ s = 23.5 − 62.5 GeV ISR energy range. From the relation R 2 = 4B (in natural units) [24,39], the Gaussian radius of the substructure is about R LHC ≈ 0.9 fm (TeV energies) and R ISR ≈ 0.6 fm (few 10 GeV energies). The analysis of the proton substructure and the determination of its contribution to the total cross-section is detailed in the next subsection, as well as in Appendix C.
The values R LHC and R ISR are strikingly similar to the radii of an effective diquark (R d ) and quark (R q ), respectively, that were independently obtained to characterize the substructures inside the protons in Ref. [35]. It turned out that the unitarized quarkdiquark model of elastic pp scattering (called the Real extended Bialas-Bzdak or ReBB model) predicted the √ s dependence of the total cross-section, the dip position and even certain scaling properties of the differential cross-section of elastic pp scattering at √ s = 13 TeV with a reasonably good accuracy, based on its tuning at ISR energies and the TOTEM data set at √ s = 7 TeV. So, the hypothesis about a proton substructure gains a larger weight and evidence in our analysis and, thus, definitely deserves more detailed investigationsthat, however, go beyond the scope of the model-independent approach elaborated in this work. In particular, an important question is whether the observed two distinct scales R LHC and R ISR correspond to the dressed diquark and the dressed quark, respectively, or simply represent a single substructure whose size grows with energy, remains open.
A dynamical model for the elastic amplitude based upon a two-scale structure of the proton was previously proposed also in Refs. [40][41][42]. In this model, while the first scale was associated with the confinement radius R c 1 fm and can be attributed to the proton "shell", the second semi-hard scale r 0 ≈ 0.3 fm originates due to non-perturbative inter- actions of gluons and characterizes an effective gluonic "spot", or a cloud, around each of the valence quarks. Despite somewhat different values of the physical scales adopted in this model, it has appeared to predict the energy dependence of the total and elastic cross sections quite accurately, at least, in a parameter-dependent way. In our current work, however, instead of reviewing various possible model interpretations, we study the model-independent properties of the elastic scattering data and search for C-odd (or Odderon) effects. We employ our model-independent imaging method to sharpen the picture of the proton as can be "seen" by elastic scattering measurements at different energies.

Odderon effect and the difference between B(t|pp) and B(t|pp)
We can observe in Fig. 11, that in pp collisions B(t) starts in the cone region with nearly constant values. However, after the cone region, B(t) increasea shortly then falls sharply, to cut through the B(t) = 0 line and to reach a deep minimum with B(t) 0 values. Then this function B(t) starts to increase, it cuts through the B(t) = 0 line second time, from below, to approach slowly its asymptotic value, B tail .
Remarkably, for pp collisions, the B(t) functions behave apparently in a qualitatively different way. In this case, after the cone region, B(t) approaches zero and it may marginally cross zero, but not so deeply and sharply as in the case of pp collisions. Taking into account that the published error on B in pp collisions is about 0.5 GeV −2 [43], the error on the extrapolated B(t) function can be estimated also. The latter appears to be similar or even larger as compared to the error of B at the optical point t = 0. So it seems to us, that the crossing of the B(t) function below zero is within errors and thus it is likely not a significant effect in any of the pp collision data.
To clarify this point more, we compare the B(t) functions for pp collisions at √ s = 7 and 13 TeV with that of the pp collisions at √ s = 1.96 TeV, see Fig. 13. The nuclear slope in pp collisions becomes clearly and significantly negative in an extended |t| region, starting from the diffractive minimum (dip) and lasting to the subsequent diffractive maximum (bump). These dip and bump structures are clearly visible in the corresponding data sets, as visualized in Figs. 2 and 3 as well. In contrast, for pp collisions in the Tevatron energy range, we do not find any dip and bump structure, that would correspond to a |t|-region where B(t) were negative.
We have performed further tests to cross-check if, within the errors of the analysis, the dip and bump structure is indeed absent in pp collisions at √ s = 1.96 TeV, or not. First of all, one can directly inspect Appendix B to see that any reasonably smooth (e.g. spline) extrapolation of the pp data points would lack a diffractive minimum structure at √ s = 1.96 TeV. One may argue that we do not see the minimum because the values of α and R were fixed. However, these numbers specify the approximate Lévy shape only, dσ/dt ∝ exp (−(R 2 |t|) α ), that decreases monotonically, so their variation cannot cause diffractive minima or maxima, as also apparent on the dashed lines of partial fits described in Appendices C and D. In any case, we have also tested numerically that changing α in the region of fixed 0.8 -1.0 does not qualitatively change the behaviour of B(t). Within the allowed range of variation of the essential Levy expansion parameters c i , we find that the diffractive minimum is lacking in elastic pp collisions at 1.96 TeV.
This lack of diffractive minimum as well as the lack of the subsequent diffractive maximum in elastic pp collisions is contrasted to the strong diffractive minimum and maximum (the dip and bump structure) in elastic pp collisions at all investigated energies, hence it indicates a rather evident C-odd contribution to the elastic scattering amplitude, the so called Odderon effect, shown also in Fig. 13.

D. Evidence for proton substructure from Lévy fits at large |t|
Lévy fits to the tails of the differential cross-section of pp elastic scattering data from √ s = 23.5 GeV to 13 TeV are shown on a summary plot in Fig. 6 and detailed in Appendix C. Note, that all the tails of the differential cross-sections are nearly linear on a log-linear plot, indicating a nearly Gaussian substructure of the proton. As was discussed earlier, two distinct sizes of such a substructure are seen, given the two different values of the slope in the ISR energy range of √ s = 23.5 -62.5 GeV, and in the LHC energy range of √ s = 7 -13 TeV. After the dip-bump structure, the differential cross-section of elastic pp collisions can be described by a simple A exp (−(|t|R 2 ) α ) form, with α = 0.9 ± 0.1 value. Thus, for illustration this plot was done for fixed value of α = 0.9.
We found that in the 23.5 ≤ √ s ≤ 62.5 GeV range, a proton substructure of nearly constant size (within errors) was present, with a characteristic Lévy length scale of R ISR = 0.3±0.1 fm, and with the corresponding contribution to the total cross-section σ ISR = 0.3 +0.3 −0.1 mb, where the quoted errors take into account also the errors coming from the variation of the value of α between 0.8 and 1.0, see also the ISR plots of Appendix C for details. Fig. 14 indicates the TOTEM preliminary elastic scattering data at √ s = 13 TeV with their fourth-order Lévy expansion fits. A power-law tail would show up as a straight line on this plot, but apparently it does not yet show up on the currently available |t|-range that extends up to about t max = 4 GeV 2 . Although a straight line fit to the tail of this distribution is perhaps possible starting from |t| ≥ 2 GeV 2 , these points are getting close to the end of the TOTEM acceptance for this data set, and the error bars are getting large. To clarify the existence of such a possible power-law tail, more data at larger values of |t| would be desirable. In contrast, the data in the well measurable |t|-range are sufficient to demonstrate the nearly exponential behaviour of the differential cross-section in the tail region that follows the dip and bump structure. The existence of the nearly Gaussian substructure is thus an obvious and rather robust feature of the data (for more details, see Appendix C). This claim is also supported by the nearly energy-and t-independent values of the slope-parameter B(t) in the tail regions, see Figs. 11 and 12.

E. Results for ρ(t)
By reconstructing the elastic scattering amplitude from the data, we have also found the t-dependent ratio of its real to the imaginary parts in pp collisions, the ρ-parameter. Such a result is illustrated in Fig. 15 and indicates that the ρ-parameter is significantly t-dependent. This dependence is initially nearly linear in t.
However, at the 7 and 13 TeV LHC energies, ρ(t) starts to diverge to minus infinity which corresponds to a zero point or node of the imaginary part of the elastic scattering amplitude. In order to illustrate this point, we show the real and imaginary parts of the elastic amplitude in Figs. 16 and 17.
Let us now analyze in detail the nuclear phase φ(t), the argument of the complex elastic scattering amplitude T el , that is traditionally measured in units of π, as described in the next subsection. In this subsection, let us investigate in detail if we can identify the Odderon effects in the t-dependence of the nuclear phase φ(t). As we clarify in Appendices A and B, such a phase can be reconstructed mostly from the differential cross-section at low-momentum transfers squared. In Appendix B we have demonstrated that ρ(t = 0) cannot be reliably extracted from the studied elastic pp collision data, in particular, due to a significant lack of the low-|t| data points. However, in Appendix A we demonstrated that ρ(t) can be extracted from Lévy  fits to elastic pp scattering, with the exception of √ s = 44.7 GeV and 52.8 GeV data, where the confidence level of our fits is not in the statistically acceptable domain.
We have evaluated the t-dependent nuclear phase for pp collisions in the both ISR and  LHC energy ranges. Unfortunately, the ISR range is a cumbersome one, where the Reggeon as well as the Pomeron and Odderon contributions mix with each other. In this energy range, for some of the pp collision data sets, we found that φ(t) reaches π value at the same (or close) |t| values above 1 GeV 2 while our sensitivity studies indicate that in this region our method to reconstruct the nuclear phase has increasing systematic difficulties. So, we cannot at present reliably test if the points where φ(t) = π coincide or not at low energies in pp and pp collisions. We can however make a statement that the pp collisions at the ISR energy range of 23.5 ≤ √ s ≤ 62.5 GeV, the nuclear phase does not reach π for |t| ≤ 1 GeV 2 . Fortunately, a similar analysis in the TeV energy range gave interesting results. Fig. 18  Here, for illustration we also kept the corresponding principal values for the nuclear phase satisfying 0 < φ 2 (t) < π, also shown in units of π. The latter are indicated by thinner and discontinuous curves explicitly seen at large |t| > |t dip |.
indicates that the phase φ(t) reaches the π value at 7 and 13 TeV simultaneously at |t 0 | (pp) ≈ 0.45 ± 0.05 GeV 2 in pp collisions, while φ = π is reached at a rather different value of |t 0 | (pp) ≈ 0.70 ± 0.05 for 1.96 TeV pp collisions. This is an important qualitative feature of φ(t), that indicates a significant Odderon contribution, that apparently cannot be attributed to an s-dependent effect. This subtle Odderon signature is clearly illustrated in Fig. 18 and cannot be directly seen in the differential cross sections. Thus, the √ s independence of the crossing point of the nuclear phase φ(t) = π in the TeV energy range for pp collisions and its strong shift in pp collisions in the TeV energy range, where the Reggeon contributions are apparently negligible, can be considered as a second reliable signature of the Odderon.

IV. DISCUSSION
Our experience is that many of our results, for example, for ρ(t) and also for the nuclear phase φ(t) in the |t| > 1 GeV 2 region are very sensitive to the precise details of the fits, so interpretation of the TOTEM results at √ s = 13 TeV depends critically and sensitively on the currently preliminary TOTEM data points and their error bars. Certain features of our analysis, for example, the behaviour of B(t) in the range where the slopes could directly be evaluated from the data are more robust and stable.
Thus, we feel strongly motivated to warn the astute reader against the over-interpretation of model results that indicate certain features of the elastic scattering data correctly only on the qualitative level, but fail miserably on a confidence level test. Actually, the Odderon effects that we discuss in detail in this work are due to some robust and model independent features of the data, but we have investigated other more subtle effects too that we do not emphasize in this work.
Our key point is that the significance of the new Odderon effects can be revealed only if the data sets are final, with published statistical and systematic errors, and if they can be correctly and faithfully represented by theoretical calculations. So we recommend to carefully evaluate the confidence levels of all subsequent theoretical analysis of final TOTEM data in future analyses and to determine the significance of the presence of novel effects like the Odderon contributions. Our method, presented in this work, allows for such sensitivity and significance analysis, however, given that the TOTEM data at √ s = 13 TeV are preliminary, the evaluation of the systematic errors of the our fit results would most likely also be premature at present.
Nevertheless, we may warn the careful readers that descriptions of possible Odderon effects or the lack of them, based on data analysis with zero confidence levels might have apparently been over-interpreted recently: the significance of the interpretation of fits that do not describe the data in a statistically acceptable manner is not particularly well defined. Based on our experience, we recommend against the over-interpretation of the data in terms of models that do not have a confidence level of at least CL ≥ 10 −3 %, but in a final analysis, we strongly recommend to rely on descriptions that have a confidence level of at least CL > 0.1% before the model results can be interpreted. The authors should also check if their optimization procedure has converged or not and test if the error matrix is accurate and the estimated distance to the real minimum is sufficiently small.
In our studies, we have found that small variations in the fit range or in the values of the fit parameters do not change the following robust features of the data, that we highlight below.

A. Robust qualitative features
While searching for the differences between the differential cross-sections of elastic pp and pp collisions, that could exhibit a C-odd contribution, or an Odderon effect in 13 TeV pp elastic scattering, we have evaluated the shadow profile functions at 7 and 13 TeV pp collisions using a novel imaging method, the model-independent Lévy expansion.
We have compared the shadow profile functions of pp and pp collisions at various energies. We have found that the shadow profiles saturate at the LHC energies of 7-13 TeV: for small values of the impact parameter, a P (b) 1 region opens up. With increasing the collision energies, the protons become blacker, but not edgier, and larger, confirming the BnEL effect, that was reported in Ref. [35].
We see a significant difference between the shadow profile functions P (b) of protons and anti-protons in the TeV region, but from the current analysis we cannot determine uniquely, if this difference is an Odderon effect, or, an effect of saturation that is apparent also in pp collisions with an increase of collision energy. We would need pp and pp elastic scattering data at exactly the same collision energies, that can be realized these days only by running the LHC accelerator at energies close to the Tevatron energy scales of √ s = 1.8 -1.96 TeV. In Subsection II B, we have analyzed the dependence of the B-slope parameter on the four-momentum transfer squared t in pp as well as in pp reactions. We have found, that the B(t) functions indicate an Odderon effect very clearly. Surprisingly, we have identified a |t| region after the dip and the bump structure, where a clear-cut evidence is seen in our analysis for a proton substructure of two distinct sizes in two experimentally probed GeV and TeV energy ranges. In every case, such a substruc-ture is characterized by a Lévy exponent of α = 0.9 ± 0.1. Possible evidence for a new substructure inside the proton at LHC energies was, as far as we know, first pointed out by Dremin in Ref. [16]. In this paper, we explored this possibility in detail and identified the characteristic Lévy exponent as α = 0.9 ± 0.1, characterized the substructure with an approximate Gaussian and Lévy scales at the ISR and LHC energies, respectively, and determined the corresponding contribution to the total pp cross-section as well, as described in Subsection II B and in Appendix C. Based on the Gaussian sizes found in Subsection II B, it is tempting to note that they are strikingly similar to a dressed quark and a dressed diquark, found to describe elastic pp scattering in an earlier, model dependent analysis [35]. Their presence seems to provide a phenomenological support for the quark-diquark picture of the proton, which is deeply related to the solution of the confinement problem in QCD proposed recently by Brodsky and collaborators in Ref. [44]. This, however, does not exclude the possibility for a single substructure growing with energy in such a way that the substructure grows only in between the ISR and the LHC energies, but remains constant between √ s = 23.5 and 62.5 GeV, then it grows but again remains of constant size between √ s = 7 and 13 TeV. More measurements would be strongly desirable to justify the quark-diquark picture or to investigate the growth of a single substructure (a dressed quark) with increased colliding energy. Presently this is feasible only by running the LHC accelerator at decreased energies, varying the collision energies between √ s = 900 GeV and 7 TeV, to identify the transition energy. Another measurement at the desinged top LHC energy of 14 TeV is already planned and approved, as far as we know, by the TOTEM Collaboration at CERN LHC.

B. Highlighted results
Let us highlight some of the important points of our study: 1. We have found a solid, stable and clear-cut evidence for a proton substructure, with two different sizes extracted for two distinct (GeV and TeV) energy ranges that are similar to the sizes of a dressed quark and a dressed diquark, respectively, as discussed in Ref. [35], and as also derived from QCD in Ref. [44]. Fig. 13, even without the quantitative results, demonstrates the existence of such a proton substructure, corresponding to the second extended plateau in the large-|t| region (besides the usual elastic cone region at low |t|) with a nearly exponential contribution to the differential cross-section of elastic pp scattering. We noticed that such a plateau corresponds to a dressed quark-scale substructure in the lower √ s = 23.5 -62.5 GeV energy range, while it resembles a larger, dressed diquark-scale substructure in the 7 -13 TeV energy range, see Fig. 6. This, however, does not exclude the possibility for a single substructure growing with energy in such a way that the substructure grows only stepwise: within the resolution errors, it remains of a constant size between the √ s = 23.5 and 62.5 GeV ISR energies, then it grows but its size remains constant within the experimental resolution error of about 0.1 fm in the energy region between √ s = 7 and 13 TeV. In short we observed two significantly different in size substructures inside the protons at ISR and LHC energies but the physical interpretation of this substructures is an open question, and requires more experimental and theoretical investigations.
2. At each energy and for each investigated data set, the scattering amplitude of elastic pp and pp collisions was described by our new Lévy series expansion method. With the help of the elastic scattering amplitude, we have reconstructed values for the total cross-section and for the differential cross-section of elastic scattering. The published values of the total cross-sections were reproduced within errors and the fits to the differential cross-section looked fine. In case of several data sets they have also passed the more stringent tests of mathematical statistics, namely the confidence level of most of our fits was not unacceptable from the point of view of mathematical statistics, either, with confidence levels CL > 0.1 %.
3. For all those data-sets, where the confidence level of the fit was not unacceptable from the point of view of mathematical statistics, we found that the exponent α was significantly less than unity in pp collisions, the deviation being a more than 5σ effect. Given that exponential behaviour corresponds to the case of α = 1, and fits with CL > 0.1 % represent the data undoubtedly, we find that the differential cross-section at low |t| is apparently non-exponential in 23.5, 30.7, 62.5 GeV and in 13 TeV pp elastic scattering data. It is quite remarkable, that the corresponding values of the non-exponentiality are α = 0.88 ± 0.01, 0.89 ± 0.02, 0.90 ± 0.01, 0.90 ± 0.01, when rounded up to two decimal digits. This implies that an average value of α = 0.89±0.02 is consistent with all the measurements in a very broad energy range from 23.5 GeV to 13 TeV. The energy independence of this α = 0.89 ± 0.02 value calls for a physics interpretation, and for further studies.
4. Signals of non-exponentiality in the cone region are also indicated in pp elastic scattering data at all energies, where we have been able to describe all the analyzed datasets with an α = 0.9 fixed value, from √ s = 53 GeV to 1.96 TeV. At a first sight, difference of α from unity may imply non-linearity of the Regge trajectories. However, it is known that the extrapolation of the Regge trajectories from masses to negative values of momentum transfer t depends on the assumed analyticity of the elastic scattering amplitude. It is worth to mention here that our amplitudes are singular at t = 0 for α < 1, and this behaviour does not allow for an easy analytic continuation and, hence, an interpretation of non-exponentiality is not straightforward in terms of Regge theory.
5. This non-exponential nature of the differential cross-section in the low-|t| region implied that the slope parameter B = B(t) is a function that is strongly dependent on the |t| range, where it is determined. We have evaluated B(t) both numerically and analytically in the whole |t| region, where it is defined. The extrapolation of B(t) to the optical point of t = 0 turned out to be model-dependent, not only due to the fact that there is a Coulomb effect that induces CNI terms and modifies the slope at very small |t|, but also due to the analytic result that lim t→0 B(t) = ∞, i.e. our Lévy expansion method is non-analytic at the optical point. This closely corresponds to the physical picture that we allow for the underlying source contributions that may have infinite root mean square, which is typical for a Lévy stable source distribution. (For a typical example, one may consider a Lorentzian distribution, that is a symmetric, Lévy stable distribution with Lévy index of stability α L = 1, corresponding to our non-exponential parameter α = 0.5. 6. The failure of the leading-order Levy fits for fixed α = 0.9 at √ s = 7 and 13 TeV, while their success at any lower energy, have indicated that not only the size of the proton increases with an increase of collision energies, but also the shape of the protons changes. These results can be explained due to an emergence of the saturated plateaux P (b) ≈ 1 in the small b ≤ 0.4 − 0.5 fm region at √ s = 7 and 13 TeV which were not seen at lower energies. In our analysis, the data points from the CNI region, i.e. at |t| < 0.01 GeV 2 , were not included neglecting any possible Coulomb effect in the Levy expansion.
7. The lack of detailed data in the very low-or very large-|t| regions in pp collisions prevented us to determine precisely the ρ and the B parameters of this case.
8. The analysis of the four-momentum transfer squared and the center of mass energy dependent nuclear slope, B(t, s) in Fig. 11 not only confirms the existence of a proton substructure (corresponding to the existence of large regions in t where B(t) is approximately but not exactly constant), but also indicates a sharp difference between the B(t) functions of pp and pp collisions, when comparing Figs. 11 and 12. This is a clear-cut and significant Odderon effect. For pp collisions, B(t) starts with positive values, then it cuts sharply through the B(t) = 0 line and returns above it shortly but very significantly, corresponding to the dip and bump structure in the differential cross-section for elastic scattering in each of the investigated data sets. In contrast, in pp collisions, B(t) approaches zero but within the errors of the analysis it does not cross it at √ s = 1.96 TeV. At two of the lower energy scales of √ s = 630 and 53 GeV, B(t) apparently crosses zero and develops minima whose values are 2 and 1 GeV −2 , respectively. However, taking into account that the published error on B in pp collisions is about 0.5 GeV −2 [43], the error on the extrapolated B(t) function can be estimated to be similar or larger as compared to the error of B at the optical point of t = 0. So it seems to us, that the crossing of the B(t) function below zero is within errors likely not a significant effect in any of the pp collision data. In the only data set that we could access in the TeV energy range, where the complicated Reggeon contributions are already negligible, B(t) does not cross zero, so the diffractive minimum and maximum, the dip-bump structure is lacking in these √ s = 1.96 TeV pp elastic collisions. Such a behaviour is in sharp contrast to the pp differential cross-sections at all energies. Apparently, this is a clear-cut Odderon effect, as illustrated in Fig. 13. 9. In addition, we have also found a surprisingly clear Odderon effects in the t-dependent nuclear phase φ(t). Our analysis indicates that this phase reaches π value at the very different, high LHC energies of 7 and 13 TeV in elastic pp collisions at the same value of |t 0 | ≈ 0.45 GeV 2 , that suggests that this |t 0 | value has a negligibly small dependence on the collision energy, √ s. However, in elastic pp collisions at √ s = 1.96 TeV, such a crossing point in the nuclear phase φ(t) = π is at a very different location from pp collisions. 10. We have also found a weak difference between the shadow profile functions P (b) belonging to pp collisions at √ s = 1.96 TeV as compared to that of pp collisions at 7 and 13 TeV. However, we also found a significant evolution of the shadow profile functions from 23.5 to 62.5 GeV and from 7 TeV to 13 TeV in elastic pp collisions. 11. In order to clarify if the difference between the shadow profiles of pp and pp collisions occurs due to the change of the colliding system type or due to the change of t he center of mass energy of the collisions, as well as to clarify differences of the other observables like B(t) and φ(t), we strongly recommend to run the LHC measurements at lower energies, preferably as close to √ s = 1.96 TeV, as reasonably achievable, to measure the difference between pp and pp collisions exactly and to clarify the Odderon contribution without any possible energy evolution and extrapolation effects.
12. We recommend extreme care before drawing big conclusions, given that we see the sensitivity of some of the details like φ(t) at large |t| for tiny details in the data and in changing some of the higher order coefficients of the fits. Possibly these tiny details differ in some papers that may apparently draw big, but contradicting, and not particularly well founded conclusions about the existence or non-existence of the Odderon effects. When looking for a robust conclusion about the Odderon contribution, we recommend to look at the summary plot of the Lévy fits to the slope of the differential cross-sections of elastic scattering, as indicated in Fig. 13.

V. SUMMARY AND CONCLUSIONS
In summary, we conclude that we have found clear-cut and post-factum rather obvious differences between the differential cross-sections of elastic pp and pp collisions, indicating a C-odd contribution: the Odderon effect. This corresponds to a small difference of the tdependent slope parameters between pp collisions at 13 TeV and 7 TeV collision energies at the LHC as compared to a large change of the t-dependence of the slope parameter B(t) in 1.96 TeV pp collisions. Another characteristic Odderon signal is the difference is between the existence of a diffractive minimum and maximum in both 13 and 7 TeV elastic pp scattering, corresponding to two distinct crossing-points of the B(t) functions with the B(t) = 0 line, as contrasted to the monotonically decreasing differential cross-section of elastic pp collisions, with a t-dependent elastic slope that is B(t) > 0 significantly, a function that never crosses the B(t) = 0 line. These Odderon signals, the change in B(t) and the disappearance of the diffractive minimum as well as the diffractive maximum, when changing from pp reactions to pp reactions, are rather obvious, stable and clear-cut effects. Once they are identified, they can be directly seen on the data sets, when one plots the 13 TeV and 7 TeV differential cross-section of elastic pp scattering on the same plot with the differential cross-section of elastic pp scattering, as illustrated in Fig. 19.
We have confirmed such Odderon effects with a more refined and subtle analysis, that indicated the lack of energy dependence of the crossing point |t 0 | of the nuclear phase φ(t) both in the TeV region φ(t) = π at the same value of |t 0 |(pp) ≈ 0.45 ± 0.05 GeV 2 both at 7 and 13 TeV. When evaluating the crossing point of the nuclear phase for pp collisions at √ s = 1.96 TeV, a significantly different value of |t 0 |(pp) = 0.70 ± 0.05 GeV 2 was obtained for the position of this point. The difference between |t 0 |(pp) and |t 0 |(pp) is apparently a clear but subtle Odderon effect, that cannot be obviously obtained by directly inspecting Fig. 19, but it supports the same conclusion about the presence of Odderon effects in the few TeV elastic scattering data.
As a by-product, but perhaps even more importantly, we have also found a clear-cut evidence for a proton substructure, as shown by the presence of the second, nearly exponential region in the differential cross-sections of elastic pp collisions at large |t|. At the ISR region, 23.5 ≤ √ s ≤ 62.5 GeV, from the asymptotic value of the t-dependent slope parameter of B ISR ≈ 2 GeV −2 a substructure with a Gaussian radius of R ISR,G ≈ 0.6 fm, while at the LHC energies of 7 ≤ √ s ≤ 13 TeV, a substructure with a different Gaussian radius of R LHC,G ≈ 0.9 fm is identified. Apparently, the size of this structure found at the ISR and LHC matches reasonably well the size of the dressed quarks and diquarks, respectively, as found recently in a unitarized Bialas-Bzak model analysis of 7 TeV elastic pp scattering [35]. These results may provide a phenomenological support for the quark-diquark picture of hadron confinement as obtained recently by Brodsky using AdS/QCD techniques [44], although they may not exclude other reasonable interpretations.
Our analysis indicates that the proton substructure contributes to the total pp crosssections with σ ISR = 0.3 +0. 3 −0.1 mb at ISR, and σ LHC ≈ 8.2 +7.9 −4.7 mb at the LHC energies. We have also found that this substructure can be better characterized by a Lévy source with α = 0.9 as compared to a Gaussian source (corresponding to α = 1). Using the characteristic Lévy length scale R L , we find that these substructures of the protons are characterized of of R ISR,L = 0.3 ± 0.1 fm at ISR, and R LHC,L = 0.5 ± 0.1 fm at the LHC energies.
From the analysis of the cone region, we clearly demonstrated that the shape of the protons actually changes in the 7 -13 TeV energy range, corresponding to an opening of a new channel, as clearly demonstrated by the appearance of a saturated P (b) ≈ 1 region in the shadow profile functions in the TeV energy range.
Finally, based on our experience with precision description of the differential cross-sections of elastic pp and pp collisions let us warn the careful readers against over-interpreting reasonably looking fit results in cases when the fitted function does not represent the data with a statistically not unacceptable confidence level.
We hope that this data analysis method of Lévy series expansion, detailed for the first time in this manuscript for a positive definite function, may find several important applications in the future, in a broad range of quantitative sciences. Essentially this method is able to characterize the deviations from Fourier-transformed and symmetric Lévy stable source distributions. Given the ubiquity of Lévy distributions in Nature, we hope that our new method will be relevant in several areas of human knowledge, that extend far beyond the science of physics.    Some of the features are common for all of these plots.
1) The parameter α was found to be within the errors independent of √ s, so it was fixed to a collision energy independent constant value α = 0.9 in these cases. It is an interesting and open problem to search for a physics interpretation of the universality of the α value of elastic pp collisions, but this topic goes beyond the scope of the model-independent approach  FIG. 25: Second-order Lévy expansion fits of the differential cross-section data on pp elastic scattering at √ s = 53 GeV [46], with fixed α = 0.9.
followed in the current manuscript. As an extra bonus, fixing α = 0.9 has resulted in a reasonably good reproduction of the √ s dependence of the total pp cross-sections, achieved without an additional tuning.
2) These pp elastic scattering data were less detailed as compared to the corresponding  Third-order Lévy expansion fits of the differential cross-section data on pp elastic scattering at √ s = 630 GeV [48], with fixed α = 0.9.
pp data, presented in Appendix A. The "dip" region was covered in all the cases, however, at very low |t| as well as at large |t| values, the acceptance and thus the t range were rather limited. This prevented us from a reliable analysis of the ρ parameter of elastic scattering for these collisions even at t = 0.
3) We have tested if the third-order and fourth-order Lévy expansions give similar results for these data sets or not. Within errors, the fourth-order expansion parameters were found to be consistent with zero, so we have fixed them to zero and investigated the third-order Lévy expansion results. In some cases, e.g. for √ s = 53 and 1960 GeV, the second-order Lévy fits were employed as providing better confidence levels for these rather scarce data sets. 4) We have checked that ρ(t = 0) changes well outside the Minuit indicated errors, if we change the order of the expansion from a third-order to a fourth-order Lévy expansion. Usually high precision data at low values of |t| are needed to reconstruct ρ reliably, so this limitation is not unexpected. However, it indicates the sensitivity of our method. In contrast to ρ(t), the B(t) slope parameters and the P (b) shadow profiles were found to be stable with respect to increasing the order of the Lévy expansion from the third to the fourth order. 5) We have found that a third-order Lévy expansion provides a reasonable overall description for all of these four data sets: Minuit has converged, error matrix accurate, and the fit quality is not unacceptable. The confidence level is sufficiently large, CL > 0.1 % for all the four data sets. We concluded, that with the exception of the overly sensitive parameter ρ, we can interpret the physical meaning of the fit parameters, as they represent the data. Note that the |t| ranges of the various data sets are rather different, but it is clear that the Lévy expansion reproduces the elastic pp differential cross-section in the respective s and |t| ranges of these four data-sets.
The parameters of the Lévy expansion are shown on the corresponding figures, together with the extracted values of the total cross-section, and the measures of the fit quality. The value of the ρ parameter is also shown, but the fluctuation of these values as a function of √ s and their correlation with the Lévy exponent α also indicates that even the value of ρ(t = 0) cannot be reliably determined from the these pp elastic scattering data, due to lack of data points in the sufficiently low |t| region.
In this Appendix, we detail the Lévy fits, dσ/dt = A exp (−(R 2 t) α ) to the elastic scattering data of pp collisions for seven different data sets, at √ s = 23.5, 30.7, 44.7, 52.8, 62.5 GeV as well as at 7 and 13 TeV. These fits correspond to the zeroth order, leading terms of the Lévy expansions, with all the expansion coefficients set to zero or c i = 0, as detailed in Subsection II A. Some of the features are common for all of these plots.
In the ISR energy domain fits were performed in the 2.5 < |t| < 5 GeV 2 region. In each case, Minuit has converged, error matrix was accurate, and the confidence level of the fits was acceptable, with CL > 0.1 %. In each case, the value of the parameter α was, within errors, consistent with 0.9 so we have fixed its value to 0.9 in each case. Let us note, however, that the fixed value of this parameter α can be varied between 1.0 and 0.8 without changing the acceptability of the fit, so the criteria of CL > 0.1 % was satisfied.
The good quality of these fit results indicates, that a substructure is present inside the protons: the data after the dip and bump region have the same structure and quality, as a usual low-|t| elastic scattering data below the dip region, used to be parametrized by a nearly exponential shape. However, the related cross-sections and length-scales are smaller than that of the protons. This conclusion is rather obvious after one compares the results of the Lévy fits to the tail or large-|t| regions, as presented in Appendix C with the results of Lévy fit results in the cone or low-|t| region, as detailed in Appendix D below. These comparisons were summarized clearly in Fig. 6 for the tail and in Fig. 7 for the cone region. It is quite remarkable, that the proton substructure in the ISR region has a size that is apparently independent of the change of the energy in the region of a few tens of GeV, but this substructure appears to be different (in size) from a substructure that is emerged in pp collisions in the 7 -13 TeV energy range, as evidenced from the parallel dashed lines in Fig. 6. This observation and the lack of parallel behaviour between the ISR and the LHC data fits suggests that two substructures, characterized by different Lévy source radii but similar Lévy index of stability α are identified at the ISR and at the LHC energies.
We have also studied the stability of these fits. Increasing the value of α decreased the Lévy scale R and the contribution to the total cross-section. Taking into account a covariation of the fit parameters with the fixed value of α we find that in each of these data sets, we are allowed to vary the value of α in a reasonable range of 0.8 to 1.0. We found that in each case the same-size substructure of the proton was present, with a characteristic Lévy length scale of R ISR = 0.3 ± 0.1 fm, and a contribution to the total cross-section with σ ISR = 0.3 +0. 3 −0.1 mb, where the quoted errors take into account the errors coming from a variation of the value of α as well. Apparently, the same structure is seen in these reactions, within the systematic error of the analysis, as indicated in Figs. 29, 30, 31, 32, and 33. However, the same type of analysis reveals a different-size proton substructure, when one performs a similar analysis in the TeV energy range. Fits to the differential crosssections of elastic pp scattering data at √ s = 7 TeV and preliminary TOTEM data at √ s = 13 TeV are shown in Figs. 34, and 35, respectively. This substructure corresponds to a Lévy scale of R LHC = 0.5 ± 0.1 fm, and to a contribution to the total cross-section with σ LHC (7 TeV) = 6.1 +3.3 −2.6 mb at 7 TeV, which is within the errors the same as the preliminary value of σ LHC (13 TeV) = 10.2 +5.9 −4.7 mb at 13 TeV. It is clear that these scales of R and the contributions to the total cross-section at the level of σ LHC ≈ 8.2 +7.9 −4.7 mb, that emerge at the TeV scale colliding energies, are significantly larger than the corresponding cross-sections and scales in the ISR energy range of 23.5 -62.5 GeV. The errors of R LHC and σ LHC include  the estimated systematic uncertainties as well, and those are larger than the statistical errors shown on the figures. Indeed, they also take into account the fact, that similar results were obtained for α = 0.8 and α = 1.0 fixed values as well, where the increasing or decreasing of α value has resulted in a decrease or increase of R (and σ), respectively. The emergence of the proton substructure with two different sizes in elastic pp collisions  FIG. 32: Zeroth-order Lévy fits to the tail of the differential cross-section data on pp elastic scattering at √ s = 52.8 GeV [45], with fixed α = 0.9, A and R as free fit parameters.
at the ISR energy range of √ s = 23.5 -62.5 GeV and at the LHC energy range of √ s = 7 -13 TeV, respectively, is clearly demonstrated on the summary plot that shows the data and the Lévy fit results to the tails of the differential cross-sections of elastic pp collisions on the same plot, also indicating with dashed lines the extrapolation of the fit results outside the fitted domain, see Fig. 6.  The change of the slope of these lines from ISR to the LHC energies is so much obvious and striking to the naked eye, that it is rather surprising that such a change was not reported in the literature before, at least, to the best of our knowledge. Perhaps, the reason for this effect is the expectation, that perturbative QCD effects should start to be visible soon after the dip-bump structure which could result in a power-law decrease of the differential elastic  FIG. 36: Zeroth-order Lévy fits to the cone region of the differential cross-section data on pp elastic scattering at √ s = 23.5 GeV [45], with fixed α = 0.9, A and R being free fit parameters.
cross-section. Indeed, such a power-law fit was performed and published by the TOTEM Collaboration in Ref. [49] when analysing the first set of 7 TeV elastic scattering data. In particular, the TOTEM Collaboration has reported that for |t|-values larger than ∼1.5 GeV 2 , the differential cross-section of elastic pp collisions exhibits a power law behaviour, with an exponent of −7.8 ± 0.3 (stat) ± 0.1 (syst).   Given that the new (but still preliminary) TOTEM data at 13 TeV are much more detailed and cover an extended |t| range after the dip and bump structure, we can easily test the presence (or not) of the perturbative QCD predicted power-law tails in these elastic scattering data: if dσ/dt ∝ |t| n then the power-law tail emerges as a straight line on a log-log plot. In order to illustrate this point, we have prepared a log-log plot that includes the TOTEM preliminary dσ/dt data at √ s = 13 TeV, and the Lévy expansion fit results. As clearly seen in this Fig. 14, the tail is not a straight line on this plot, so the |t|-range is apparently not yet in the domain of perturbative QCD, with a possible exception for data points close to the end of the acceptance region at large |t|. These data points, however, have rather big error bars.
Appendix D: Lévy fits to elastic pp collisions at small |t| In this Appendix, we describe the zeroth-order Lévy fits, dσ/dt = A exp (−(R 2 t) α ) to the low-|t| or cone region of elastic pp scattering data for seven different data sets, at √ s = 23.5, 30.7, 44.7, 52.8, 62.5 GeV as well as at 7 and 13 TeV.
Traditionally, an exponential behaviour is assumed for this kinematic region, with α = 1 fixed fits. However, recently the TOTEM collaboration demonstrated a non-exponential behaviour in √ s = 8 TeV pp collisions, that was found to be significant, corresponding to a more than 7σ effect [25]. These results were highlighted and some of their implications and ramifications were detailed also in Ref. [51]. In this Appendix, we re-examine the other, already published data sets using the zeroth-order Lévy fits, where the selected low-|t| regions are very similar to the 8 TeV analysis of TOTEM, namely we fit a |t| region that is similar to the TOTEM analysis at 8 TeV, when |t| is measured in units of t dip . Our fit region is thus 0.04 |t dip | ≤ |t| ≤ 0.4 |t dip |. We have performed a detailed investigation of this region and found that all the data sets were described with a good confidence level with the zeroth-order Lévy fit, with only three free parameters, A, R and α. Within the errors, in the ISR region, the values of the parameter α were in the region of 0.90 ± 0.02 in these fits. FIG. 40: Zeroth-order Lévy fits to the cone region of the differential cross-section data on pp elastic scattering at √ s = 62.5 GeV [45], with fixed α = 0.9, A and R being free fit parameters.
Hence, we have fixed the value of α to 0.9, that allowed us to demonstrate the trends more clearly and also to compare the results with the similar analysis in the large |t| region.   FIG. 41: Zeroth-order Lévy fits to the cone region of the differential cross-section data on pp elastic scattering at √ s = 7 TeV, with fixed α = 0.9, A and R being free fit parameters.
All the data in the ISR energy range were fitted successfully with a CL ≥ 0.1 % for the α = 0.9 fixed case as well. The Lévy radii R kept on increasing monotonically with increasing √ s. The best fit parameters as well as the fit quality measures are indicated on Figs. 36,37,38,39,40. However, the fixed α = 0.9 fits failed both on the √ s = 7 TeV final and on the √ s = 13 TeV preliminary data, as indicated on Figs. 41 and 42. These fit results, although on a qualitative level apparently may be called as "reasonable", are in fact statistically not acceptable. This implies that in the TeV energy region, some new mechanism starts to work, that changes not only the Lévy scale but also the shape of the proton. Actually, from the analysis of the shadow profile function, we know that such a new effect in the TeV energy range corresponds to the saturation of the shadow profile functions P (b) in the region of b ≤ 0.4 − 0.5 fm. This result, detailed in Subsection III B is supported strongly and independently by the fits detailed above.