1 Introduction

Elastic scattering of protons is a process mediated by the strong and the electromagnetic interactions – the weak interaction is commonly neglected since its carriers are heavy compared to the small momentum transfers, |t|, typical of elastic scattering. In this context, the strong interaction is traditionally called ‘nuclear’ or ‘hadronic’ and the electromagnetic one ‘Coulomb’. In quantum-theory description, each of these interactions is described by a scattering amplitude, nuclear \(\mathcal{A}^\mathrm{N}(t)\) and Coulomb \(\mathcal{A}^\mathrm{C}(t)\). Moreover, the combined scattering amplitude receives a third contribution reflecting Feynman diagrams with both strong and electromagnetic exchanges. This term, together with the complex character of the scattering amplitudes, describes the effects of Coulomb–nuclear interference (CNI) in the differential cross-section. Since the Coulomb amplitude is known, measuring the CNI gives access to the phase of the nuclear amplitude, which is necessary for a complete understanding of the interaction but not directly observable in the pure hadronic differential cross-section. The CNI effect is most pronounced in the t-region where the two amplitudes have similar magnitudes, i.e. – for typical LHC centre-of-mass energies of a few TeV – near \(|t| \sim 5 \times 10^{-4}\,\mathrm{GeV}^{2}\). Thus the experimental sensitivity to the nuclear phase, \(\arg \mathcal{A}^\mathrm{N}(t)\), is limited to a region at very small |t|, making difficult any conclusions on the functional form of the phase.

Fig. 1
figure 1

Left Schematic view of the RP stations on both sides of IP5 with two proton tracks from an elastic event. Right Schematic view of the silicon detector positions in an RP station with a track traversing the overlap zone between top and horizontal detectors, providing detector alignment information

In the analyses of past experiments – see e.g. [18] (ISR), [9, 10] (\(\mathrm{S}\mathrm{\bar{p}}\mathrm{p}S\)), [11] (Tevatron) – a simplified interference formula was used. This so-called simplified West–Yennie (SWY) formula [12] is based on restrictive assumptions on the hadronic amplitude, implying in particular a purely exponential modulus and a constant phase for all t (see the discussion in Sect. 6.1.4). As a representative quantity, the phase value at \(t=0\), or equivalently

$$\begin{aligned} \rho \equiv \cot \, \arg \mathcal{A}^\mathrm{N} (0) = \frac{\mathfrak {R}\mathcal{A}^\mathrm{N} (0)}{\mathfrak {I}\mathcal{A}^\mathrm{N} (0)} \end{aligned}$$
(1)

was traditionally quoted. An interesting aspect of \(\rho \) is its predictive power in extrapolating the total cross-section to higher centre-of-mass energies via dispersion relations [13].

The present article discusses the first measurement of elastic scattering in the CNI region at the CERN LHC by the TOTEM experiment. The data have been collected at \(\sqrt{s} = 8\,\)TeV with a special beam optics (\(\beta ^{*}=1000\,\)m) and cover a |t|-interval from \(6\times 10^{-4}\) to 0.2 GeV\(^{2}\), extending well into the interference region. In order to strengthen the statistical power and thus enable a cleaner identification of the interference effects, the analysis also exploits another, complementary data set with higher statistics [14], taken at the same energy, but with different beam optics (\(\beta ^{*}=90\,\)m), and thus covering a different |t|-range: \(0.027< |t| < 0.2\,\mathrm{GeV^2}\). The isolated analysis of the latter data set has excluded a purely exponential behaviour of the observed elastic cross-section with more than \(7\,\mathrm{\sigma }\) confidence. The new data in the CNI region allow to study the source of the non-exponentiality: nuclear component, CNI effects or both. In order to explore the full spectrum of possibilities, an interference formula without the limitations of SWY is needed. In the present study the more general and complex interference formulae of Cahn [15] and Kundrát–Lokajíček (KL) [16] are used, offering much more freedom for the choice of the theoretically unknown functional forms of the hadronic modulus and phase. Since the data cannot unambiguously determine all functional forms and their parameters, the results of this study, still representatively expressed in terms of \(\rho \), become conditional to the choice of the model describing the hadronic amplitude. This choice has implications on the behaviour of the interaction in impact parameter space. In particular, the functional form of the hadronic phase at small |t| determines whether elastic collisions occur predominantly at small or large impact parameters (centrality vs. peripherality). It will be shown that both options are compatible with the data, thus the central picture still prevalent in theoretical models is not a necessity.

Section 2 of this article outlines the experimental setup used for the measurement. The properties of the special beam optics are described in Sect. 3. Section 4 gives details of the data-taking conditions. The data analysis and reconstruction of the differential cross-section are described in Sect. 5. Section 6 presents the study of the Coulomb–nuclear interference together with the functional form of the hadronic amplitude. The values of \(\rho \) and \(\sigma _\mathrm{tot}\) are determined.

2 Experimental apparatus

The TOTEM experiment, located at the LHC Interaction Point (IP) 5 together with the CMS experiment, is dedicated to the measurement of the total cross-section, elastic scattering and diffractive processes. The experimental apparatus, symmetric with respect to the IP, is composed of a forward proton spectrometer (Roman Pots, RPs) and the forward tracking telescopes T1 and T2. A complete description of the TOTEM detector instrumentation and its performance is given in [17, 18]. The data analysed here come from the RPs only. An RP is a movable beam-pipe insertion capable of approaching the LHC beam to a distance of less than a millimetre, in order to detect protons with scattering angles of only a few microradians. The proton spectrometer is organised in two RP stations: one on the left side of the IP (LHC sector 45) and one on the right (LHC sector 56), see Fig. 1 (left). Each RP station, located between 215 and 220 m from the IP, is composed of two units: “near” (215 m from the IP) and “far” (220 m). A unit consists of 3 RPs, one approaching the outgoing beam from the top, one from the bottom, and one horizontally. Each RP houses a stack of 5 “U” and 5 “V” silicon strip detectors, where “U” and “V” refer to two mutually perpendicular strip orientations. The sensors were designed with the specific objective of reducing the insensitive area at the edge facing the beam to only a few tens of micrometers. Due to the 5 m long lever arm between the near and the far RP units the local track angles can be reconstructed with a precision of about \(10\,\upmu \)rad. A high trigger efficiency (>99 %) is achieved by using all RPs independently. Since elastic scattering events consist of two collinear protons emitted in opposite directions, the detected events can have two topologies, called diagonals: 45 bottom–56 top and 45 top–56 bottom.

Fig. 2
figure 2

Trigger rates as a function of time from the beginning of the run (24 October 2012, 23:00 hours). The T2 rate (blue) is roughly proportional to luminosity, while the RP rate (red) is in addition sensitive to beam-halo level. The grey bands represent periods of uninterrupted data taking

This article uses a reference frame where x denotes the horizontal axis (pointing out of the LHC ring), y the vertical axis (pointing against gravity) and z the beam axis (in the clockwise direction).

3 Beam optics

The beam optics relates the proton kinematical states at the IP and at the RP location. A proton emerging from the interaction vertex \((x^*\), \(y^*)\) at the angle \((\theta _x^*,\theta _y^*)\) (relative to the z axis) and with momentum \(p\,(1+\xi )\), where p is the nominal initial-state proton momentum, is transported along the outgoing beam through the LHC magnets. It arrives at the RPs in the transverse position

$$\begin{aligned} x(z_\mathrm{RP})&= L_x(z_\mathrm{RP})\, \theta _x^*\ +\ v_x(z_\mathrm{RP})\, x^*\ +\ D_x(z_\mathrm{RP})\, \xi ,\\y(z_\mathrm{RP})&= L_y(z_\mathrm{RP})\, \theta _y^*\ +\ v_y(z_\mathrm{RP})\, y^*\ +\ D_y(z_\mathrm{RP})\, \xi \quad \end{aligned}$$
(2)

relative to the beam centre. This position is determined by the optical functions, characterising the transport of protons in the beam line and controlled via the LHC magnet currents. The effective length \(L_{x,y}(z)\), magnification \(v_{x,y}(z)\) and dispersion \(D_{x,y}(z)\) quantify the sensitivity of the measured proton position to the scattering angle, vertex position and momentum loss, respectively. Note that for elastic collisions the dispersion terms \(D\,\xi \) can be ignored because the protons do not lose any momentum. The values of \(\xi \) only account for the initial state momentum offset and variations, see Section 4 in [14]. Due to the collinearity of the two elastically scattered protons and the symmetry of the optics, the impact of \(D\,\xi \) on the reconstructed scattering angles is negligible compared to other uncertainties.

The data for the analysis presented here have been taken with a new, special optics, conventionally labelled by the value of the \(\beta \)-function at the interaction point, \(\beta ^{*} = 1000\,\)m, and specifically developed for measuring low-|t| elastic scattering. It maximises the vertical effective length \(L_{y}\) at the RP position \(z = 220\,\mathrm{m}\) and minimises the vertical magnification \(|v_{y}|\) at \(z = 220\,\)m (Table 1). This configuration is called “parallel-to-point focussing” because all protons with the same angle in the IP are focussed on one point in the RP at 220 m. It optimises the sensitivity to the vertical projection of the scattering angle – and hence to |t| – while minimising the influence of the vertex position. In the horizontal projection the parallel-to-point focussing condition is not fulfilled, but – unlike in the \(\beta ^{*} = 90\,\)m optics used for previous measurements [1922] – the effective length \(L_{x}\) at \(z = 220\,\)m is non-zero, which reduces the uncertainty in the horizontal component of the scattering angle.

Table 1 Optical functions for elastic proton transport for the \(\beta ^{*} = 1000\,\)m optics. The values refer to the right arm, for the left one they are very similar

4 Data taking

The results reported here are based on data taken in October 2012 during a dedicated LHC proton fill (3216) with the special beam properties described in the previous section.

The vertical RPs approached the beam centre to only 3 times the beam width, \(\sigma _{y}\), resulting in an acceptance for |t|-values down to \(6 \times 10^{-4}\,\mathrm{GeV}^{2}\). The exceptionally close distance was possible due to the low beam intensity in this special beam operation: each beam contained only two colliding bunches and one non-colliding bunch for background monitoring, each with \(10^{11}\) protons. A novel collimation strategy was applied to keep the beam halo background under control. As a first step, the primary collimators (TCP) in the LHC betatron cleaning insertion (point 7) scraped the beam down to \(2\,\sigma _{y}\); then the collimators were retracted to \(2.5\,\sigma _{y}\), thus creating a \(0.5\,\sigma _{y}\) gap between the beam edge and the collimator jaws. With the halo strongly suppressed and no collimator producing showers by touching the beam, the RPs at \(3\,\sigma _{y}\) were operated in a background-depleted environment for about 1 h until the beam-to-collimator gap was refilled by diffusion, as diagnosed by the increasing RP trigger rate (Fig. 2). When the background conditions had deteriorated to an unacceptable level, the beam cleaning procedure was repeated, again followed by a quiet data-taking period. The beam cleaning at \(1.5\,\mathrm{h}\) from the beginning of the run employed only vertical collimators and led to a quickly increasing background rate, see Fig. 2. Therefore, the following beam cleaning operations were also performed in the horizontal plane. Altogether there were 6 beam cleaning interventions until the luminosity had decreased from initially \(1.8\times 10^{27}\) to \(0.4\times 10^{27}\,\mathrm{cm}^{-2}\,\mathrm{s}^{-1}\) at which point the data yield was considered as too low. During the 9 h long fill, an integrated luminosity of \(20\,\upmu \mathrm{b}^{-1}\) was accumulated in 6 data sets corresponding to the calm periods between the cleaning operations.

Due to an anti-collision protection system, the top and the bottom pots of a vertical RP unit could not approach each other close enough to be both at a distance of \(3\,\sigma _{y} = 780\,\upmu \)m from the beam centre. Therefore a configuration with one RP diagonal (45 top–56 bottom) at \(3\,\sigma _{y}\) (“close diagonal”) and the other (45 bottom–56 top) at \(10\,\sigma _{y}\) (“distant diagonal”) was chosen. The distant diagonal provides a systematic comparison at larger |t|-values. The horizontal RPs were only needed for the data-based alignment and therefore placed at a safe distance of \(10\,\sigma _{x} \approx 7.5\) mm, close enough to have an overlap with the vertical RPs (Fig. 1, right).

The events collected were triggered by a logical OR of: inelastic trigger (at least one charged particle in either arm of T2), double-arm proton trigger (coincidence of any RP left of IP5 and any RP right of IP5) and zero-bias trigger (random bunch crossings) for calibration purposes.

In the close and distant diagonals a total of 190 and 162 k elastic event candidates have been tagged, respectively.

5 Differential cross-section

The analysis method is very similar to the previously published one [14]. Section 5.1 covers all aspects related to the reconstruction of a single event. Section 5.2 describes the steps of transforming a raw t-distribution into the differential cross-section. The t-distributions for the two diagonals are analysed separately. After comparison (Sect. 5.3) they are finally merged (Sect. 5.4).

5.1 Event analysis

The event kinematics are determined from the coordinates of track hits in the RPs after proper alignment (see Sect. 5.1.2) using the LHC optics (see Sect. 5.1.3).

5.1.1 Kinematics reconstruction

For each event candidate the scattering angles and vertex positions of both protons (one per arm) are first determined separately by inverting the proton transport equation (2), assuming \(\xi = 0\):

$$\begin{aligned} \theta _x^{*\mathrm L,R}&= {v_x^\mathrm{N} x^\mathrm{F} - v_x^\mathrm{F} x^\mathrm{N}\over v_x^\mathrm{N} L_x^\mathrm{F} - v_x^\mathrm{F} L_x^\mathrm{N}} ,\quad \theta _y^{*\mathrm L,R} = {1\over 2} \left( {y^\mathrm{N}\over L_y^\mathrm{N}} + {y^\mathrm{F}\over L_y^\mathrm{F}} \right) ,\nonumber \\ x^{*\mathrm L,R}&= {L_x^\mathrm{N} x^\mathrm{F} - L_x^\mathrm{F} x^\mathrm{N}\over L_x^\mathrm{N} v_x^\mathrm{F} - L_x^\mathrm{F} v_x^\mathrm{N}} , \end{aligned}$$
(3)

where the N and F superscripts refer to the near and far units, L and R to the left and right arm, respectively. This one-arm reconstruction is used for tagging elastic events, where the left and right arm protons are compared.

Once a proton pair has been selected, all four RPs are used to reconstruct the kinematics of the event, optimising the angular resolution (see Sect. 5.1.4):

$$\begin{aligned} \theta _x^*= & {} { \sum {v_x^i}^2 \sum L_x^i x^i - \sum L_x^i v_x^i \sum v_x^i x^i \over \sum {v_x^i}^2 \sum {L_x^i}^2 - \sum L_x^i v_x^i \sum v_x^i L_x^i } ,\nonumber \\&\theta _y^* = {1\over 4} \sum {y^i\over L_y^i} , \end{aligned}$$
(4)

where the sums run over the superscript i representing the four RPs of a diagonal.

Eventually, the scattering angle, \(\theta ^*\), and the four-momentum transfer squared, t, are calculated:

$$\begin{aligned} \theta ^* = \sqrt{{\theta _x^*}^2 + {\theta _y^*}^2},\quad t = - p^2 ({\theta _x^*}^2 + {\theta _y^*}^2) , \end{aligned}$$
(5)

where p denotes the beam momentum.

5.1.2 Alignment

TOTEM’s usual three-stage procedure [18] for correcting the detector positions and rotation angles has been applied: a beam-based alignment prior to the run followed by two offline methods. First, track-based alignment for relative positions among RPs, and second, alignment with elastic events for absolute position with respect to the beam – repeated in 15 min time intervals to check for possible beam movements.

The offline procedure has been extended further to improve the vertical alignment. The new steps exploit the fact that elastic events with their two collinear protons relate the alignments in the left and right arm with an uncertainty of \(20\,\mathrm{\upmu m}\). Furthermore, the horizontal RPs in the right arm recorded a hit distribution usable for vertical alignment in addition to the standard technique based on the vertical RPs, see Fig. 3.

Fig. 3
figure 3

Hit scatter plot in the right-far unit, corresponding to a period of 15 min. The black dots represent track hits in the vertical and horizontal RPs. Track hits close to the sensor edges are removed because of possible bias due to acceptance effects. The green histogram shows the horizontal profile of hits in the vertical RPs, the dashed green line interpolates the profiles between the top and bottom RPs. Similarly, the red histogram gives the vertical profile of the hits in the horizontal RP and the dashed red line its extrapolation to the beam region. The blue dashed line indicates the vertical centre of symmetry of the hits in the vertical RPs (see [18] for more details). The crossing of the dashed lines represents the position the beam centre (the two vertical-alignment results are averaged)

Exploiting all the methods, the alignment uncertainties have been estimated to \(30\,\mathrm{\upmu m}\) (horizontal shift), \(70\,\mathrm{\upmu m}\) (vertical shift) and \(2\,\mathrm{m rad}\) (rotation about the beam axis). Propagating them through Eq. (4) to reconstructed scattering angles yields \(0.28\,\mathrm{\upmu rad}\) (\(0.19\,\mathrm{\upmu rad}\)) for the horizontal (vertical) angle. RP rotations induce a bias in the reconstructed horizontal scattering angle:

$$\begin{aligned} \theta _x^* \rightarrow \theta _x^* + c \theta _y^* , \end{aligned}$$
(6)

where the proportionality constant c has a mean of 0 and a standard deviation of 0.005.

5.1.3 Optics

It is crucial to know with high precision the LHC beam optics between IP5 and the RPs, i.e. the behaviour of the spectrometer composed of the various magnetic elements. The optics calibration has been applied as described in [23]. This method uses RP observables to determine fine corrections to the optical functions presented in Eq. (2).

The residual errors induce a bias in the reconstructed scattering angles:

$$\begin{aligned} \theta _x^* \rightarrow (1 + d_x)\, \theta _x^*,\quad \theta _y^* \rightarrow (1 + d_y)\, \theta _y^*. \end{aligned}$$
(7)

For the two-arm reconstruction, Eq. (4), the biases \(d_x\) and \(d_y\) have uncertainties of 0.34 and \(0.25\,\mathrm{\%}\), respectively, and a correlation factor of \(-0.89\). These estimates include the effects of magnet harmonics. To evaluate the impact on the t-distribution, it is convenient to decompose the correlated biases \(d_x\) and \(d_y\) into eigenvectors of the covariance matrix:

$$\begin{aligned} \begin{pmatrix} d_x\\ d_y \end{pmatrix} = \eta _1 \underbrace{\begin{pmatrix} +0.338\,\mathrm{\%} \\ -0.234\,\mathrm{\%} \end{pmatrix}}_\mathrm{mode\ 1} \ +\ \eta _2 \underbrace{\begin{pmatrix} -0.053\,\mathrm{\%} \\ -0.076\,\mathrm{\%} \end{pmatrix}}_\mathrm{mode\ 2} \end{aligned}$$
(8)

normalised such that the factors \(\eta _{1,2}\) have unit variance.

5.1.4 Resolution

Statistical fluctuations in \(\theta _y^*\) are mostly due to the beam divergence and can be studied by comparing the angles reconstructed from the left and right arm. As illustrated in Fig. 4, the distributions show only minimal deviations from a Gaussian shape. By dividing their standard deviation by a factor of 2, one can estimate the resolution of the two-arm reconstruction (Eq. (4)) of elastic events, see Fig. 5, bottom. Moreover, measurements of beam emittances [24] indicate that the vertical divergences of the two beams can be considered as equal with a tolerance of about \(25\,\mathrm{\%}\). Exploiting this fact, one can de-convolute the distribution of \(\theta _y^{*\mathrm R} - \theta _y^{*\mathrm L}\) in order to obtain the beam-divergence distribution, used for the acceptance corrections discussed in Sect. 5.2.3.

Fig. 4
figure 4

Difference between vertical scattering angles reconstructed in the right and left arm, for the diagonal 45 top–56 bottom. Red: data from run start (0 to \(1\,\mathrm{h}\) from the beginning of the run). Blue: data from run end (7–\(8\,\mathrm{h}\)), scaled by 0.1. The solid lines represent Gaussian fits

In the horizontal projection, a more complex procedure is used since the one-arm reconstruction, Eq. (3), is strongly influenced by the detector resolution. First, the horizontal beam divergence is estimated from the standard deviation of reconstructed vertices, \(\sigma (x^*)\):

$$\begin{aligned} \sigma ^\mathrm{bd}(\theta _x^*) = {\sigma (x^*) \sqrt{2}\over \beta ^*}. \end{aligned}$$
(9)

It increases from 0.75 to \(0.9\,\mathrm{\upmu rad}\) over the time of the fill. Subtracting this component from the standard deviation of \(\theta _x^{*\mathrm R} - \theta _x^{*\mathrm L}\), one determines the (mean) spatial resolution of the sensors in each diagonal: \(10.7\,\mathrm{\upmu m}\) (45 top–56 bottom) and \(12.1\,\mathrm{\upmu m}\) (45 bottom–56 top). These results have been verified to be time independent. Finally, the beam divergence and sensor resolution components can be propagated through Eq. (4) to estimate the \(\theta _x^*\) resolution for elastic events, as plotted in Fig. 5, top.

Fig. 5
figure 5

Angular resolution for the two-arm reconstruction, Eq. (4), as a function of time (from the beginning of the run). The blue (red) dots correspond to the diagonal 45 bottom–56 top (45 top–56 bottom). The grey bands indicate regions of uninterrupted data-taking, whereas in the remaining periods the beam cleaning procedure described in Sect. 4 was performed

5.2 Differential cross-section reconstruction

For a given t bin, the differential cross-section is evaluated by selecting and counting elastic events:

$$\begin{aligned} {\mathrm{d}\sigma \over \mathrm{d}t}(\hbox {bin}) = \mathcal{N}\, \mathcal{U}(\mathrm{bin})\, \mathcal{B}\, {1\over \Delta t} \sum \limits _{t\, \in \, \mathrm{bin}} \mathcal{A}(\theta ^*, \theta_y^*)\ \mathcal{E}(\theta_y^*) , \end{aligned}$$
(10)

where \(\Delta t\) is the width of the bin, \(\mathcal{N}\) is a normalisation factor and the other symbols stand for various correction factors: \(\mathcal{U}\) for unfolding of resolution effects, \(\mathcal{B}\) for background subtraction, \(\mathcal{A}\) for acceptance correction and \(\mathcal{E}\) for detection and reconstruction efficiency.

5.2.1 Event tagging

Table 2 The elastic selection cuts. The superscripts R and L refer to the right and left arm. The \(\alpha \theta _x^*\) term in cut 3 absorbs the effects of residual optics imperfections, \(\alpha \) is of the order of \(0.1\,\mathrm{\upmu m/\upmu rad}\). The right-most column gives a typical RMS of the cut distribution

The cuts used to select the elastic events are summarised in Table 2. Cuts 1 and 2 require the reconstructed-track collinearity between the left and right arm. Cut 3 ensures that the protons come from the same vertex (horizontally). The correlation plots corresponding to these cuts are shown in Fig. 6. Thanks to the very low beam divergence, the collinearity cuts are very powerful, and consequently other conceivable cuts (cf. Table 2 in [20]) bring no significant improvement.

Fig. 6
figure 6

Correlation plots for the event selection cuts presented in Table 2, showing events with diagonal topology 45 top–56 bottom. The solid black lines delimit the signal region within \(\pm 4\,\mathrm{\sigma }\)

Since a Monte-Carlo study shows that applying the three cuts at the \(3\,\mathrm{\sigma }\) level would lead to a loss of about \(0.5\,\mathrm{\%}\) of the elastic events, the cut threshold is set to \(4\,\mathrm{\sigma }\).

The tagging efficiency has been studied by applying the cuts also at the \(5\,\mathrm{\sigma }\)-level. This selection has yielded \(0.3\,\mathrm{\%}\) more events in every |t|-bin. This kind of inefficiency only contributes to a global scale factor, which is irrelevant for this analysis because the normalisation is taken from a different data set (cf. Sect. 5.2.6).

Fig. 7
figure 7

Distributions of discriminator 1, i.e. the difference between the horizontal scattering angle reconstructed from the right and the left arm. Solid curves: data from diagonal 45 top–56 bottom, the different colours correspond to various combinations of the selection cuts (see numbering in Table 2). Dashed curves: data from anti-diagonal RP configurations, obtained by inverting track coordinates in the left arm. The vertical dashed lines represent the boundaries of the signal region (\(\pm 4\,\mathrm{\sigma }\))

5.2.2 Background

As the RPs were very close to the beam, one may expect an enhanced background from coincidence of beam halo protons hitting detectors in the two arms. Other background sources (pertinent to any elastic analysis) are: central diffraction and pile-up of two single diffraction events.

The background rate (i.e. impurity of the elastic tagging) is estimated in two steps, both based on distributions of discriminators from Table 2 plotted in various situations, see an example in Fig. 7. In the first step, diagonal data are studied under several cut combinations. While the central part (signal) remains essentially constant, the tails (background) are strongly suppressed when the number of cuts is increased. In the second step, the background distribution is interpolated from the tails into the signal region. The form of the interpolation is inferred from non-diagonal RP track configurations (45 bottom–56 bottom or 45 top–56 top), artificially treated like diagonal signatures by inverting the coordinate signs in the arm 45; see the dashed distributions in the figure. These non-diagonal configurations cannot contain any elastic signal and hence consist purely of background which is expected to be similar in the diagonal and non-diagonal configurations. This expectation is supported by the agreement of the tails of the blue solid and dashed curves in the figure. Since the non-diagonal distributions are flat, the comparison of the signal-peak size to the amount of interpolated background yields the estimate \(1 - \mathcal{B} < 10^{-4}\).

5.2.3 Acceptance correction

The acceptance of elastic protons is limited by two factors: sensor coverage (relevant for low \(|\theta ^*_y|\)) and LHC beam aperture (at \(|\theta ^*_y| \approx 100\,\mathrm{\upmu rad}\)). Moreover, there is a region in the kinematic parameter space where elastic protons may interact with the horizontal RPs leading to uncertain detection efficiency. To avoid this region, an additional fiducial cut has been adopted: \(-50< \theta _x^* < 80\,\mathrm{\upmu rad}\). In the far vertical RPs, this restriction corresponds to about \(-2.3< x < 3.7\,\mathrm{mm}\). All acceptance related cuts are visualised in Fig. 8.

Fig. 8
figure 8

Distribution of scattering angle projections \(\theta _y^*\) vs. \(\theta _x^*\). The upper (lower) part comes from the diagonal 45 bottom–56 top (45 top–56 bottom). The red horizontal lines represent cuts due to the LHC apertures, the blue horizontal lines cuts due to the sensor edges. The vertical magenta lines delimit the fiducial region with detection efficiency not affected by the horizontal RPs. The dotted circles show contours of constant scattering angle \(\theta ^*\) as indicated in the middle of the plot (values in micro-radians). The parts of the contours within acceptance are emphasized in thick black

The correction for the above limitations includes two contributions – a geometrical correction \(\mathcal{A}_\mathrm{geom}\) reflecting the fraction of the phase space within the acceptance and a component \(\mathcal{A}_\mathrm{fluct}\) correcting for fluctuations around the vertical acceptance limitations:

$$\begin{aligned} \mathcal{A}(\theta ^*, \theta _y^*) = \mathcal{A}_\mathrm{geom}(\theta ^*)\ \mathcal{A}_\mathrm{fluct}(\theta _y^*). \end{aligned}$$
(11)

The fiducial cuts in \(\theta _x^*\) have been given sufficient margin from the region with uncertain efficiency to render the respective fluctuation correction negligible.

The calculation of the geometrical correction \(\mathcal{A}_\mathrm{geom}\) is based on the azimuthal symmetry of elastic scattering, experimentally verified for the data within acceptance. As shown in Fig. 8, for a given value of \(\theta ^*\) the correction is given by:

$$\begin{aligned} \mathcal{A_\mathrm{geom}}(\theta ^*) = { \hbox {full circumference}\over \hbox {arc length within acceptance} }. \end{aligned}$$
(12)

The correction \(\mathcal{A}_\mathrm{fluct}\) is calculated analytically from the probability that any of the two elastic protons leaves the region of acceptance due to the vertical beam divergence. The beam divergence distribution is modelled as a Gaussian with the spread determined by the method described in Sect. 5.1.4. This contribution is sizeable only close to the acceptance limitations. Data from regions with corrections larger than 2.5 are discarded. The uncertainties are related to the resolution parameters. For the lowest |t| bin their relative values are: vertical beam divergence: \(2\,\mathrm{\%}\), left-right asymmetry: \(1\,\mathrm{\%}\), and non-Gaussian shape: \(1\,\mathrm{\%}\).

Figure 9 shows an example of the t-dependence of the acceptance correction for the diagonal reaching lower |t|-values. Since a single diagonal cannot cover more than half of the phase space, the minimum value of the correction is 2. The very low |t| data points with the full correction larger than 10 are discarded to avoid biases. At the high-|t| end all data points are kept.

Fig. 9
figure 9

Full acceptance correction, \(\mathcal{A}\), for diagonal 45 bottom–56 top. The points give the mean value per bin, the error bars indicate the standard deviation. The abrupt changes in the shape correspond to acceptance cuts as indicated by the arrows

5.2.4 Inefficiency corrections

Since the overall normalisation will be determined from another dataset (see Sect. 5.2.6), any inefficiency correction that does not alter the t-distribution shape does not need to be considered in this analysis (trigger, data acquisition and pile-up inefficiency discussed in [20, 22]). The remaining inefficiencies are related to the inability of a RP to resolve the elastic proton track.

One such case is when a single RP does not detect and/or reconstruct a proton track, with no correlation to other RPs. This type of inefficiency, \(\mathcal{I}_{3/4}\), is evaluated by removing the RP from the tagging cuts (Table 2), repeating the event selection and calculating the fraction of recovered events. A typical example is given in Fig. 10, showing that the efficiency decreases gently with the vertical scattering angle. This dependence originates from the fact that protons with larger \(|\theta _y^*|\) hit the RPs further from their edge and therefore the potentially created secondary particles have more chance to induce additional signal. Since the RP detectors cannot resolve multiple tracks (non-unique association between “U” and “V” track candidates), a secondary particle track prevents from using the affected RP in the analysis.

Fig. 10
figure 10

Single-RP uncorrelated inefficiency for the near bottom RP in the right arm. The rapid drop at \(\theta _y^* \approx 8\,\mathrm{\mu rad}\) is due to acceptance effects at the sensor edge. The red lines represent a linear fit of the efficiency dependence on the vertical scattering angle (solid) and its extrapolation to the regions affected by acceptance effects (dashed)

Another source of inefficiency are proton interactions in a near RP affecting simultaneously the far RP downstream. The contribution from these near-far correlated inefficiencies, \(\mathcal{I}_{2/4}\), is determined by evaluating the rate of events with high track multiplicity (\(\gtrsim \)5) in both near and far RPs. Events with high track multiplicity simultaneously in a near top and near bottom RP are discarded as such a shower is likely to have started upstream from the RP station and thus be unrelated to the elastic proton interacting with detectors. The outcome, \(\mathcal{I}_{2/4} \approx 1.5\,\mathrm{\%}\), is compatible between left/right arms and top/bottom RP pairs and compares well to Monte-Carlo simulations (e.g. Section 7.5 in [25]).

The full correction is calculated as

$$\begin{aligned} \mathcal{E}(\theta _y^*) = {1\over 1 - \left( \sum \limits _{i\in \mathrm RPs} \mathcal{I}^i_{3/4}(\theta _y^*) + 2 \mathcal{I}_{2/4} \right) }. \end{aligned}$$
(13)

The first term in the parentheses sums the contributions from the four RPs of a diagonal and increases from about 16 to \(18\,\mathrm{\%}\) from the lowest to the highest \(|\theta _y^*|\). These values are higher than in the previous analyses (e.g. Section 5.2.4 in [14]) due the contribution from the far RPs in the left arm. The reconstruction efficiency in these pots is decreased by showers initiated by beam halo protons in the horizontal RPs upstream (closer to the beam in the left arm than in the right one).

5.2.5 Unfolding of resolution effects

Due to the very small beam divergence, the correction for resolution effects can be safely determined by the following iterative procedure.

  1. 1.

    The differential cross-section data are fitted by a smooth curve.

  2. 2.

    The fit is used in a numerical-integration calculation of the smeared t-distribution (using the resolution parameters determined in Sect. 5.1.4). The ratio between the smeared and the non-smeared t-distributions gives a set of per-bin correction factors.

  3. 3.

    The corrections are applied to the observed (yet uncorrected) differential cross-section yielding a better estimate of the true t-distribution.

  4. 4.

    The corrected differential cross-section is fed back to step 1.

As the estimate of the true t-distribution improves, the difference between the correction factors obtained in two successive iterations decreases. When the difference becomes negligible, the iteration stops. This is typically achieved after the second iteration.

The final correction is negligible (\(\mathcal{U} \approx 1\)) for all bins except at very low |t| where the rapid cross-section growth occurs, see Fig. 11.

Fig. 11
figure 11

Unfolding correction for the close diagonal (45 bottom–56 top). The vertical dashed line indicates the position of the acceptance cut due to sensor edges

For the uncertainty estimate, the uncertainties of the \(\theta _x^*\) and \(\theta _y^*\) resolutions (accommodating the full time variation) as well as fit-model dependence have been considered, each contribution giving a few per-mille for the lowest-|t| bin.

5.2.6 Normalisation

The normalisation \(\mathcal{N}\) is determined by requiring the same cross-section integral between \(|t| = 0.014\) and \(0.203\,\mathrm{GeV^2}\) as for dataset 1 published in [22]. This publication describes a measurement of elastic and inelastic rates at the same collision energy of \(\sqrt{s} = 8\,\mathrm{TeV}\). These rates can be combined using the optical theorem in order to resolve the value of the luminosity which consequently allows for normalisation of the differential cross-section. The leading uncertainty of \(\mathcal{N}\), \(4.2\,\mathrm{\%}\), comes from the rate uncertainties in [22].

5.2.7 Binning

At very low |t|, where the cross-section varies the fastest (\(\approx \) \(0.001\,\mathrm{GeV^2}\)), a fine binning is used. In the middle of the |t| range (\(\approx \) \(0.03\,\mathrm{GeV^2}\)), the bin width is chosen to give about \(1\,\mathrm{\%}\) statistical uncertainty. This rule is abandoned at higher |t| (above \(0.07\,\mathrm{GeV^2}\)) in favour of bins with a constant width of \(0.01\,\mathrm{GeV^2}\) to avoid excessively large bins.

5.2.8 Systematic uncertainties

Besides the systematic uncertainties mentioned at the above analysis steps, the beam momentum uncertainty needs to be considered when the scattering angles are translated to t, see Eq. (5). The uncertainty was estimated to \(0.1\,\mathrm{\%}\) in Section 5.2.8 in [14] which is further supported by a recent review [26].

Two different methods are used to propagate the systematic effects into the t-distribution. The first is based on a Monte-Carlo simulation which uses a fit of the final differential cross-section data to generate the true t-distribution. In parallel, another t-distribution is built, introducing one of the above mentioned systematic effects at \(1\,\mathrm{\sigma }\) level. The difference between the two t-distributions gives the systematic effect on the differential cross-section. The second method is similar, however using numerical integration techniques instead of Monte-Carlo simulations. Both methods are formally equivalent to evaluating

$$\begin{aligned} \delta s_{q}(t) \equiv \frac{\partial (\mathrm{d}\sigma /\mathrm{d}t)}{\partial q}\ \delta q , \end{aligned}$$
(14)

where \(\delta q\) corresponds to a \(1\,\mathrm{\sigma }\) bias in the quantity q responsible for a given systematic effect.

The Monte-Carlo simulations show that the combined effect of several systematic errors is well approximated by linear combination of the individual contributions from Eq. (14).

5.3 Systematic cross-checks

Compatible results have been obtained by analysing data subsets of events from different bunches, different diagonals and different time periods – in particular those right after and right before the beam cleanings.

5.4 Final data merging

Finally, the differential cross-section histograms from both diagonals are merged. This is accomplished by a per-bin weighted average, with the weight given by inverse squared statistical uncertainty. The statistical and systematic uncertainties are propagated accordingly. For the systematic ones, the correlation between the diagonals is taken into account. For example the vertical (mis-)alignment of the RPs of one unit is almost fully correlated; thus the effect on the differential cross-section is opposite for the two diagonals and consequently its impact is strongly reduced once the diagonals are merged.

The cross-section values can be found in Table 3 and visualised in Fig. 12. The figure clearly shows a rapid cross-section rise below \(|t| \lesssim 0.002\,\mathrm{GeV^2}\), which will later be interpreted as an effect due to electromagnetic interaction.

Table 3 The elastic differential cross-section as determined in this analysis. The three left-most columns describe the bins in t. The representative point gives the t value suitable for fitting [27]. The other columns are related to the differential cross-section. The six right-most columns give the leading systematic biases in \(\mathrm{d}\sigma /\mathrm{d}t\) for \(1\sigma \)-shifts in the respective quantities, \(\delta s_q\), see Eqs. (14) and (15). The two contributions due to optics correspond to the two vectors in Eq. (8)
Fig. 12
figure 12

Differential cross-section from Table 3 with statistical (bars) and systematic uncertainties (bands). The yellow band represents all systematic uncertainties, the brown one all but normalisation. The bands are centred around a data fit including both nuclear and Coulomb components (Eqs. (25), (18) and (17) with \(N_b=3\)). Inset: a low-|t| zoom featuring cross-section rise due to the Coulomb interaction

The final systematic uncertainties, except the \(4.2\,\mathrm{\%}\) coming from the normalisation, are summarised in Fig. 13 where their impact on the differential cross-section is shown. The leading uncertainties include normalisation, optics imperfections, beam momentum offset and residual misalignment. The vertical misalignment is the dominant systematic effect in the very-low |t| region. The leading effects are quantified in Table 3 and can be used to approximate the covariance matrix of systematic uncertainties:

$$\begin{aligned} \mathsf{V}_{ij} = \sum _{q} \delta s_{q}(i)\ \delta s_{q}(j), \end{aligned}$$
(15)

where i and j are bin indices (row numbers in Table 3) and the sum goes over the leading error contributions q (six right-most columns in the table).

Fig. 13
figure 13

Impact of t-dependent systematic effects on the differential cross-section. Each curve corresponds to a systematic error of \(1\,\mathrm{\sigma }\), cf. Eq. (14). The two contributions due to the optics correspond to the two vectors in Eq. (8). The envelope is determined by summing all shown contributions in quadrature for each |t| value

Let us emphasize that the systematic effects with linear t dependence (see Fig. 13) cannot alter the non-purely-exponential character of the data. This is the case for the effects of normalisation, beam momentum and to a large degree also of optics-mode 2. For the beam momentum, this can also be understood analytically: changing the value of p would yield a scaling of t, see Eq. (5), and consequently also scaling of the \(b_n\) parameters in Eq. (17). However, the non-zero \(b_2, b_3\) etc. parameters (reflecting the non-exponentiality) cannot be brought to 0 (as in purely-exponential case).

6 Coulomb–nuclear interference

The Coulomb–nuclear interference (CNI) can be used to probe the nuclear component of the scattering amplitude. Since the CNI effects are sensitive to the phase of the nuclear amplitude, both modulus and phase can be tested.

For the modulus, a relevant question is whether the earlier reported non-exponentiality of the differential cross-section [14] can be attributed solely to the nuclear component or whether Coulomb scattering gives a sizeable contribution. Concerning the phase, several parametrisations with different physics interpretations will be tested; for each of them the \(\rho \) parameter (representative for the phase value at \(t = 0\) according to Eq. (1)) will be determined.

Section 6.1 outlines the theoretical concepts needed to describe the CNI effects. Section 6.2 provides details on fitting procedures used to analyse the data. Sections 6.3 and 6.4 discuss the fit results for two relevant alternatives in the description of the nuclear modulus: either exponential functions with exponents linear in t (called “purely exponential”) or exponential functions with higher-degree polynomials of t in the exponent (called “non-exponential”).

6.1 Theoretical framework

The amplitude describing elastic scattering of protons may be expected to receive three contributions, each corresponding to one of the following sets of Feynman diagrams.

  • Containing QED elements only. This amplitude can be obtained by perturbative calculations, see Sect. 6.1.1.

  • Containing QCD elements only. This amplitude is not directly calculable from the QCD lagrangian, Sects. 6.1.2 and 6.1.3 will propose several phenomenologically motived parametrisations.

  • Containing both QED and QCD elements. This contribution can neither be directly calculated from the Lagrangians, nor can ad hoc parametrisations be used – this amplitude is correlated with the previous two. Section 6.1.4 will introduce several interference formulae attempting to calculate the corresponding effects.

6.1.1 Coulomb amplitude

The Coulomb amplitude can be calculated from QED (e.g. [28]), using empirical electric \(\mathcal{F}_\mathrm{E}\) and magnetic \(\mathcal{F}_\mathrm{M}\) form factors of the proton. It can be shown (e.g. Section 1.3.1 in [29]) that, at low |t|, the effect of both form factors can be described by a single function \(\mathcal{F}\):

$$\begin{aligned} {\mathrm{d}\sigma ^\mathrm{C}\over \mathrm{d}t} = {4\pi \alpha ^2\over t^2}\,\mathcal{F}^4 ,\ \mathcal{F}^2 = {\mathcal{F}_\mathrm{E}^2 + \tau \mathcal{F}_\mathrm{M}^2\over 1 + \tau } ,\ \tau = {|t|\over 4m^2} , \end{aligned}$$
(16)

where \(\alpha \) is the fine-structure constant and m represents the proton mass.

Fig. 14
figure 14

Illustration of nuclear-phase forms. The dashed lines correspond to predictions by theoretical models ([31] and references therein). The solid lines are typical examples of phases used in this report, all with the same value of \(\rho = 0.10\). The peripheral example corresponds to the parameter values in Eq. (22)

6.1.2 Nuclear amplitude—modulus

At \(|t| \gtrsim 0.02\,\mathrm{GeV^2}\) the effects due to the Coulomb interaction are not expected to be large (c.f. Fig. 15 or [30]). Thus, the measured cross-section can be attributed – to a large extent – to the nuclear component. Following Table 3 and our previous publication [14] with high-precision data for \(|t| < 0.2\,\mathrm{GeV^2}\), the nuclear modulus will be parametrised as

$$\begin{aligned} \left| \mathcal{A}^\mathrm{N}(t) \right| = \sqrt{s\over \pi } {p\over \hbar c} \sqrt{a} \exp \left( {1\over 2} \sum \limits _{n = 1}^{N_b} b_n\, t^n \right) , \end{aligned}$$
(17)

where \(N_b\) is the number of free parameters in the exponent. Consistently with [14]Footnote 1, the parameter \(b_1\) gives the forward diffractive slope and a the intercept of the differential cross-section at \(t=0\). This parametrisation is also compatible with a number of theoretical models (see e.g. [31]).

Since the calculation of CNI may, in principle, involve integrations (e.g. Eq. (25)), it is necessary to extend the nuclear amplitude meaningfully to \(|t| > 0.2\,\mathrm{GeV^2}\). Therefore the parametrisation Eq. (17) is only used for \(|t| < 0.2\,\mathrm{GeV^2}\) while at \(|t| > 0.5\,\mathrm{GeV^2}\) the amplitude is fixed to follow a preliminary cross-section derived from the same data set as in [14] which features a dip-bump structure similar to the one observed at \(\sqrt{s} = 7\,\mathrm{TeV}\) [32]. In order to avoid numerical problems, the intermediate region \(0.2< |t| < 0.5\,\mathrm{GeV^2}\) is modelled with a continuous and smooth interpolation between the low and high-|t| parts. It will be shown that altering the extended part of the nuclear amplitude (\(|t| > 0.2\,\mathrm{GeV^2}\)) within reasonable limits has negligible impact on the results presented later on.

6.1.3 Nuclear amplitude—phase

The following phase parametrisations are considered.

  1. (a)

    A constant phase is obviously the simplest choice:

    $$\begin{aligned} \arg \mathcal{A}^\mathrm{N}(t) = {\pi \over 2} - \arctan \rho = \hbox {const}. \end{aligned}$$
    (18)

    It leads to a strict proportionality between the real and the imaginary part of the amplitude at all t.

  2. (b)

    The standard phase parametrisation,

    $$\begin{aligned} \arg \mathcal{A}^\mathrm{N}(t) =&{\pi \over 2} - \arctan \rho + \arctan \left( \frac{|t|-|t_{0}|}{\tau }\right) \\&- \arctan \left( \frac{-|t_{0}|}{\tau }\right) \, , \end{aligned}$$
    (19)

    describes the main features of many theoretical models – almost imaginary amplitude in the forward direction (\(\rho \) small) while almost purely real at the diffraction dip. The parameter values \(t_0 = - 0.50\,\mathrm{GeV^2}\) and \(\tau = 0.1\,\mathrm{GeV^2}\) have been chosen such that the shape is similar to a number of model predictions, see Fig. 14.

  3. (c)

    The parametrisation by Bailly et al. [33]:

    $$\begin{aligned} \arg \mathcal{A}^\mathrm{N}(t) = {\pi \over 2} - \arctan {\rho \over 1 - {t\over t_\mathrm{d}}} \end{aligned}$$
    (20)

    where \(t_\mathrm{d} \approx -0.53\,\mathrm{GeV^2}\) gives the position of the diffractive minimum at \(8\,\mathrm{TeV}\) (preliminary result derived from the \(\beta ^* = 90\,\mathrm{m}\) data [14]). This phase has a behaviour qualitatively similar to the model of Jenkovszky et al., see Fig. 14.

  4. (d)

    Another parametrisation was proposed in [16]:

    $$\begin{aligned} \arg \mathcal{A}^\mathrm{N}(t) = {\pi \over 2} - \arctan \rho - \zeta _1 \left( - {t\over 1\,\mathrm{GeV^2}} \right) ^\kappa \mathrm{e}^{\nu t}.\nonumber \\ \end{aligned}$$
    (21)

    As shown in Fig. 14, it features a peak at \(t = -\kappa / \nu \) and for asymptotically increasing |t| it returns to its value at \(t=0\). Due to a potentially rapid variation at low |t|, this functional form can yield an impact-parameter-space behaviour that is qualitatively different from the one obtained with the above parametrisations. In order to ensure fit stability, the parameters

    $$\begin{aligned} \zeta _1 = 800 ,\quad \kappa = 2.311 ,\quad \nu = 8.161\,\mathrm{GeV^{-2}} \end{aligned}$$
    (22)

    have been fixed to example values maintaining the desired impact-parameter behaviour at \(\sqrt{s} = 8\,\mathrm{TeV}\), using a method detailed in [34]. This parametrisation with one free parameter will be denoted as peripheral phase in what follows.

Figure 14 shows on the same plot a comparison of phase predictions by several models to typical examples of parametrisations proposed above.

It should be noted that the nuclear phase has a strong influence on the amplitude behaviour in the space of impact parameter b (for a detailed discussion see e.g. Section 3 in [35]). A particularly decisive feature is the rate of phase variation at low |t|. Looking at Fig. 14 one can see that the constant, standard and Bailly phases are essentially flat at low |t|, thus leading to qualitatively similar pictures in the impact parameter space: elastic collisions being more central (preferring lower values of b) than the inelastic ones. Conversely, the peripheral phase parametrisation can yield a description with the opposite hierarchy, which is argued to be more natural by some authors (e.g. Section 4 in [36]). An impact-parameter study of the presented data will be given at end of Sect. 6.4.

6.1.4 Coulomb–nuclear interference formulae

The simplified West–Yennie formula (SWY) [12] was derived in the framework of perturbative quantum field theory by evaluating the lowest-order Feynman diagrams that comprise both nuclear and Coulomb interactions. In this approach, the interference is reduced to an additional phase between the Coulomb and nuclear amplitudes. Moreover, several approximations were used in the derivation. First, in order to avoid integrating over off-mass-shell contributions to the nuclear amplitude (essentially unknown), a very slow variation of the nuclear amplitude phase was assumed: \(\arg \mathcal{A}^\mathrm{N} \approx \hbox {const}\). Then, in order to obtain a closed-form expression, the exponential slope of the nuclear modulus was assumed constant (i.e. only the \(b_1\) parameter is non-zero in the parametrisation Eq. (17)) which is formally incompatible with the existence of the diffractive minimum. The original formula did not contain the electromagnetic form factor \(\mathcal{F}\), which was added later by hand:

$$\begin{aligned}&{\mathrm{d}\sigma \over \mathrm{d}t}^\mathrm{C+N} = {\pi (\hbar c)^2 \over s p^2} \left| {\alpha s\over t} \mathcal{F}^2 \mathrm{e}^{\mathrm{i}\alpha \Phi (t)} + \mathcal{A}^\mathrm{N} \right| ^2 ,\nonumber \\&\Phi (t) =- \left( \log {b_1 |t|\over 2} + \gamma \right) , \end{aligned}$$
(23)

where \(\alpha \) is the fine-structure constant and \(\gamma \doteq 0.577\) the Euler constant. Despite the many limitations, the formula has been extensively used in past data analyses. For backward-comparison reasons it is also considered in this report.

Fig. 15
figure 15

Illustration of the effects due to the Coulomb interaction, using the KL formula. With the Cahn formula the plot looks identical. For the SWY formula, the picture is similar, however it misses the effects at \(|t| \gtrsim 0.02\,\mathrm{GeV^2}\). The curves show a response of the interference formula to different nuclear phases with a purely exponential nuclear modulus. The solid curves correspond to phases of the same shape (constant) but different values of \(\rho \): the maximal response can be seen at \(|t| \lesssim 0.01\,\mathrm{GeV^2}\). Conversely, the dashed lines correspond to phases with fixed \(\rho \) but various shapes (the same examples as in Fig. 14): the response may be sizeable (e.g. in the peripheral case) also at \(|t| \gtrsim 0.02\,\mathrm{GeV^2}\)

The approach of Cahn [15] uses an impact parameter formalism and is based on the additivity of eikonals. The first part of his derivation does not impose any limit on the nuclear amplitude, leading to the formula (Eq. (30) in [15]):

$$\begin{aligned} {\mathrm{d}\sigma \over \mathrm{d}t}^\mathrm{C+N}&= {\pi (\hbar c)^2 \over s p^2} \left| -{\alpha s\over q^2} \mathcal{F}^2 + \mathcal{A}^\mathrm{N}\, [1 - \mathrm{i}\alpha G(-q^2)] \right| ^2 ,\nonumber \\ G(-q^2)&= - \int \limits _0^\infty \mathrm{d}q'^2 \log {q'^2\over q^2} {\mathrm{d}\phantom {q'^2}\over \mathrm{d}q'^2} \mathcal{F}^2(-q^2)\nonumber \\&\quad + {1\over \pi } \int \mathrm{d}^2 q' {\mathcal{F}^2(-q'^2)\over q'^2} \left[ {\mathcal{A}^\mathrm{N}\left( -[\mathbf {q} - \mathbf {q'}]^2\right) \over \mathcal{A}^\mathrm{N}(-q^2)} - 1 \right] , \nonumber \\ \end{aligned}$$
(24)

where \(t=-q^2\), \(\mathbf {q}'\) is a two-dimensional vector and \(q'^2 = |\mathbf {q}'|^2\). The second part of the article gives simplified formulae for nuclear amplitudes with purely-exponential modulus and constant phase and is, thus, of limited interest for the present analysis.

Kundrát and Lokajíček (KL) [16] transformed the formula of Cahn, Eq. (24), into a form better suited for practical applications and added the kinematic limits on the momentum transfer:Footnote 2

$$\begin{aligned} {\mathrm{d}\sigma \over \mathrm{d}t}^\mathrm{C+N}&= {\pi (\hbar c)^2 \over s p^2} \left| {\alpha s\over t} \mathcal{F}^2 + \mathcal{A}^\mathrm{N}\, [1 - \mathrm{i}\alpha G(t)] \right| ^2 ,\nonumber \\ G(t)&= \int \limits _{-4p^2}^0 \mathrm{d}t'\, \log {t'\over t} {\mathrm{d}\phantom {t'}\over \mathrm{d}t'} \mathcal{F}^2(t')\nonumber \\&\quad - \int \limits _{-4p^2}^0 \mathrm{d}t' \left( {\mathcal{A}^\mathrm{N}(t') \over \mathcal{A}^\mathrm{N}(t)} - 1 \right) { I(t, t')\over 2\pi } , \nonumber \\ I(t, t')&= \int _0^{2\pi } \mathrm{d}\phi \ {\mathcal{F}^2(t'')\over t''} ,\nonumber \\ t''&= t + t' + 2\sqrt{t\, t'} \cos \phi . \end{aligned}$$
(25)

A slightly different variant proposed in Eq. (22) in [37] was considered, too:

$$\begin{aligned} {\mathrm{d}\sigma \over \mathrm{d}t}^\mathrm{C+N} = {\pi (\hbar c)^2 \over s p^2} \left| {\alpha s\over t} \mathcal{F}^2 + \mathcal{A}^\mathrm{N}\, \mathrm{e}^{- \mathrm{i}\alpha G(t)} \right| ^2. \end{aligned}$$
(26)

The interference formula by Cahn, Eq. (24), and the KL formula, Eq. (25), are very similar by construction and therefore they give practically identical interference effects.

Since the quantities G in Eqs. (24) and (25) are complex, the interference effects in these treatments are generally more feature-rich than with the SWY formula, Eq. (23), where the interference is reduced to a single additional phase \(\Phi \).

By analysing Eqs. (23), (24) and (25), one can conclude that in the region where the nuclear amplitude dominates (\(|t| \gtrsim 0.003\,\mathrm{GeV^2}\)), the effects due to the Coulomb interaction are of the order of \(\alpha \) or the ratio \(|\mathcal{A}^\mathrm{C}| / |\mathcal{A}^\mathrm{N}|\). In both cases, the magnitude of the interference effects can be expected at a percent level, as shown in Fig. 15. The figure also shows that the effects at different |t| probe different parts of the nuclear phase: maximum sensitivity to \(\rho \) lies at very low |t| while at higher |t| the effects are sensitive to phase values at slightly higher |t|. It can also be observed that for the constant, standard and Bailly phase the effects are very similar and rather mild at higher |t|. This can be understood from a very limited variation of the phase at low |t|, which is the region contributing most to the integral in Eq. (24) or (25). On the contrary, the higher |t| response to peripheral phases can have various forms, often similar to the deviation of the reconstructed cross-section from pure-exponential, see the top plots in Figs. 16 and 17.

Fig. 16
figure 16

Visualisation of the fit results from Table 4 obtained with \(N_b=1\). The continuous (dashed) lines correspond to fits with Cahn or KL (SWY) formula. Note that the fits with constant nuclear phase largely overlap. Top: fits compared to differential cross-section data in a relative reference frame, see the vertical axis label. The reference is identical to the one in [14]. Bottom left: t-dependence of the nuclear phase as extracted from the fits. Bottom right: the effects induced by the Coulomb interaction for each of the fits

Fig. 17
figure 17

Visualisation of the fit results from Table 5 obtained with Cahn or KL formula and \(N_b=3\). The solid lines correspond to fits with different nuclear phases. Top fits compared to differential cross-section data in a relative reference frame, see the vertical axis label. The reference is identical to the one in [14]. Bottom left t-dependence of the nuclear phase as extracted from the fits. Bottom right the effects induced by the Coulomb interaction for each of the fits

6.2 Analysis procedure

In addition to using the data from Table 3, one might consider including the \(\beta ^* = 90\,\mathrm{m}\) data [14] which benefit from much smaller uncertainties. However, due to the limited reach, \(|t| \gtrsim 0.03\,\mathrm{GeV^2}\), they have essentially no sensitivity to the \(\rho \) parameter, cf. Fig. 15. Furthermore, due to possible systematic tensions between the data sets, the inclusion of the \(\beta ^* = 90\,\mathrm{m}\) data may have a deteriorating impact on the \(\rho \) determination. Therefore, the value of \(\rho \) was determined from the \(\beta ^* = 1000\,\mathrm{m}\) data only. For other parameters to which both data sets have non-negligible sensitivity (e.g. \(b_i\) in Eq. (17)), both data sets should give compatible results. This was verified for all the fits that will be presented later on. Since the \(\beta ^* = 90\,\mathrm{m}\) data yield much lower uncertainties, both data sets have been used for determining all parameters except \(\rho \). In practice, a series of two fits was performed:

  • Step 1: fit of \(\beta ^* = 1000\,\mathrm{m}\) data with \(\rho \) free,

  • Step 2: fit of \(\beta ^* = 1000\) and \(90\,\mathrm{m}\) data with \(\rho \) fixed from the preceding step.

The standard least-squares method was used for all the fits. In particular, minimising

$$\begin{aligned} \chi ^2&= \Delta ^\mathrm{T}\mathsf{V}^{-1} \Delta ,\quad \Delta _i = \left. {\mathrm{d}\sigma \over \mathrm{d}t}\right| _{\mathrm{bin}\ i} - {\mathrm{d}\sigma ^\mathrm{C+N}\over \mathrm{d}t}\left( t^\mathrm{rep}_{\mathrm{bin}\ i}\right) ,\nonumber \\ \mathsf{V}&= \mathsf{V}_\mathrm{stat} + \mathsf{V}_\mathrm{syst} , \end{aligned}$$
(27)

where \(\Delta \) is a vector of differences between the differential cross-section data and a fit function \(\mathrm{d}\sigma ^\mathrm{C+N}/\mathrm{d}t\) evaluated at the representative point \(t^\mathrm{rep}\) of each bin [27]. The minimisation is repeated several times, and the representative points are updated between iterations. The CNI effects are calculated using the computer code from [31]. The covariance matrix \(\texttt {V}\) has two components. The diagonal of \(\texttt {V}_\mathrm{stat}\) contains the statistical uncertainty squared from Table 3 and from the Table 3 in [14]. \(\texttt {V}_\mathrm{syst}\) includes all systematic uncertainty contributions except the normalisation, see Eq. (15) and Eq. (14) in [14]. For improved fit stability, the normalisation uncertainty is not included in the \(\chi ^2\) definition. Instead, the uncertainty is propagated for each fit parameter. For this purpose, the fit is repeated with \(-1\,\mathrm{\sigma }\), \(0\,\mathrm{\sigma }\) and \(+1\,\mathrm{\sigma }\) biases independently in: global normalisation (\(1\,\mathrm{\sigma } = 4.2\,\mathrm{\%}\)), \(90\,\mathrm{m}\) data normalisation (\(0.08\,\mathrm{\%}\)) and \(1000\,\mathrm{m}\) data normalisation (\(0.25\,\mathrm{\%}\)). This gives a sample of 27 fit results, from which one can estimate the propagated normalisation uncertainty of a parameter as \((\hbox {max} - \hbox {min})/2\), where “max” (“min”) is the greatest (smallest) value in the sample. This normalisation uncertainty is, at the end, added quadratically to the uncertainty reported by the fit with no bias.

The fits have shown low sensitivity to several of the choices presented above, summarised in the following list.

  • Choice of the form factor in Eq. (16). The options considered in [31] have been tested, none of them giving any significant difference with respect to the default choice [38].

  • Extension of the modulus of the nuclear amplitude to the unobserved |t| region, see the last paragraph in Sect. 6.1.2. No effect was observed when the high-|t| part was altered (both shape and normalisation) nor when the size of the transition region was changed.

  • Use of the Cahn or KL formula. Only the latter will be used in what follows to represent both of them.

  • The two variants of the KL formula, Eqs. (25) and (26). The latter will be used below.

  • Fits with constant, standard and Bailly phase are practically indistinguishable. This can be expected from Fig. 15 showing that the corresponding CNI effects are very similar. Therefore, in the remainder of this article, these phases will be treated as a single family represented by the constant phase.

One of the goals of this study is to probe the origin of the differential cross-section non-exponentiality reported earlier [14]. Therefore, the following two classes of fits were considered.

  • Section 6.3: fits with purely exponential nuclear modulus, that is \(N_b=1\) in Eq. (17). In this case, the non-exponentiality can come from the CNI effects only.

  • Section 6.4: fits with nuclear modulus flexible enough to describe the non-exponentiality without the CNI effects. Here, the non-exponentiality may be due to the nuclear modulus, CNI effects or both.

For each of these nuclear modulus cases, the following two phase parametrisations were considered:

  • constant phase, Eq. (18), as a representative of the central-phases family,

  • peripheral phase, Eq. (21) with parameters fixed to the values in Eq. (22) to represent peripheral behaviour in the impact parameter space.

In each case, the fit results are used to calculate the total cross-section via the optical theorem:

$$\begin{aligned} \sigma _\mathrm{tot}^2 = {16\pi \, (\hbar c)^2\over 1 + \rho ^2}\, a. \end{aligned}$$
(28)

Note that unlike all previous total cross-section determinations at LHC, in this article all the ingredients come consistently from a single analysis.

6.3 Fits with purely exponential nuclear modulus

The goal of this section is to test whether the data are compatible with a purely exponential nuclear modulus, i.e. \(N_b=1\) in Eq. (17). In other words, the non-exponentiality is forced to originate from the Coulomb-induced effects. The fit results obtained with the KL and (where applicable) SWY formulae are summarised in Table 4 and graphically shown in Fig. 16.

Table 4 Fit results with \(N_b=1\). Each column corresponds to a fit with different interference formula and/or nuclear phase

Table 4 shows that both fits with constant phase are essentially identical and have bad quality. The step-2 fit using both \(\beta ^*=1000\) and \(90\,\mathrm{m}\) data can be excluded with \(7.6\,\mathrm{\sigma }\) significance. Consequently, since the combination of \(N_b=1\) and constant phase is the only one compatible with the SWY approach, that formula is experimentally excluded even on the basis of only the low-|t| data set discussed here. This result is complementary to the observation of a diffractive minimum at \(\sqrt{s} = 8\,\mathrm{TeV}\) (to be published in a forthcoming article) which also contradicts the assumptions of the SWY formula.

Although the quality of the fit with the peripheral phase is good, this option seems disfavoured from different perspectives.

  • There are several theoretical reasons for the nuclear component not to be purely exponential, e.g. [3942]. Indeed, most elastic scattering models predict a non-exponential nuclear modulus, see e.g. [31] and references therein.

  • The value of \(\rho \) obtained in this fit may be regarded as an outlier with respect to a consistent pattern of other fits from this article and extrapolations from lower energies: e.g. [4345] and most models in [31].

Let us also recall that the good quality of this fit is possible due to the more complex KL formula where the CNI effects go beyond a simple additional phase in the traditional SWY concept.

6.4 Fits with non-exponential nuclear modulus

The aim of this section is to discuss fits with enough flexibility in the nuclear modulus to describe the non-exponentiality in the data. Since a non-exponential hadronic modulus is used, the only applicable interference formula is KL. \(N_b=2\) to 5 were considered. The optimal degree was chosen according to two criteria: reasonable \(\chi ^2/\hbox {ndf}\) and stability of fit parameters (among which \(\rho \) is one of the most sensitive). For instance, with constant phase the fit (step 1) with \(N_b=2\) yields \(\chi ^2/\hbox {ndf} = 1.07\) and \(\rho = 0.10\) while the one with \(N_b=3\) gives \(\chi ^2/\hbox {ndf} = 1.03\) and \(\rho = 0.12\). Both fits have the normalised \(\chi ^2\) reasonably close to 1, but the value of \(\rho \) changes significantly between \(N_b=2\) and 3 which is unexpected should \(N_b=2\) be sufficient. On the other hand \(N_b=4\) gives \(\chi ^2/\hbox {ndf} = 0.861\) which is unreasonably low. Therefore \(N_b=3\) was chosen.

As shown in Table 5, both fits have reasonable fit quality and remarkably consistent values of \(\rho \) (identical within the resolution) which are compared to previous determinations at lower energies in Fig. 18. Take note that the obtained parameters for the nuclear amplitude (a and \(b_i\)) are consistent between step 1 (\(\beta ^* = 1000\,\mathrm{m}\) data only) and step 2 (both \(\beta ^* = 1000\) and \(90\,\mathrm{m}\) data) of the fitting procedure as already mentioned in Sect. 6.2.

Fig. 17 shows that the level of Coulomb-induced effects is very different in the fits. It is much stronger in the case of the peripheral-phase, which can be expected as this phase features a faster variation in the low-|t| region.

Table 5 Fit results with Cahn or KL formula and \(N_b=3\)
Fig. 18
figure 18

Energy dependence of the \(\rho \) parameter. The blue (green) triangles correspond to \(\mathrm{pp}\) (\({\bar{\mathrm{p}}}\mathrm{p}\)) data from PDG [46] – note that most of these points were determined with the help of the SWY formula, shown to be inconsistent with the present data. The hollow red circle stands for the earlier indirect determination by TOTEM [21]. The filled red circle represents the two results from Table 5 which are numerically identical within the resolution. The black curve gives the preferred \(\mathrm pp\) model by COMPETE [45], obtained without using LHC data

The total cross-section results from the two fits in Table 5 are well consistent with each other and also with previous measurements [14, 22]. The slightly higher values with respect to previous analyses neglecting the Coulomb interaction are expected as long as \(\rho > 0\). This gives negative interference at low |t| and when separated leads to an increase of nuclear cross-section intercept a and thus also total cross-section via Eq. (28).

It is interesting to study the fit behaviour in the impact-parameter space. The scattering amplitude in this representation (sometimes called profile function), \(\mathcal{P}(b)\), can be obtained from the nuclear amplitude by means of Fourier-Bessel transformation (see e.g. [35]):

$$\begin{aligned}&\mathcal{P}(b) = {1\over 4 p \sqrt{s}} \int \limits _{-\infty }^0 \mathrm{d}t\,J_0\left( b\sqrt{-t}\over \hbar c\right) \,\mathcal{A}^\mathrm{N}(t) ,\nonumber \\&\hbox {normalised that } {\sigma _\mathrm{el} = 8\pi \int \limits _0^{+\infty } b\, \mathrm{d}b\, |\mathcal{P}(b)|^2} , \end{aligned}$$
(29)

where \(\sigma _\mathrm{el}\) is the integrated elastic cross-section. The profile functions for the two fits from Table 5 are shown in Fig. 19. The fit with constant nuclear phase gives a distribution peaked at \(b=0\). It corresponds to a behaviour that is more central than for the fit with peripheral phase, where the amplitude modulus reaches maximum at \(b\approx 1.2\,\mathrm{fm}\). These considerations can be extended to inelastic channels. Following Section 3 in [35], one can calculate the mean values of \(b^2\) for elastic (\(\langle b^2\rangle _\mathrm{el}\)), inelastic (\(\langle b^2\rangle _\mathrm{inel}\)) or all (\(\langle b^2\rangle _\mathrm{tot}\)) collisions:

$$\begin{aligned} \langle b^2\rangle _j&= { \int b\,\mathrm{d}b\,b^2\, h_j(b) \over \int b\,\mathrm{d}b\, h_j(b) } ,\qquad h_\mathrm{el}(b) = |\mathcal{P}(b)|^2\nonumber \\ h_\mathrm{tot}(b)&= \mathfrak {I}\mathcal{P}(b) ,\qquad h_\mathrm{inel}(b) = h_\mathrm{tot}(b) - h_\mathrm{el}(b). \end{aligned}$$
(30)

Their values reproduced in Fig. 19 indicate that the fit with constant nuclear phase leads to a picture with elastic collisions more central than the inelastic ones. The hierarchy is inverted for the fit with peripheral phase.

Fig. 19
figure 19

Square of the impact-parameter amplitude, \(\mathcal{P}\), as a function of impact parameter, b. The two lines correspond to the fits in Table 5, using the same colour code as in Fig. 17. The root-mean-squares of b in the legend are calculated from Eq. (30)

7 Summary and outlook

For the first time at LHC the differential cross-section of elastic proton–proton scattering has been measured at |t|-values down to the Coulomb–nuclear interference (CNI) region. This was made possible by a special beam optics, a novel collimation procedure and by moving the RPs to an unprecedented distance of only \(3\,\sigma \) from the centre of the circulating beam.

To fit \(\mathrm{d}\sigma /\mathrm{d}t\) in the CNI region, several interference formulae – Simplified West and Yennie (SWY), Cahn and Kundrát-Lokajíček (KL) – were explored in conjunction with different mathematical descriptions of the modulus and phase of the nuclear amplitude as a function of t. The nuclear modulus was parametrised as an exponential function with a polynomial of degree \(N_b=1\) or 3 in the exponent. These two alternatives allowed to test whether the nuclear modulus can be purely exponential or more flexibility is required. For the phase two options were considered, leading to different impact-parameter distributions of elastic scattering events: a constant phase implying a central behaviour, and another description favouring peripheral collisions. The following conclusions can be drawn.

  • Purely exponential nuclear modulus (\(N_b=1\)), constant phase: excluded with more than \(7\,\mathrm{\sigma }\) confidence. Since this is the only combination compatible with the SWY formula, the data exclude the usage of the formula.

  • Purely exponential nuclear modulus (\(N_b=1\)), peripheral phase: the data do not exclude this option which, however, is disfavoured from other perspectives.

  • Non-exponential nuclear modulus (\(N_b=3\)): both constant and peripheral phases are well compatible with the data, therefore the central impact-parameter picture prevalent in phenomenological descriptions is not a necessity.

The \(\rho \) parameter was for the first time at LHC extracted via the Coulomb–nuclear interference. In the preferred fits (\(N_b=3\)):

$$\begin{aligned} \rho = 0.12 \pm 0.03. \end{aligned}$$
(31)

The new total cross-section determination is conceptually more accurate than in all previous LHC publications since the CNI effects are explicitly treated. Moreover, the value of \(\rho \) comes from the same analysis, not from an external source, which underlines consistency. The \(\sigma _\mathrm{tot}\) values are very well consistent among all non-excluded fits and compatible with the previous measurements. As expected, the new determination yields slightly greater values relative to previous results where the negative CNI was not taken into account. Also note that if the SWY formula with purely exponential hadronic modulus is used, the total cross-section is underestimated by about \(1\,\mathrm{mb}\). A similar underestimation may occur if the non-exponentiality is not taken into account [47].

For even stronger results in the future the key point is a better distinction between the nuclear and CNI cross-section components, which can be achieved from both theoretical and experimental sides. New theory developments may narrow down the range of allowed parametrisations of the nuclear modulus and phase or better constrain the induced CNI effects. The experimental improvements include increasing statistics and reducing the lower |t| threshold. For the former, TOTEM has already upgraded the RP mechanics such that both vertical pots can be simultaneously placed very close to the beam. For the latter, TOTEM foresees an optics with extremely high \(\beta ^* \approx 2500\,\mathrm{m}\) which would allow to reach the CNI region even at Run II energies. Moreover, recent experience with the \(\beta ^* = 90\,\mathrm{m}\) optics at \(\sqrt{s} = 13\,\mathrm{TeV}\) shows that very low beam emittances can be achieved, thus possibly further reducing the RP distance from the beam.