A new method for obtaining a Born cross section using visible cross section data from $e^+e^-$ colliders

In this article, we propose a new method for obtaining a Born cross section using visible cross section data. It is assumed that the initial state radiation is taken into account in a visible cross section, while in a Born cross section this effect is ommited. Since the equation that connects Born and visible cross sections is an integral equation of the first kind, the problem of finding its numerical solution is ill-posed. Various regularization-based approaches are often used to solve ill-posed problems, since direct methods usually do not lead to an acceptable result. However, in this article it is shown that a direct method can be successfully used to numerically solve the considered equation under the condition of a small beam energy spread and uncertainty. This naive method is based on finding a numerical solution to the integral equation by reducing it to a system of linear equations. The naive method works well because the kernel of the integral operator is a rapidly decreasing function of the variable $x$. This property of the kernel leads to the fact that the condition number of the matrix of the system of linear equations is of the order of unity, which makes it possible to neglect the ill-posedness of the problem when the above condition is satisfied. The advantages of the naive method are its model independence and the possibility of obtaining the covariance matrix of a Born cross section in a simple way.


Relationship between Born and visible cross sections
The precise measurement of the inclusive cross section of an e + e − annihilation to hadrons is one of the goals of experiments carried out with electron-positron colliders. This cross section is of interest in connection with the measurement of the fine structure constant and the problem of the muon anomalous magnetic moment. The inclusive cross section of an e + e − annihilation to hadrons is considered as the sum of the cross sections for exclusive processes of an e + e − annihilation into different hadronic states. In such experiments, a visible cross section σ vis is measured, while a Born cross section σ Born is of interest. In this paper, we assume that the initial state radiation is taken into account in a visible cross section, while in a Born cross section this effect is omitted. Visible and Born cross sections are related by the following integral equation: where ε(x, s) is the detection efficiency, that depends on center-of-mass energy (c.m. energy) √ s and the energy fraction carried away by initial state radiation. The upper limit of the integral in eq. (1.1) is determined by the experimental conditions. The visible cross section of a certain process can be found as the ratio of the event yield N (s) of this process to the integral luminosity L int (s): (1. 2) The kernel function F (x, s) from eq. (1.1) has the following form: F (x, s) = βx β−1 1 + α π where α is the fine structure constant, m is the electron mass, L = ln s m 2 and β = 2α π (L−1). The function F (x, s) has an integrable singularity at x = 0. The dependence of this function on x is shown in figure 1 at √ s = 1 GeV. The values of the F (x, s) at √ s = 1 GeV are also listed in table 1. The relationship between visible and Born cross sections, as well as the form of the function F (x, s), were first obtained by E.A. Kuraev and V.S. Fadin in the work [1].
Taking into account that a Born cross section is equal to zero at the energies below the threshold, eq. (1.1) can be rewritten as follows: where s T is the square of the threshold energy. The fact that the upper integral limit form eq. (1.1) can be less than 1 − s T /s is taken into account in eq. (1.4) due to detection efficiency ε(x, s). In this work, we assume that the efficiency is always equal to unity, unless otherwise stated. It should also be noted that eq.  Table 1: The dependence of kernel function F (x, s) on x at √ s = 1 GeV. this paper. Further, it will be shown that due to the properties of the kernel, this equation can be solved numerically with good accuracy without using regularization. It should be noted that the beam energy in colliders has a spread. For instance, this spread is of the order of 1 MeV in the case of the VEPP-2000 [2,3]. It should also be noted that the beam energy is measured with some uncertainty. In order to take these effects into account, eq. (1.4) needs to be modified as follows: where σ 2 E (s) is the sum of the squares of the c.m. energy spread and the uncertainty of this energy. The c.m. energy spread is usually much greater than the uncertainty of its measurement. Therefore, below we refer to the parameter σ 2 E (s) simply as the c.m. energy spread. It should be noted that the value of this parameter depends on the experimental conditions, i.e. this parameter can be different for different c.m. energy points.

Brief overview of conventional methods for obtaining a Born cross section
Let us take a quick look at some of the techniques that are commonly used to obtain a Born cross section using visible cross section data. One of the most frequently used methods is the fit to a visible cross section using eq. (1.4) or eq. (1.5), where a Born cross section is specified in a certain model. The main disadvantage of this approach is its model dependence. In addition, the result of this method is a set of parameters with their uncertainties, on which a Born cross section depends in the considered model. However, it is customary to present a Born cross section in the form of points with error bars. In this form, a Born cross section is more convenient to compare with results obtained in other experiments. In order to obtain a Born cross section in the form of points with error bars, it is first necessary to calculate the radiative correction. This can be done using the following equation: , (1.6) where σ model Born and σ model vis are model Born and visible cross sections, respectively. In the case of using the relationship between visible and Born cross sections given in the form of eq. (1.4) or eq. (1.5), the product ε(s)(1 + δ(s)) in eq. (1.6) is calculated as a whole. In eq. (1.6), ε(s) means the average detection efficiency for each c.m. energy point, and δ(s) means the radiative correction. A model visible cross section is determined using the relationship between Born and visible cross sections given by eq. (1.4) or eq. (1.5), while the dependence of a model Born cross section on c.m. energy and model parameters is given. Further, a Born cross section σ Born in the form of points with error bars can be calculated using the following equation: , where k is the index of the c.m. energy point and σ vis (s k ) is the experimental data of a visible cross section. Next, the error propagation formula can be used in order to calculate the statistical uncertainties of a Born cross section at each c.m. energy point. The simplest way to do this is as follows: , (1.8) where ∆σ Born (s k ) is the statistical uncertainty of a Born cross section and ∆σ vis (s k ) is the statistical uncertainty of a visible cross section. However, the statistical uncertainty of a Born cross section found using eq. (1.8) is not correct. The reason is that the values of a Born cross section at different c.m. energy points are related by means of an integral equation. Consequently, the statistical uncertainties of a Born cross section at different c.m. energy points are not independent. To take into account the fact that the uncertainties are not independent, an additional concept of the radiative correction uncertainty is sometimes introduced. It should be noted, however, that the most reliable way to obtain the statistical uncertainties of a Born cross section is to find the corresponding covariance matrix. In the case of the considered method, the covariance matrix can be estimated using the Monte Carlo technique. In this approach, a visible cross section is randomly generated several times at each c.m. energy point according to its statistical uncertainties. Each time a visible cross section is generated, a Born cross section is calculated using the algorithm described above. The values of a Born cross section obtained in this way can be used to estimate the sample covariance matrix of this cross section. Another commonly used method for obtaining a Born cross section using visible cross section data is based on an iterative procedure for finding the radiative correction. At the first step of this procedure, the fit to the Born cross section σ external Born (s) obtained in some of the previous experiments is substituted into the integral from eq. (1.4) or eq. (1.5) as the zero approximation of a Born cross section. Further, the zero approximation of the radiative correction is calculated through the ratio of the integral from eq. (1.4) or eq. (1.5) to the zero approximation of a Born cross section. At the i-th step of the iterative procedure, the radiative correction and the next approximation of a Born cross section are calculated using the following equations: where the visible cross section σ (i) vis at the i-th step is calculated using eq. (1.4) or eq. (1.5) and the i-th approximation of a Born cross section σ (i) Born . As in the case of a first method, the term ε (i) (s)(1 + δ (i) (s)) is calculated as a whole. The main disadvantage of this method is that, as in the case of the first method, the covariance matrix of a Born cross section can be estimated using the Monte Carlo technique, which has a fairly high computational complexity.

Integral equation of the first kind as an example of ill-posed problem
Note that eq. (1.1) is the Fredholm integral equation of the first kind, and equations (1.4) and (1.5) can be written in the form of the Volterra integral equation of the first kind. There is a well-known fact that finding a numerical solution to Fredholm or Volterra integral equations of the first kind is an ill-posed problem. A problem is said to be well-posed in the sense of Hadamard [4,5] if: 1. a solution exists; 2. that solution is unique; 3. the solution is changes continuously with changes in the data.
Otherwise, the problem is ill-posed.
Before discussing integral equations (1.4) and (1.5) in detail, let us consider an arbitrary Fredholm integral equation of the first kind: where f (t) is the given function, q(x) is the unknown function, K(t, x) is the given kernel function, a and b are the integration limits. The fact that the problem given by eq. (1.10) is ill-posed means mainly that any arbitrarily small perturbations of the function f (t) can lead to the significant changes in the solution of the integral equation. This property of the integral equation (1.10) can be explained using the Riemann-Lebesgue lemma [6]. A corresponding explanation can be found in the book [7]. It is assumed that the dependence of the function f (t) on the variable t is measured experimentally, so that the function f (t) is always perturbed due to the presence of uncertainties in its measurement. This means that due to the ill-posedness of the problem given by eq. (1.10), the numerical solution of this equation may differ significantly from the exact function q(x). Linear integral equation can be approximately reduced to a system of linear equations. Let us consider a discrete problem that corresponds to the integral equation (1.10): Despite the fact that other options are possible, for the purposes of this work, it is sufficient to consider the case when theK is a full rank square matrix. In this case, the system of linear equations (1.11) can be naively solved as follows q =K −1 f . Suppose that the right side of eq. (1.11) is perturbed, which leads to a perturbation of the left side. If we neglect the perturbations of the matrixK, it is easy to show that the relative perturbations of q and f are related by means of the following inequality [8]: where ∆f is the perturbation of the right side, ∆q is the perturbation of the numerical solution of eq. (1.11) and cond(K) = K 2 K −1 2 is the condition number of the matrixK.
In a more general case, when the matrixK is also perturbed, the following inequality [9] holds: where ∆K is the perturbation of the matrixK. The inequality (1.13) is obtained under the assumption K −1 2 ∆K 2 < 1. Note that in this paper only the spectral norm of the matrices 1 is used. Almost always, in the case of ill-posed problems, the matrixK is ill-conditioned, i.e. cond(K) 1. Equations (1.12) and (1.13) show that even for small perturbations of the right-hand side of eq. (1.11), the perturbations in the solution of this equation can be large if the matrixK is ill-conditioned.
Another way to understand the nature of perturbations in the solution of eq. (1.11) is to consider the singular value decomposition (SVD) of the matrixK. The SVD decomposition of the matrixK can be written in the following form: (1.14) whereÛ andV are orthogonal matrices, the matrixΣ is a diagonal matrix whose diagonal elements σ 1 , σ 2 , . . . , σ n are always non-negative numbers and called singular values. Let us choose the numbering of singular values in such a way that these values decrease as the index increases. Using the SVD decomposition of the matrixK, the direct solution to eq. (1.11) can be written in the following form: where vectors u i and v i are the columns of the matricesÛ andV, respectively. Eq. (1.15) shows that terms with small singular values in the denominator can make a larger contribution to the solution than terms with large singular values. It is also important to note the well-known fact that the frequency of oscillations of the components of the vector v i increases as the index i increases, i.e. as the corresponding singular values decrease. Thus, with a large condition number cond(K) = σ n /σ 1 , high-frequency perturbations make a large contribution to the solution.
The simplest way to avoid high-frequency perturbations in the solution is to exclude the terms with small singular numbers from the sum in eq. (1.15). This approach to the solution regularization is called the truncated SVD method (TSVD). The solution q TSVD in this case is given by the following formula: The spectral norm K 2 = sup K x 2 = σmax of an arbitrary matrixK ∈ F m×n is induced by the Value σmax is the maximum singular value of the matrixK and field F is the field of real (R) or complex (C) numbers.
where n < n. However, in practice, the Tikhonov regularization method is most often used. This method consists in finding the minimization of a functional of the following form: where λ is a fixed parameter that has the meaning of the strength of the regularization. The functional L[q] is minimized by the components of the vector q. In a more general case, Tikhonov regularization method consists in minimizing a slightly more complex functional given by the following formula: whereB is some matrix. The simplest examples of such a matrix are the identity matrix or a matrix that links the vector q with the vector of numerical derivatives of the function q(x) with respect to its argument. There are a number of semi-heuristic methods for choosing the optimal value of the regularization parameter λ. One of these methods is the L-curve criterion. This criterion is considered in section (4), where the Tikhonov regularization method is discussed. It is also important to note that the minimization of the generalized functional L G [q] can be reduced to the selective SVD method (SSVD). In the SSVD method, the regularized solution q SSVD is calculated according to the formula: where φ i are some weight coefficients. In the simple case, when the matrixB in eq. (1.18) is equal to the identity matrix (i.e. in the case of eq. (1.17)), the weight coefficients can be found using the following formula 2 : (1.20) More details about the methods for solving ill-posed problems, as well as the methods of optimal choice of the regularization parameter, can be found in the works [7, 10-13].

Proposed method
The condition number depends on the properties of a particular integral operator and can indeed be large in many applications such as image processing [14][15][16][17], geophysics [18][19][20][21] and high energy physics [22][23][24][25][26]. However, as it is shown in section 2, the condition number of the matrix of the integral operator from eq. (1.4) is of order unity due to the properties of the function F (x, s). Therefore, one can hope that the ill-posedness of problem given by eq. (1.4) can be neglected. Since the matrix of the system of linear equations corresponding to integral equation (1.4) is well conditioned, the numerical solution can be obtained directly by inverting this matrix. This is the new method for obtaining a Born cross section, proposed in this paper. Further, this method will be referred to as the naive method. The conditionality of the matrices corresponding to integral equations (1.4) and (1.5) are discussed in section 2. Section 3 presents the results of applying the naive method obtained using a number of numerical experiments. Note, that when eq. (1.5) is reduced to a system of linear equations, the condition number can be large if the c.m. energy spread is comparable to or greater than the distance between c.m. energy points. Therefore, the numerical solution to this equation cannot always be found using the naive method, i.e. without using regularization techniques.
On the other hand, it is known [23] that when using various regularization techniques, such as Tikhonov regularization, the covariance matrix of the regularized solution is incorrect, since this solution is biased. For this reason, regularization techniques can hardly be used for precise obtaining of a Born cross section. By the phrase "discretization of the problem" in this work we mean the reduction of an integral equation to a system of linear equations. In order to reduce the integral eq. (1.4) to a system of linear equations, we first interpolate an unknown Born cross section, taking into account the fact that the values of a Born cross section at the c.m. energies below the threshold energy are equal to zero. The next step is to linearly express the interpolation coefficients in terms of the unknown values of a Born cross section at points with the c.m. energies √ s 1 , √ s 2 , . . . , √ s N . Assuming that the measurements of visible cross section and detection efficiency were carried out at points with these c.m. energies, we can write the following system of equations: where σ interp Born (s) is the interpolation function of a Born section, which linearly depends on the unknown values of a Born cross section at the points with the c.m. energies √ s 1 , √ s 2 , . . . , √ s N . After taking the integrals in eq. (2.1), we obtain the following system of linear equations: whereF is the matrix of this system of linear equations, σ Born is the vector of the Born cross section unknown values and σ vis is the vector of the visible cross section values at the considered c.m. energies. To obtain the system of linear equations (2.2) in the above way, it is necessary that the used interpolationL has the linearity property. i.e. the interpolation coefficients should linearly depend on the unknown values of a Born cross section or, which is the same, the following equality should be satisfied: where ψ(x) is the function linearly expressed in terms of the functions h i (x): and c i are arbitrary real coefficients. In this work, piecewise linear interpolation is used, as well as cubic spline interpolation. Both of these interpolation kinds satisfy the linearity property (2.3) with respect to the values of a Born cross section at different c.m. energy points. A package used in this work to obtain numerical solutions to equations (1.4) and (1.5) also includes the ability to alternate these kinds of interpolation at different c.m. energy intervals. This package is called ISRSolver and is discussed in more detail in section 5.
It should be noted that when the piecewise linear interpolation is used, the matrixF is a lower triangular matrix with small off-diagonal elements. An example of the matrixF obtained using piecewise linear interpolation is shown in figure 2. This matrix corresponds to the threshold energy of the e + e − → ηπ + π − process and is obtained using 50 equally spaced c.m. energy points belonging the range from 1.18 to 2.00 GeV. The condition number of the matrixF in this case is approximately equal to 1.6. Let us consider an example similar to the previous one, but corresponding to an extremely non-uniform distribution of c.m. energy points. The positions of the c.m. energy points for this example are shown in figure 3. The number of points, as in the previous case, is 50. Despite the non-uniform distribution of the c.m. energies, the condition number is still comparable to unity. In this case, the condition number is approximately 1.9.
In the case of cubic spline interpolation or mixed type interpolation, the matrixF is not a lower triangular matrix. If the interpolation function describes a Born cross section well, that is, this function does not have oscillatory outliers between the points of this cross section, then the condition number is approximately the same as in the examples described above.

Calculation and estimation of the condition number
The condition number from the previous and subsequent examples is calculated using the formula cond(F) = σ max /σ min . To find the minimum (σ min ) and maximum (σ max ) singular values, the singular value decomposition of the matrixF is performed numerically. In the ISRsolver package, the singular value decomposition of matrices is carried out using the Eigen 3 linear algebra library [27]. Finding an analytical expression for the condition number even in the case of piecewise linear interpolation of a Born cross section is a quite complicated task due to the bulky kernel function F (x, s) and the large size of the matrixF. However, in the case of piecewise linear interpolation of a Born cross section, a rough estimation of the condition number can be obtained analytically, which makes it possible to qualitatively understand the reason that the condition number of the matrixF is of the order of unity. It is known that in the case when the matrixĈ is a normal matrix, i.e. [Ĉ,Ĉ † ] = 0, the condition number can be found as the ratio of the largest eigenvalue to the smallest one. A special case of a normal matrix is a diagonal matrix. Since the off-diagonal elements of the matrixF are small, this matrix is close to a diagonal matrix, and as a consequence to a normal matrix. Based on this assumption, it is possible to obtain a rough estimation of the condition number using the formula cond(F) = λ max /λ min , where λ max is the maximum eigenvalue of the matrixF and λ min is the minimum one. If the matrixF is lower triangular matrix, its eigenvalues coincide with the diagonal matrix elements. Diagonal matrix elements in the case of piecewise linear interpolation are given by the following formula: where the index j takes values from one to N and s 0 = s T . For simplicity, the last equation is obtained under the assumption that the detection efficiency is equal to unity. Further, assuming that the c.m. energies are located equidistantly, it can be easily shown that the diagonal elements of the matrixF decrease while the index j increases. Thus, taking into account the above assumptions, we can roughly estimate the condition number by the following formula: The last formula in the case of the matrix shown in figure 2 gives an estimation of the condition number equal to 1.3, which is 0.3 less than the exact value of the condition number in this case.

Closeness of equation (1.4) to a well-posed problem
It can be shown in an alternative way that the ill-posed problem given by eq. (1.4) is close to some well-posed problem. To begin with, let us introduce a typical c.m. energy scale Γ at which a cross section changes and assume that this scale is greater than or of the order of several MeV. Let us also introduce a condition on the range 0 ÷ ∆x of the argument x, at which a cross section changes weakly: This condition leads us to the following inequality for the upper bound of the corresponding x-range: Taking into account the last inequality, eq. (1.4) can be written as follows: In this case the parameter ∆x must satisfy the following inequality: ∆x 8.4 × 10 −3 . When the parameter ∆x = 2 × 10 −5 is set, the function g(m 2 φ , ∆x) is approximately equal to 0.52, i.e. has value comparable to one. Since division by a number of the order of unity is a well-posed problem [13], let us divide eq. (2.8) by The equation obtained as the result of division can be written as the Volterra integral equation of the second kind. It is well known fact [13] that the Volterra integral equation of the second kind is a well-posed problem. It is also worth noting that the parameter ∆x can be always chosen in a such way for any reasonable parameters Γ and √ s that the coefficient g(s, ∆x) is comparable to one and the inequality (2.7) is satisfied. Thus, the problem given by eq. (1.4) is close to some well-posed problem. (1.5) Equation (1.5) can be reduced to a system of linear equations in the same way as eq. (1.4). The system of linear equations in this case has the following form:

Discretization of equation
where matrixĜ corresponds to the external integral operator from eq. (1.5), while matrix F corresponds to the internal integral operator and in fact is the same as the analogous matrix from eq. (2.2). However, due to the limits in the outer integral from eq. (1.5), the matrixĜF of the full integral operator is not a lower triangular matrix in this case even if piecewise linear interpolation is used. Since the kernel 1 √  MeV. The reason is that the distance (16.4 MeV) between the c.m. energy points in this case is much greater than the parameter σ E . As the density of c.m. energy points increases, the dependence of the condition number on the parameter becomes steeper. An example of this is the dashed curve shown in the same figure. This curve also represents the dependence of the condition number on the c.m. energy spread, but obtained at twice the density of c.m. energy points than the solid curve. Some oscillations are observed on the dashed curve when the energy spread exceeds 6 MeV. Similar oscillations are observed for the solid curve too, but they appear outside the figure starting at about 14 MeV.
In the case when the distance between some c.m. energy points is less than or of the order of the c.m. energy spread, the condition number increases dramatically. A typical dependence of the condition number on the c.m. energy spread in this case is shown in figure 4b. This dependence is obtained with the highly non-uniform distribution of c.m. energy points, which is shown in figure 3. From figures 4b and 3 it follows that if some c.m. energy points are located at a distance much less than the c.m. energy spread from each other, then the condition number grows rapidly with an increase of the c.m. energy spread and can be large even when the c.m. energy spread is of the order of several MeV. In this case, the numerical solution of eq. (1.5) at such points is strongly scattered. If the number of close energy points is small, then a large scatter is observed only at these points.

Covariance matrix
Before considering examples of obtaining a Born cross section using the naive method, let us derive the relationship between the covariance matrices of Born and visible cross sections. In this section, we denote the matrix of the integral operator asÂ. Therefore, in the case of eq. (1.4),Â =F, while in the case of eq. (1.5),Â =ĜF. The covariance matrixΛ of a visible cross section can be written as follows: whereM is the covariance matrix of a Born cross section, the functional E[.] denotes the expectation of a value in square brackets and the functional Cov[.] denotes the covariance matrix of a vector in square brackets, it is also assumed that the summation is taken over the repeated indices. Thus, the covariance matrixM of a Born cross section is expressed in terms of the covariance matrixΛ of a visible cross section and the integral operator matrix A by means of the following formula: In this paper, we assume that both covariance matricesΛ andM correspond only to statistical uncertainties. It is assumed that a visible cross section at different c.m. energy points is measured independently, so that its covariance matrixΛ is a diagonal matrix. In principle it is possible to use a non-diagonal matrixΛ. All the formulas will be the same in this case. In the first example, let us consider the Born cross section of the e + e − → ηπ + π − process within the framework of the vector meson dominance model (VMD). Since we are only interested in the test of the naive method for obtaining the Born cross section, we consider only a simple model containing only ρ → ρη and ρ → ρη intermediate states.
To describe the dependence of the e + e − → ηπ + π − Born cross section on the c.m. energy, the function and parameters from the work [28] Figure 5e shows that there is agreement between the model Born cross section and the numerical solution within a given statistical uncertainty. However, it can be seen from figure 5f that there is a systematic discrepancy of the numerical solution relative to the model Born cross section at the c.m. energies close to the threshold energy. This discrepancy is caused by the fact that interpolation poorly describes the threshold behavior of the cross section. Despite the fact that at the point with the lowest c.m. energy the relative value of the discrepancy reaches 7%, the absolute value of this discrepancy is small due to the smallness of the Born   cross section near the threshold. Figure 5f also shows that the relative discrepancy rapidly decreases with increasing c.m. energy. The discrepancy between the numerical solution and the model Born cross section, as well as the fact that the covariance matrix describes the fluctuations of the numerical solution in various numerical experiments, can be tested by considering a chi-square histogram. Let us consider the chi-square χ 2 model of the numerical solution calculated with respect to the model Born cross section. In this case, the chi-square is given as follows: where σ naive Born is the numerical solution obtained using the naive method, σ model Born is the model Born cross section. The chi-square (3.3) histogram obtained as a result of 10 5 numerical experiments with the cross section of the e + e − → ηπ + π − process is shown in figure 6a. In each numerical experiment, the visible cross section is generated according to the same covariance matrixΛ. The chi-square histogram is fitted with a chi-square probability density function (PDF). The normalization factor (amp. param. in figure 6a) and the number of degrees of freedom (NDF param. in figure 6a) are free fit parameters. The dashed curve in figure 6a corresponds to the fitting function. The fit parameters are also shown in this figure. It can be seen from the figure that the fitting function describes well the chi-square histogram. The value of the NDF parameter obtained as a result of the fit is 50.41, which is slightly different from the expected NDF value. The expected value of NDF is 50 because the cross section is obtained at 50 points. The difference between the NDF parameter in fit and its expected value is associated with the accuracy of the interpolation of the Born cross section. Indeed, if we plot the chi-square distribution with respect to the numerical solution averaged over all numerical experiments, then the deviation of the NDF parameter from the expected value disappears, because the discrepancies associated with interpolation are canceled out in the difference between the numerical solution and its mean value. The corresponding chi-square distribution is shown in figure 6b.
As the second example, let us consider the cross section of the e + e − → π + π − π 0 process. To describe the dependence of the Born cross section of this process on the c.m. energy, the vector-meson dominance model is also used. In this work, to describe the energy dependence of the Born cross section, we use the function and parameters given in the work [29]. However, we also introduced nonzero ω(1420) and ω(1650) contributions in order to make the behavior of the cross section more complicated at energies above the φmeson production energy. The model Born cross section of the e + e − → π + π − π 0 is shown in figure 7a as a solid curve. The considered cross section has two sharp peaks at the ω and φ production energies. In order for interpolation to describe them well, it is necessary to provide a high density of c.m. energy points in the regions of these peaks.    Figure 7b shows that, as in the case with the e + e − → ηπ + π − process, there is a discrepancy between the numerical solution and the model Born cross section near the c.m. energies close to the threshold energy. As noted above, this discrepancy is due to the fact that interpolation poorly describes the threshold behavior of the Born cross section. In addition, a significant relative discrepancy between the numerical solution and the model Born cross section is observed at the c.m. energy of 1.75 GeV. This discrepancy is due to the fact that interpolation does not adequately describe the sharp behavior of the Born cross section at this energy. However, it should be noted that the absolute value of this discrepancy small as well as in the case of the discrepancy near the threshold energy.
The results shown in figure 7 are obtained using mixed interpolation, i.e. piecewise linear interpolation of the Born cross section is used on some c.m. energy ranges, while cubic spline interpolation is used on other ranges. Mixed interpolation is used for the reason that in the energy ranges with abrupt changes in the Born cross section, the accuracy of piecewise linear interpolation decreases significantly. A comparison between the numerical solution and the model Born cross section, similar to that shown in figure 7, is shown in figure 8, but in the case of using piecewise linear interpolation. It can be seen from this figure that in areas with abrupt changes in the cross-section, the interpolation accuracy is indeed significantly lower than in the case of using mixed interpolation.
Similarly, as in the case of the process e + e − → ηπ + π − , the chi-square histograms of the numerical solution for the e + e − → π + π − π 0 is obtained in the cases of mixed interpolation (a) Comparison of the e + e − → π + π − π 0 model Born cross section and the Born cross section obtained using the naive method with the generated visible cross section.  and piecewise linear interpolation. Figure 9 shows the chi-square histograms obtained using mixed interpolation of the Born cross section. Figure 9a shows This discrepancy is associated with a significant discrepancy between the numerical solution and the model Born cross section, caused by the fact that piecewise linear interpolation has low accuracy in the case of sharply varying cross sections. In contrast, the histogram shown in figure 10b (chi-square with respect to averaged numerical solution) corresponds to a chi-square distribution with 200 degrees of freedom.      : Chi-square distribution for the e + e − → π + π − π 0 Born cross section obtained using the naive method. The numerical solution is obtained using mixed interpolation.  Figure 10: Chi-square distribution for the e + e − → π + π − π 0 Born cross section obtained using the naive method. The numerical solution is obtained using piecewise linear interpolation.

Numerical experiments with equation (1.5)
Next, let us consider obtaining a numerical solution to eq. (1.5) using the naive method. As  Figure 11a also shows that error bars of the numerical solution (square roots of the diagonal elements of the corresponding covariance matrix) are also significantly larger than in figure 5e. Since the matrixĜF is a full rank matrix, the covariance matrix of the numerical solution obtained using the naive method is correct [23] in the sense that the square roots of its diagonal elements are correct 68.27% confidence intervals for the Born cross section in the absence of any model knowledge about its behavior. It should also be noted that with an increase in the value of the parameter σ E , the scatter of the points of the numerical solution increases.   We can reduce the scatter of the points of the numerical solution, assuming that it should be more or less smooth. This idea is a core of various regularization methods such as Tikhonov regularization method. However, the assumption of smoothness leads to an additional systematic uncertainty associated with the fact that the numerical solution in this case is biased (regularization error). Figure 11b shows a comparison of the model Born cross section and the regularized numerical solution of eq. (1.5) obtained using the Tikhonov regularization method under the conditions described above. The numerical solution is presented in this figure as points with error bars. The model Born cross section is shown in figure 11b as a solid curve. Figure 11b shows that the points of the numerical solution are close to the curve of the model Born cross section, and their scatter is significantly less than in figure 11a. However, a comparison of figures 11b and 5e shows that the error bars shown in figure 11b are significantly less than the error bars shown in figure 5e. The reason is that since the regularized numerical solution is biased, the square roots of the diagonal elements of its covariance matrix do not represent the correct 68.27% confidence intervals for the Born cross section points. As already noted in section 1.4, this fact significantly limits the applicability of regularization methods for precise obtaining a Born cross section.
In section 2.4, it is shown that if the c.m. energy spread is much less than the distance between c.m. energy points, then the condition number of the matrixĜF is of the order of unity. This fact allows us to hope that in the case of small values of the c.m. energy spread the effects associated with the ill-posedness of the problem given by eq. (1.5) will be small, and the considered problem can be solved using the naive method. Let us consider a numerical experiment similar to the previous one, but with the parameter σ E equal to 2 MeV. The condition number of the matrixĜF in this case is approximately equal to 2.0. A comparison of the model Born cross section and the numerical solution of eq. (1.5), obtained using the naive method, is shown in figure 12a. It can be seen from this figure that the numerical solution and the model Born cross section are in good agreement, and the error bars of the numerical solution are of the same order of magnitude as in figure 5e. Figure 12b shows a comparison of the model Born cross section and the numerical solution obtained using the Tikhonov regularization method. This numerical solution, like the numerical solution discussed in the previous paragraph, corresponds to the parameter σ E equal to 2 MeV. As in the example shown in figure 11b, the error bars of the numerical solution are small compared to the similar error bars shown in figure 5e. It should be noted that the value of the elements of the covariance matrix depends on the choice of the regularization parameter. In the case of the figures 11b and 12b, the value of the regularization parameter is chosen using the L-curve criterion. The L-curve criterion, as well as the Tikhonov regularization method, is considered in detail in section 4.
Let us consider again the naive method for solving eq. (1.5) in the case when the parameter σ E is equal to 2 MeV. As in the examples from section 3.2.1, it is possible to plot the chi-square histogram of the numerical solution obtained using the naive method. This histogram is shown in figure 13. As in the previous examples, the histogram is fitted using the chi-square distribution function. The NDF parameter obtained with fit is 50.55 ± 0.03, which is consistent with the expected number of degrees of freedom equal to 50. However, as in the case of the example of the chi-square distribution from figure 6a, there is a slight systematic deviation of this parameter from the expected number of degrees of freedom, associated mainly with interpolation accuracy. Figure 14 shows the ratio between the numerical solution (σ E = 2 MeV) obtained using the naive method and the model Born cross section, averaged over 10 5 numerical experiments. In each numerical experiment, the visible cross section is generated using the same covariance matrixΛ. It can be seen from the figure that at points with the c.m. energies close to the minimum c.m. energy, there is a significant relative deviation of the numerical solution from the model Born cross section. As in the examples described in section 3.2.1, this deviation is due to the fact that interpolation poorly describes the threshold behavior of the Born cross section. Since the cross section at these c.m. energy points is small, the absolute values of the deviation of the numerical solution from the model Born cross section are also small. Figure 14 also shows that at the point with the maximum c.m. energy, a relative deviation of the numerical solution with respect to the model cross section is observed. This deviation is approximately equal to 0.5% and is associated with the lack of experimental   It should be noted that the naive method works better with smooth cross sections. In the case of sharp cross sections, it is necessary to have a sufficient density of c.m. energy points in order to provide the desired interpolation accuracy. On the other hand, as the density of c.m. energy points increases, the distance between them decreases. At a certain density of c.m. energy points, the distance between them can become of the order of the c.m. energy spread. As discussed above, this leads to a significant scatter in a numerical solution, which is due to the ill-posedness of the problem.
In the example with the process e + e − → π + π − π 0 from section 3.2.1, the density of c.m. energy points is such that the distance between some of these points is comparable to 1 MeV. On the other hand, the authors were unable to provide good interpolation accuracy when using a lower density of c.m. energy points. Thus, the naive method is not applicable for obtaining the cross section of the process e + e − → π + π − π 0 with the c.m. energy spread  In this example, we consider a numerical experiment in which the Born cross sections of the e + e − → ηπ + π − process are found as a numerical solution to eq. (1.4) using the naive method at 50 equally spaced points in the c.m. energy range from 1.18 to 2.00 GeV. Figure 15a shows  Figure 15d shows the ratio of the numerical solution to the model Born cross section, averaged over 10 6 numerical experiments. It can be seen from the figure that the relative deviation at low energies became less than in the similar example with the detection efficiency equal to unity (figure 5f). A significant relative deviation is observed only at the point with the minimum c.m. energy and is approximately equal to 0.7%. At the rest of the energy points, the relative deviation is small. In the case of using the detection efficiency equal to one, the relative deviation in the points with the lowest energy is approximately equal to 7%. The reason for the decrease in the relative difference between the numerical solution and the model Born cross section is that when the function F (x, s) is multiplied by the detection efficiency, which decreases with increasing x, the kernel of the integral operator decreases faster than in the case with the detection efficiency equal to unity.

Tikhonov regularization
In this section, we will briefly consider the application of the Tikhonov regularization method to equations (1.4) and (1.5). The Tikhonov regularization method, unlike the naive method, can be used to find a solution to eq. (1.5) in the case when the parameter σ E is greater or comparable to the distances between c.m. energy points. However, the application of the Tikhonov regularization method is essentially limited by the fact that the diagonal elements of the covariance matrix for the numerical solution of eq. (1.5) (or eq. (1.4)) do not represent the correct 68.27% confidence intervals for the Born cross section points. More details about this fact can be found in work [23].
In this work, the Tikhonov regularization method is applied to eq. (1.5) by means of minimizing the functional of the following form: Comparing equations (4.1) and (1.18), we see that in the second term of the functional L G [σ Born ], the matrixB has been replaced by the matrixD, which in this work denotes the matrix of the c.m. energy derivative operator. The regularization functional L G [σ Born ] given by eq. (4.1) also differs from the standard functional of the Tikhonov regularization by the presence of the chi-square term instead of the term Â σ Born −σ vis 2 . The replacement of the term Â σ Born − σ vis 2 by the chi-square term is performed to identify the minimization of the functional L G with the minimization of the chi-square under the condition λ D σ Born 2 on the smoothness of the numerical solution. In fact, the results obtained using term Â σ Born − σ vis 2 are very similar to those obtained using the chi-square term and can be obtained in a similar manner.  Figure 15: Results of the numerical experiment with the e + e − → ηπ + π − cross section obtained using the naive method with a detection efficiency equal to unity.
but, as noted above, its diagonal elements do not represent 68.27% confidence intervals for the Born cross section at different c.m. energies since the regularized numerical solution σ Tikh Born is biased. Indeed, the expectation of the regularized solution σ Tikh Born can be written in the following form [23]: where the last term is the bias of the regularized solution and the matrixÎ is the identity matrix. It is seen from the last equation that at λ = 0, the bias term disappears.

Numerical experiments and L-curve criterion
Let us consider a numerical experiment with the e + e − → ηπ + π − cross section. We assume that there are 50 equally spaced c.m. energy points in the range from 1.18 to 2.00 GeV. We also assume that the visible cross section is generated at these points using some diagonal covariance matrix in the same way as it is done in section 3. Next, eq. (4.2) is solved using the generated visible cross section for different values of the regularization parameter λ. We also assume that in this numerical experiment the parameter σ E is 20 MeV.  Figure 16 shows that as the regularization parameter increases, the dependence of the numerical solution on the c.m. energy becomes more regular. However, for too large values of the regularization parameter, the numerical solution is systematically lower than the model Born cross section. Therefore, the problem of the optimal choice of the regularization parameter arises in a natural way. There are several semi-heuristic methods to solve this problem. In this work, the L-curve criterion is used. A more detailed description of the methods for choosing the optimal regularization parameter can be found in the works [7,11,23]. Consider a graph in which the values of the first term from functional (4.1) are plotted along the abscissa, and the values of the second term divided by the regularization parameter are plotted along the ordinate. Both the first and second terms depend on the regularization parameter, therefore the coordinates of the points on this graph depend on this parameter. The continuous dependence of these coordinates on the regularization parameter is called an L-curve. An example of an L-curve is shown in figure 17a. There are two typical parts that can be distinguished in the L-curve. The first part of the L-curve is almost vertical and corresponds to small almost constant values of the first term of the functional (4.1). The value of the regularization parameter in this part of the L-curve is insufficient to suppress the effects of ill-posedness of the problem, therefore the value of the term D σ      large values of the first term of the functional (4.1). In this part of the L-curve, the value of the regularization parameter is too large, which leads to a systematically underestimated numerical solution in comparison with the exact Born cross section. Therefore, the optimal regularization parameter corresponds to the region of the L-curve located between the two considered parts. In this region, the curvature of the L-curve reaches a maximum. Thus, the optimal value of the regularization parameter can be estimated as the value of this parameter at which the curvature of the L-curve reaches its maximum. An example of the curvature of the L-curve is shown in figure 17b.
It should be noted that very often the L-curve is brought to a log-log scale in order to emphasize its typical features. The L-curve curvature is also calculated using the L-curve plot at this scale. It was checked that in the case of the considered problem, the L-curve criteria in linear and log-log scales give approximately the same results. However, in the case of the log-log scale, an additional local maximum of the curvature appears near small values of the regularization parameter, which can be falsely interpreted as the maximum curvature at which the optimal value of the regularization parameter is achieved. For this reason, in this work, the L-curve is plotted on a linear scale and its curvature is calculated on the same scale.
The maximum of the L-curve curvature shown in figure 17 is achieved at a regularization parameter of approximately 0.78 GeV 2 nb −2 . A comparison of the numerical solution corresponding to this parameter with the model Born cross section is shown in figure 11b. This figure shows that the regularized numerical solution is in good agreement with the model Born cross section, although the error bars of this solution do not represent the correct 68.27% confidence intervals.

ISRSolver package
Examples of numerical solutions to equations (1.4) and (1.5) considered in this paper were obtained using the ISRSolver package. This package is written using the C++ programming language and also has a Python API for calling some functions. The package includes tools for obtaining a numerical solution using both the naive method and the Tikhonov regularization method. Tools for verifying these numerical solutions are also included. It should be noted that the ISRSolver can be used as a library, for example, to find the convolution of the kernel with an arbitrary function or to implement conventional methods for obtaining a Born cross section. The source code of the package can be found in the repository https://github.com/sergeigribanov/ISRSolver.

Summary
In this work, we considered the problem of finding a numerical solution to eq. (1.4) or eq. (1.5) using visible cross-section data. The problems given by these equations are illposed. However, due to the fact that the condition number of the matrix of the integral operator from eq. (1.4) is comparable to unity, the problem specified by this equation can be solved with a good accuracy using the naive method, i.e. without using regularization. The condition number in this case is of the order of unity due to the fact that the kernel of the integral operator is a function that rapidly decreases with increasing x. The work also shows that the naive method can be used to find a numerical solution to eq. (1.5) at values of the c.m. energy spread that are small compared to the distances between c.m. energy points. Otherwise, the numerical solution found using the naive method has a large scatter due to the ill-posedness of the considered problem. The naive method consists in the direct solving of a system of linear equations that approximately describes the original integral equation. In the case of this work the integral operator matrix is a full rank square matrix. Therefore, according to the work [23], this covariance matrix describes the variability of a Born cross section.
If the c.m. energy spread is greater or comparable to the distances between c.m. energy points, the Tikhonov regularization method can be used to solve eq. (1.5). However, the applicability of this method is limited by the fact that the diagonal elements of the covariance matrix of the regularized numerical solution do not represent 68.27% confidence intervals [23] for a Born cross section at different c.m. energies, since this numerical solution is biased.
The advantages of the naive method are its model independence and the possibility of obtaining the covariance matrix of a Born cross section in a simple way. The main proposal of this work is that the naive method can be used to find a Born cross section using visible cross section data in cases where it is possible, i.e. when there is a sufficient density of c.m. energy points to interpolate a cross section, and the c.m. energy spread is small or not taken into account. Otherwise, some of the conventional methods should be used, for example, a model-dependent fit to a visible cross section using eq. (1.4) or eq. (1.5).