Going to the light front with contour deformations

We explore a new method to calculate the valence light-front wave function of a system of two interacting particles, which is based on contour deformations combined with analytic continuation methods to project the Bethe-Salpeter wave function onto the light front. In this proof-of-concept study, we solve the Bethe-Salpeter equation for a scalar model and find excellent agreement between the light-front wave functions obtained with contour deformations and those obtained with the Nakanishi method frequently employed in the literature. The contour-deformation method is also able to handle extensions of the scalar model that mimic certain features of QCD such as unequal masses and complex singularities. In principle the method is suitable for computing parton distributions on the light front such as PDFs, TMDs and GPDs in the future.

We explore a new method to calculate the valence light-front wave function of a system of two interacting particles, which is based on contour deformations combined with analytic continuation methods to project the Bethe-Salpeter wave function onto the light front. In this proof-of-concept study, we solve the Bethe-Salpeter equation for a scalar model and find excellent agreement between the light-front wave functions obtained with contour deformations and those obtained with the Nakanishi method frequently employed in the literature. The contour-deformation method is also able to handle extensions of the scalar model that mimic certain features of QCD such as unequal masses and complex singularities. In principle the method is suitable for computing parton distributions on the light front such as PDFs, TMDs and GPDs in the future.

I. INTRODUCTION
Understanding the quark-gluon structure of hadrons is a major goal in strong interaction studies. Ongoing and future experiments at the LHC, Jefferson Lab, RHIC, the Electron-Ion Collider, COMPASS/AMBER and other facilities aim to establish a three-dimensional spatial imaging of hadrons and measure structure observables such as the spin and orbital angular momentum distributions inside hadrons and their longitudinal and transverse momentum structure. These properties are encoded in PDFs (parton distribution functions), GPDs (generalized parton distributions) and TMDs (transverse momentum distributions), see e.g. [1][2][3][4][5] and references therein, whose matrix elements are illustrated in Fig. 1. We denoted the field operators generically by Φ(z), T denotes time ordering, O is some operator (usually containing a Wilson line), ∆ is the momentum transfer, and P f,i = P ± ∆/2 are the final and initial momenta of the hadron. A common feature of parton distributions is that they are defined on the light front, i.e., the partons (quarks and gluons) inside the hadron are probed at a lightlike separation. In light-front coordinates this amounts to z + = 0, which by a Fourier transform translates to an integration over q − in momentum space, where q is the relative momentum between the probed partons. In principle, the various quantities of interest can then be derived from Eq. (1) as indicated in Table I: After taking the forward limit ∆ = 0, TMDs follow from an integration over q − and PDFs from another integration over the transverse momentum q ⊥ ; for ∆ = 0 the same steps lead to generalized TMDs and GPDs [4,5]. * gernot.eichmann@tecnico.ulisboa.pt † eduardo.b.ferreira@tecnico.ulisboa.pt ‡ stadler@uevora.pt In practice, taking Eq. (1) to the light front is difficult in Euclidean formulations such as lattice QCD and continuum functional methods. To this end, lattice QCD has seen major recent progress in calculating quasi-PDFs, pseudo-PDFs and other quantities connected to Eq. (1) from the path integral which allow one to reconstruct the desired parton distributions; see e.g. [6][7][8][9][10][11][12][13][14][15][16].
Functional methods, on the other hand, are usually formulated in momentum space. Here the matrix element (1) needs to be constructed in a consistent manner from the elementary n-point correlation functions along the lines of Refs. [17][18][19][20][21][22][23][24][25]. However, the integration over q − is not straightforward unless the analytic structure of the integrands is fully known, which is only the case in simple models where one may employ residue calculus, Feynman parametrizations or similar methods.
In this respect, the Nakanishi integral representation has proven very efficient in recent years [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40]. Here the idea is to recast the hadronic amplitudes that appear in the integrands in terms of a weight function with a denominator that absorbs the analytic structure. There remain however questions regarding the formulation of generalized spectral representations for gauge theories, and the singularity structure of the remaining parts of the integrands (propagators, vertices, etc.) must still be known explicitly which poses practical limitations.
The combination of such techniques, occasionally together with a reconstruction using moments, has found widespread recent applications in the calculation of parton distributions and related quantities [20][21][22][23][24][25].

arXiv:2112.04858v2 [hep-ph] 18 Feb 2022
The goal of the present work is to explore a new technique to compute correlation functions on the light front directly. It is based on contour deformations and analytic continuations, and in principle it does not rely on explicit knowledge of the analytic structure of the integrands except for certain kinematical constraints. Instead of the hadron-to-hadron correlator (1), we consider the simpler case of the vacuum-to-hadron amplitude shown in Fig. 1, This is the generic form of a two-body Bethe-Salpeter wave function (BSWF) for a hadron carrying momentum P , which can be dynamically calculated from its Bethe-Salpeter equation (BSE). Taking this object onto the light front by integrating over q − gives the valence light-front wave function (LFWF), cf. Table I, which will be the central object of interest in this work to establish and test the method. In light-front quantum field theory, the LFWFs are the coefficients of a Fock expansion and thus acquire a probability interpretation [57][58][59][60][61][62][63]. By integrating out also the transverse momentum one obtains the parton distribution amplitude (PDA). In practice we employ a scalar model, namely the massive version of the Wick-Cutkosky model [27,64,65], which encapsulates the relevant features that are also present in more general situations and useful for testing the method. In particular, this model allows for detailed comparisons with the Nakanishi method and its application to LFWFs [32,36], and the results from both methods will turn out to be in excellent agreement. Moreover, the contour-deformation technique is general and not restricted to LFWFs, so it can be applied to hadronto-hadron transition amplitudes as in Eq. (1) and more general theories such as QCD in the future.
The article is organized as follows. In Sec. II we establish the main formalism using Minkowski conventions and work out the LFWF for a simple monopole amplitude. In Sec. III we derive the general expression for the LFWF in a Euclidean metric and analyze the singularity structure of the integrand. In Sec. IV we calculate the LFWF dynamically from its Bethe-Salpeter equation using contour deformations, and we discuss the corresponding results. Sec. V deals with generalizations to unequal masses and complex propagator poles. We conclude in Sec. VI. Two appendices provide details on the Nakanishi representation and the general properties of the singularities that appear in the integrands.

A. Definitions
We begin with some basic definitions. In Minkowski conventions, one may define the light-front components of a four-vector p µ by p ± = p 0 ± p 3 such that A scalar product of two four-vectors then becomes which implies p 2 = p + p − − p 2 ⊥ . For an onshell particle with p 2 = m 2 and p 0 > 0, this entails p + + p − > 0, p + p − = p 2 ⊥ + m 2 > 0 and therefore p ± > 0. The fourmomentum integral in light-front variables reads We now consider the BSWF Ψ(z, P ) defined in Eq. (2) for a bound state of two valence particles. We restrict ourselves to a scalar bound state of two scalar particles, although the generalization to more general theories such as QCD is straightforward. P is the onshell total momentum (P 2 = M 2 ) of the bound state and z is the relative coordinate between the two constituents. The BSWF in momentum space is obtained by taking the Fourier transform, and differs from the BS amplitude Γ(q, P ) by attaching external propagator legs: Here, G 0 (q, P ) is the product of the particle propagators, e.g., for tree-level propagators where q is the relative momentum between the constituents (cf. Fig. 1) and the two particle momenta are given by The variable ε ∈ [−1, 1] is an arbitrary momentum partitioning parameter and equal momentum partitioning corresponds to ε = 0.

3
The light-front wave function (LFWF) ψ(q + , q ⊥ ) is defined as the Fourier transform of the BSWF restricted to z + = 0, which amounts to an integration over q − in momentum space: We used Eqs. (4), (5) and (6), and N is a normalization factor 1 with mass dimension 1. We now introduce a variable α ∈ [−1, 1] through such that the four-momenta in Eq. (9) become The usual notation in terms of the longitudinal momentum fraction ξ ∈ [0, 1], then follows from the identification In the following we set ε = 0 for equally massive constituents (m 1 = m 2 = m), and we work with the variable α instead of ξ to allow for compact formulas where the symmetry α → −α is manifest. Correspondingly, we write the LFWF as The PDA φ(α) and distribution function u(α) then follow from an integration over d 2 k ⊥ : The decay constant f is proportional to the integrated BSWF; integrating φ(α) over α, one finds 1 We introduced N to account for different possible normalizations of the LFWF employed in the literature. Independently of this, the BS amplitude satisfies a canonical normalization condition, which is not important in what follows but implies that Γ in the scalar theory has mass dimension 1. Here we consider Γ to be dimensionless, such that the mass dimension enters through the factor N and Ψ has dimension −4.

B. Light-front wave function for monopole
It is illustrative to work out the LFWF for the case where the BS amplitude is a simple monopole, which is easily evaluated using residue calculus. In the rest frame of the bound state one has where P + = P − = M and thus q + = αM/2. The boundstate mass M must be below the two-particle threshold (M < 2m), and to eliminate the mass parameter m from the equations we define with t ∈ [−1, 0], w ∈ R and x > 0. This entails As a result, the LFWF in Eq. (10) and BSWF in (7-8) take the form where w = w ± are the propagator pole locations corresponding to q 2 1,2 = m 2 − i and w 0 is the pole from the BS amplitude: The residues of the BSWF at the three poles, multiplied with 2πi, are Using Eq. (23) this yields with One may verify that R + + R − + R 0 = 0, i.e., the sum of the residues vanishes as it should. From Eq. (23), the propagator poles at w = w + (w − ) in the complex w plane lie below (above) the real axis. The displacement of the w 0 pole depends on sign(α); for α > 0 it lies below the real axis and for α < 0 above. This is sketched in Fig. 2(a) for α > 0, where the horizontal tracks show the movement of the poles for w ± ± λ and w 0 + λ with a parameter λ > 0. Thus, when integrating w along the real axis and closing the contour at complex infinity, one picks up the residue R − for α > 0 and R − + R 0 = −R + for α < 0. The combined result for the LFWF is therefore which is shown in Fig. 3 for exemplary values of t and γ.
and the distribution function u(α) reads For γ = 1 + t and thus B = 0, and suppressing the prefactors, these quantities reduce to Strictly speaking, with the i factors as in Eq. (23) these results are only valid for real values of α with |α| < 1. For |α| > 1, all poles in the complex w plane in Fig. 2(a) move either to the upper or lower half plane, shifted by i , such that the closed contour is empty and the resulting LFWF is zero. Thus, the LFWF only has support for −1 < α < 1.
On the other hand, the i prescription has its origin in the imaginary-time boundary conditions, which translates to analogous boundary conditions for the q − and w integration: This suggests an integration path like in Fig. 2(b), which starts at a finite imaginary part Im w = −∞ · ε below all poles in the integrand and ends at Im w = +∞ · ε above all poles. For −1 < α < 1 the results from (a) and (b) are identical, but for |α| > 1 or for complex values of w ± and w 0 they are different. In particular, only (b) leads to an analytic function but (a) does not. Substituting |α| → √ α 2 in Eq. (27) yields the result of (b) for any α, x, t, γ ∈ C, which is an analytic function in the complex α 2 plane and plotted in Fig. 4.
By contrast, for complex values of α (with t, x, γ real) Eq. (23) implies For Im α = 0, the i factors become irrelevant and for x + 1 > 0 and x + γ > 0 the three poles always lie in the  (33). To obtain the correct light-front wave function outside R , one would need to deform the Euclidean integration contour. same half plane. Also for Im α = 0 and |α| > 1 the three poles fall in the same half plane. Therefore, the LFWF from (a) vanishes everywhere except for α 2 ∈ R + and is thus not an analytic function.
In the following we adopt the interpretation (b) for the i prescription, since this is what generates an analytic function and will allow us to perform analytic continuations. As long as the integrand vanishes sufficiently fast at complex infinity, one can then equivalently perform a Wick rotation and integrate w along a Euclidean integration path from −i∞ to +i∞.
To facilitate the numerical treatment, instead of working out the poles in the complex w plane we consider the respective pole positions in the complex x plane by solving Eq. (23) for x: When integrating w ∈ (−i∞, i∞), the pole positions turn into branch cuts in x ∈ C as sketched in the left panel of Fig. 5. They are defined by These cuts separate different regions in the complex x plane with different values of the integral. The region corresponding to the proper i prescription is shown in blue and in the following we refer to it as R . A calculation of the LFWF inside this region returns the correct result in Eq. (27).
Vice versa, for the alignment displayed in Fig. 5, x behind the first cut corresponds to the situation where the w 0 pole has moved to the wrong side of the Euclidean integration contour in the complex w plane, thereby not summing the correct residues and giving the wrong value of the integral. If we wanted to know the LFWF for x behind the first cut, we would need to deform the Euclidean contour in w (right panel in Fig. 5). We refrain from doing so in what follows but instead restrict the calculations to the region R .
Below we will transform Eq. (22) to hyperspherical variables, in which case the resulting branch cuts form more complicated curves but still define a corresponding R region. The LFWF can then be calculated numerically inside that region. In turn, this region may not include the whole positive real axis in x, which requires contour deformations in the complex x plane when calculating the light-front distributions in Eq. (16) or when solving for the BSWF Ψ(q, P ).

A. Euclidean conventions
The goal in the following is to transfer Eq. (10) to a Euclidean metric (+, +, +, +). A Euclidean four-vector is defined by In the Euclidean metric the distinction between upper and lower components becomes irrelevant, and scalar products of four-vectors pick up minus signs: The light-front variables in Eq. (3) are independent of the metric, so one has p ± = −ip 4 ± p 3 and The relation (4) turns into We also take the opportunity to remove the various signs and i factors appearing in the Minkowski quantities. In the Euclidean notation, a scalar tree-level prop- We now drop the label 'E' and continue to work with Euclidean conventions.

B. Light-front wave function
To work out the LFWF in Euclidean kinematics, we start from Eq. (15): (38) Because the light-front variables are independent of the metric, the formula has the same form as earlier except that the integration over q − now proceeds from −i∞ to +i∞, and the factor i in the denominator comes from the Euclidean definition of the BSWF Ψ(q, P ) as mentioned above. The latter is given by where the corresponding BS amplitude Γ(q, P ) is the dynamical solution of the BSE which we will discuss in Sec. IV. For the monopole example (18) it is given by We note that for a scalar bound state with scalar constituents, all quantities Ψ(q, P ), Γ(q, P ) and G 0 (q, P ) are Lorentz-invariant. Like in Eq. (11), we define a four-momentum k through q = k + αP/2, where for now we set ε = 0 for equally massive constituents. Because the BSWF is Lorentzinvariant, it can only depend on the Lorentz invariants k 2 , k · P and P 2 = −M 2 together with the momentum partitioning α. We express them through the dimensionless variables x, ω and t: The self-consistent domain of the BSE solution in Sec. IV is x > 0 and ω ∈ [−1, 1], where ω =k ·P is the cosine of a four-dimensional angle (a hat denotes a unit four-vector). The Lorentz invariants q 2 and q · P then become In Euclidean conventions the four-momenta are conveniently expressed in hyperspherical variables, which are more closely related to the Lorentz invariants of the system. In a general moving frame of the total momentum P , this amounts to When the BSE is solved in the moving frame, the natural domain of these variables is x > 0, ψ ∈ [0, 2π) and z, y, Z ∈ [−1, 1]. Because we can always choose P ⊥ = 0, q ⊥ = k ⊥ is automatic. On the other hand, the condition k + = 0 in light-front kinematics entails so that the Lorentz-invariant variable x = k 2 /m 2 = k 2 ⊥ /m 2 = q 2 ⊥ /m 2 assumes the meaning of the squared transverse momentum. Now let us transform the integration over dq − in Eq. (38) to Euclidean variables. From Eq. (37) together with P ⊥ = 0, q + = αP + /2 and P + P − = M 2 , we have (45) and comparison with Eqs. (42) and (20) gives 6. Exemplary cut configurations in the complex √ x plane for α = 0.6, γ = 4 and three different values of √ t. The three lines correspond to the propagator cuts √ x ± (blue, red) and the cut √ x 0 from the monopole amplitude (orange) in Eq. (53). The region R shown in blue connects the origin For given values of Here we also made the dependence of the BSWF on the Lorentz invariants x, ω, t and α explicit. The endpoints of the Euclidean integration contour are ω = ±∞, but the actual contour depends on the singularities of the integrand which we will determine below. What does the constraint k + = 0 mean for the BSWF? Eqs. (43)(44) entail (48) and from P + = −iP 4 + P 3 one finds ω = zP + /M . This does not pose any additional constraints on ω; e.g., for P + > 0 and z ∈ [−1, 1] also the domain ω ∈ R remains the same as before. Then, because the BSWF is frameindependent, we can calculate it in any frame and instead of Eq. (43) we may as well work in the rest frame of the bound state: Formally, the conditions k + = 0 and k 2 ⊥ = m 2 x are not meaningful in this frame, but the Lorentz invariance of the BSWF implies that the result must be identical to that obtained in a moving frame if those conditions are imposed. As a consequence, the final expression for the LFWF becomes Note that this formula no longer makes any reference to light-front variables but only depends on Lorentzinvariant quantities. From Eqs. (9) and (12) one can see that in practice it amounts to evaluating the BSWF for a general momentum partitioning parameter α and subsequently integrating over ω.
The distribution amplitude φ(α) and distribution function u(α) in Eq. (16) then follow accordingly: For general theories and bound states with arbitrary spin, Ψ is not Lorentz invariant but still covariant; thus one can expand it in a Lorentz-covariant tensor basis with Lorentz-invariant coefficients. The operations in light-front kinematics (k + = 0, taking the γ + component, etc.) then affect the tensor basis but still leave the dressing functions invariant. In this way Eq. (50) can also be generalized to matrix elements of the form P f |T Φ(z) O Φ(0)|P i taken on the light front, which enter in the definition of parton distributions such as PDFs, GPDs and TMDs.
In going from Eq. (31) to (50) we effectively changed the integration variable from w to ω through Eq. (46). This also changes the branch cuts in the complex x plane and the corresponding R regions, which are illustrated in Fig. 6 and discussed in detail in the following. This is also the reason that allowed us to drop the constant term ∝ αt from Eq. (46) in the integration path over ω, as this would only lead to different branch cuts and thus a different R region. The difficulty in evaluating Eq. (50) is the singularity structure of Ψ. On one hand, the propagators in Eq. (39) (and generalizations thereof) induce singularities in G 0 ; on the other hand, the solution of the BSE may produce singularities in the BS amplitude Γ.
To determine the singularity structure in the complex √ x plane, we consider again the tree-level propagators in Eq. (39) and the monopole amplitude (40). The momenta of the propagators are given in Eq. (12), which together with (41) entails where λ = ±1 denotes the two solution branches in each case. For real values of ω, Eq. (53) describes branch cuts in the complex √ x plane which are shown in Fig. 6.
They separate eight regions, which relate to the four distinct possibilities of picking up the three pole residues when integrating over ω: picking up none of the residues is equivalent to picking up all three of them and gives zero; picking up one of them is equivalent to picking up the other two with the opposite direction of the integration contour. Without loss of generality, we can restrict ourselves to values of √ t = iM/(2m) in the upper right quadrant. The region R corresponding to the proper i prescription is then the one shown in blue and connects the origin For the monopole example and √ x ∈ R , one can easily verify that the numerical integration in Eq. (50) in combination with (39)(40) and (52) reproduces the earlier result for ψ(α, x, t) obtained in Minkowski space, Eq. (27). Note in particular that the vanishing of the LFWF at the endpoints α = ±1 is automatic. The branch cuts in Fig. 6 can also cross the real √ x axis, in which case it is no longer straightforward to compute the light-front distributions (51) numerically and one must deform the integration contour in in which case no contour deformation is necessary. If √ t moves to the imaginary axis but does not satisfy this constraint, the branch cuts become increasingly circular. Imaginary √ t with 0 < Im √ t < 1 corresponds to the physical situation 0 < M < 2m for bound states. In this case the only integration contour starting from the origin also coincides with the imaginary axis, which makes a numerical integration impossible. The strategy is then to compute the quantities of interest for complex values of √ t and take the limit Re √ t → 0 in the end.

IV. DYNAMICAL CALCULATION OF LIGHT-FRONT WAVE FUNCTIONS
A. Bethe-Salpeter equation In general, the BS amplitude Γ(q, P ) = Γ(x, ω, t, α) is not available in a closed form but emerges dynamically from the solution of the homogeneous BSE, which is illustrated in Fig. 7 and reads Γ(q, P ) = d 4 q (2π) 4 K(q, q , P ) G 0 (q , P ) Γ(q , P ) . (55) G 0 (q, P ) is the propagator product from Eq. (39) and K(q, q , P ) is the one-boson exchange kernel for a scalar particle with mass µ, Because g in the scalar theory is dimensionful, it is convenient to define a dimensionless coupling constant c and the mass ratio β: Details on the solution of the scalar BSE can be found in Ref. [66]; the practical complication in our present case is the appearance of the momentum partitioning α = 0.
To this end, we set q = k + αP/2 as before, and with q = k + αP/2 also for the loop momentum we have q − q = k − k . We work in the rest frame defined by Eq. (49) and Writing Γ, K and G 0 in terms of Lorentz invariants, the kernel becomes and the propagator product from Eqs. (39) and (52) which carries the α dependence is × 1 x Shifting the integration from d 4 q to d 4 k , the integration measure is and the BSE becomes One can see that the mass m drops out from the equation so that only c and β remain parameters. In particular, c only multiplies the right-hand side of the BSE and thus its eigenvalue spectrum, so that only t, α and β enter as external parameters.
The homogeneous BSE is an eigenvalue equation as it has the formal structure where the eigenvalues λ i depend on √ t = iM/(2m) and the mass ratio β. Because α is just a momentum partitioning parameter, the eigenvalues are independent of α which we also confirmed in our numerical calculations. On the other hand, the dependence of the eigenvalues

on
√ t determines the physical spectrum: the BSE has solutions for λ i (t i ) = 1, which determine the masses M i of the ground and excited states. Because the coupling strength c in the scalar model is just an overall parameter that is not constrained by anything, one can always tune it to produce physical solutions by setting c = 1/λ i (t i ). In the following it is therefore not relevant where a bound state appears; it is sufficient to know that for appropriate values of c they may appear for imaginary √ t in the interval 0 < Im √ t < 1, which corresponds to M < 2m below the threshold as illustrated in Fig. 8.

B. Singularities and contour deformations
The practical difficulties in solving the BSE, in particular in view of extracting the LFWF, arise from its singularity structure: The BSE must be solved for √ x inside the region R shown in Fig. 6, because only this region returns the correct result for the LFWF (50) when integrating over real ω ∈ (−∞, ∞).
The BSE is an integral equation, where the amplitude Γ(x, ω, t, α) is fed back during the iteration. This means it must be solved along a path √ x that coincides with the integration path in √ x . The natural path is the straight line between √ x = 0 and √ x → ∞, but because of the observations above this path must be deformed into the complex plane to lie entirely within R .
The BSE kernel (59) has a pole at After integrating over ω and y, the pole turns into a branch cut in the complex √ x plane. For ω, ω , y ∈ [−1, 1], also the variable Ω is in the interval Ω ∈ [−1, 1]. The task of picking the correct residues in analogy to Sec. II then translates into avoiding the branch cuts in the √ x integration. The resulting cuts for different values of β > 0 are shown in Fig. 9: For a given point √ x, they are confined to a region bounded by the circle with radius | √ x| and a line at √ x = Re √ x (see [66] for a detailed discussion). Because the same happens for every point along the integration path (blue dashed line in Fig. 9), avoiding all possible cuts for different β values implies that once the path has reached a particular value √ x 1 , it can only proceed if both the real part of √ x and its absolute value do not decrease -otherwise one would turn back into a region populated by branch cuts from the previous point √ x 1 . This limits the possible contours in √ x on which the BSE is solved: both Re √ x and | √ x| must never decrease along such a contour. The dashed curve in Fig. 9 is an example for a contour satisfying these constraints.
The singularities arising from the propagators in (60) produce the same cuts √ x ± as in Eq. (53), except that the integration variable ω does not span the full real axis but only the interval ω ∈ [−1, 1]. These are already taken care of by choosing a path inside the region R . In particular, Eqs. (53) and (64) become identical if one replaces Ω → ω, β → 1 and √ x → A ± . The resulting cuts have analogous shapes as in Fig. 9, where the two circles have radii |1±α| | √ t|. Therefore, a path that ensures safe passage is the one that connects the origin with the point (1 + |α|) √ t and then turns back to the real axis by increasing its absolute value as shown in Fig. 9.
If we are not interested in the LFWF but only the BSWF, then only the interval ω ∈ [−1, 1] is relevant for the integration. From the discussion in Appendix B, in that case Eq. (54) relaxes to in which case no contour deformations are necessary for any value of α. For α = 0, this reduces to Im √ t < 1 as discussed in [66]. Note that the parameter γ only applies to the monopole example, whereas in the BSE solution the BS amplitude is calculated dynamically. The BS amplitude Γ(x, ω, t, α) may dynamically generate singularities in the course of the iteration. The ω dependence does not cause any trouble because the interval ω ∈ [−1, 1] is free of singularities. However, in principle the equation can generate singularities in the complex √ x plane, like the monopole amplitude (40) does (the corresponding cut in Fig. 6 is the yellow curve for √ x 0 ).
As explained in Appendix B, a singularity q 2 /m 2 = −u 2 does not affect the contour deformation if arg( √ t) < arg(iu) .
As long as this condition is satisfied, a contour deformation connecting the origin √ x = 0 with the point √ x = (1 + |α|) √ t and returning to the real axis as described above is always possible. For real singularities (imaginary iu) the calculation is therefore feasible for any √ t ∈ C in the upper right quadrant, whereas for complex singularities the calculable domain in √ t is restricted to a segment bounded by the first singularity √ t = iu as shown in Fig. 10. In other words, one does not need to know the actual singularity locations as long as they are restricted to the region (66). The same statement also holds for propagators with complex singularities which we will analyze in Sec. V B.

C. Bethe-Salpeter amplitude
By implementing the contour deformations described above, the BSE (62-63) can be solved for any √ t ∈ C. Fig. 11 shows the (inverse of the) largest eigenvalue λ 0 corresponding to the ground state. As Re √ t becomes smaller, one can see the formation of the branch point at the threshold √ t = i. The ground state is determined by the conditions Re 1/λ 0 = c and Im 1/λ 0 = 0. Depending | ) x, ω, t, α Γ( | on the coupling parameter c and the mass ratio β, it is either a bound state (with β = 4 in Fig. 11 this happens for 6 c 11), a tachyon (larger c), or a virtual state on the second Riemann sheet (smaller c), as shown in [66] by solving the corresponding scattering equation.
In the following we compare our results with those obtained using a Nakanishi representation, which has been frequently used in the calculation of light-front quantities [26][27][28][29][30][31][32][33][34]36]. The central task in that case is the calculation of the Nakanishi weight function for a given interaction model, from where all further quantities (BSWF, LFWF, etc.) are obtained; see Appendix A for details. In Fig. 11 one can see that the eigenvalues obtained from both methods are in excellent agreement.
The typical shape of the solution for the BS amplitude Γ(x, ω, t, α) is shown in Fig. 12 for a fixed value of α and t. The amplitude falls off like a monopole ∝ 1/x in the x direction, whereas the ω dependence is very weak and difficult to see in the 3D plot. From Eqs. (59)(60)(61)(62) it follows that the amplitude is invariant under a combined flip α → −α and ω → −ω, because together with ω → −ω this operation leaves the BSE invariant. Also the α dependence is rather modest, as shown in Fig. 13 for exemplary values of x, ω and t. This suggests that the α dependence of the LFWF (50) is largely carried by the propagator product G 0 entering in the BSWF Ψ = G 0 Γ.
On the other hand, the amplitude can not be entirely independent of ω because then the LFWF would no longer vanish at the endpoints α = ±1. In that case, Eq. (50) reduces to with ω 0 = ∓(1+x+4t)/(4 √ xt). If Γ were independent of ω, this integral would not converge; moreover, by Cauchy integration it can only vanish if the amplitude has sin- Re Γ gularities for some ω ∈ C that lie on the same (upper or lower) half plane as the pole at ω = ω 0 . This requirement coincides with the conditions defining the region R in Fig. 6 for α = ±1, so the vanishing of the LFWF at the endpoints is automatic as long as √ x ∈ R .

D. Light-front wave function
We now proceed with the calculation of the LFWF according to Eq. (50). At this point, the solution for the BS amplitude Γ(x, ω, t, α) and thus the BSWF Ψ(x, ω, t, α) has been determined numerically for α ∈ [−1, 1], ω ∈ [−1, 1] and √ x along a contour inside the region R (cf. Fig. 6). The additional complication is that in the integration to obtain the LFWF one needs to know the dependence on ω over the whole real axis and not just inside the interval ω ∈ [−1, 1]. To analytically continue the ω dependence to the entire real axis, we employ the Schlessinger point method (SPM) [67] which has found widespread recent applications as a high-quality tool for analytic continuations into the complex plane, see e.g. [66,[68][69][70][71][72][73][74][75]. It amounts to a continued fraction    which is simple to implement by an iterative algorithm. Given n input points ω i with i = 1 . . . n and a function whose values f (ω i ) are known, one determines the n coefficients c i and thereby obtains an analytic continuation of the original function for arbitrary values ω ∈ C. The continued fraction can be recast into a standard Padé form in terms of a division of two polynomials.
In our case, f (ω) is the BS amplitude Γ(x, ω, t, α) for fixed values of x, t and α. The input points ω i lie inside the interval ω i ∈ [−1, 1], and the analytic continuation is performed for ω ∈ R. The LFWF (50) is finally obtained by integrating the resulting BSWF over ω ∈ R.
The resulting LFWF for an exemplary value of √ t is shown in Fig. 14. The falloff in x, the symmetry in α and the vanishing at the endpoints α = ±1 are clearly visible. We did not employ any polynomial expansion to facilitate the calculation; the result is the plain integral from Eq. (50). In Fig. 14 and also the subsequent plots we employed an additional SPM step analogous to Eq. (68) to transform the LFWF, which is obtained along a complex contour in x, to the real axis x ∈ R + . Figs. 15 and 16 give a detailed view of the α dependence of the LFWF (for selected values of x) and its x dependence (for selected values of α), respectively. The curves correspond to three values of √ t ∈ C, two of which lie below the threshold (Im √ t = 1) and one above. The results agree very well with the Nakanishi method which is applicable below the threshold. The region above the threshold is unphysical since there are no physical poles on the first sheet (the scalar model also does not produce resonances but instead virtual states on the second sheet [66]), but one can see in the plots that the LFWF is well-defined and calculable for any value of √ t.

13
) A necessary condition for our strategy to be applicable is that the BS amplitude does not have singularities for ω ∈ R. We checked that this is indeed not the case. Although the SPM cannot produce branch cuts but only poles, the resulting pole structure indicates the existence of cuts connecting the origin with the point ω = √ t for |ω| > 1. However, we note that even if the BS amplitude had singularities for real ω, this would not invalidate the approach because one would merely need to rotate the integration contour in Eq. (50) accordingly. In that case the region R in Fig. 6 would also change and one would need to solve the BSE along a modified path in √ x.
We close this section with some remarks on the numerical stability. The two relevant variables in the BSE (62) are the radial variable √ x, which takes values along the deformed contour shown in Fig. 9, and the angular variable ω ∈ [−1, 1]. For each of the three segments in √ x we employed a Gauss-Legendre quadrature with N x integration points in total. For the ω direction we employed a Chebyshev quadrature with N ω points. The necessary number of integration points depends on the external variables √ t and β, where β is the mass ratio from Eq. (57). If √ t comes close to the imaginary axis, one needs more integration points since the poles in the integrand and resulting cuts move closer to the imaginary axis and the integration path (left panel in Fig. 6).
Similarly, for small values of β the cuts from the kernel in Fig. 9 become increasingly circular and for β = 0 the points along the integration path itself become the branch points. For not too small values of Re √ t 0.1 and β 1, we find that N x = 96 and N ω = 64 is sufficient to reach sub-permille precision for the BSE eigenvalues in Fig. 11. For smaller √ t or β these numbers increase, and the numerical stability is typically more sensitive to N x than N ω . The same values are sufficient to achieve numerical convergence for the LFWF, which is also more sensitive to N x than N ω . Finally, for the analytic continuation in ω using the SPM we found an optimal value N SPM = 24 for the number of input points to achieve agreement with the Nakanishi method.

E. Light-front distributions
The remaining task is to compute the distribution amplitude φ(α) and distribution function u(α) as given in Eq. (51). This amounts to an integration of the LFWF over the transverse momentum variable x, which does not need to lie on the real axis since one can integrate along the same deformed contour on which the LFWF was obtained.
The resulting distribution amplitude φ(α) is shown in Fig. 17, again for three values of √ t ∈ C and compared to the results obtained with the Nakanishi method. It inherits the same properties as the LFWF; once again, the figure shows the plain numerical result without any expansion in moments. Also here the contour deformation method is in good agreement with the results from the Nakanishi method.
Although we did not employ them in our calculations, one usually defines Mellin moments ξ m for the reconstruction of hadronic distribution functions through the momentum fraction ξ, cf. Eq. (14): For the plots we used the normalization dα φ(α) = 2, which ensures ξ 0 = 1 for the zeroth moment. The first moment ξ = 1/2 is the mean value of the distribution which is centered around ξ = 1/2 or α = 0. All results presented so far have been obtained for √ t ∈ C in the upper right quadrant. As discussed above in connection with Figs. 6 and 8, it is numerically not possible to calculate the LFWF and PDA directly for imaginary √ t corresponding to real masses 0 < M < 2m, because in that case the only allowed integration path would coincide with the branch cuts.
To compute the PDA in the physical region, we employ the SPM from Eq. (68) to analytically continue the results for complex √ t to the imaginary axis in √ t. Here we employed a Chebyshev expansion where U n (α) are the Chebyshev polynomials of the second kind, and performed the analytic continuation for the Chebyshev moments φ n . We chose an input region Re √ t ≥ 0.1 for the SPM and changed the number of input points from N = 10 to N = 50 in steps of ∆N = 2. The resulting PDA at √ t = 0.5i is shown in Fig. 18, where the mean value is the average over the different input points and the error bands are the 1σ and 2σ deviations. Also here the results match with the Nakanishi method. The analogous result for the distribution function u(α) is shown in Fig. 19 and displays a somewhat larger error from the analytic continuation.
The SPM is reliable for imaginary values of √ t sufficiently below the threshold, whereas above the threshold it would attempt to continue to the second Riemann sheet and thus the results deteriorate as one moves closer to the threshold. We emphasize, however, that an analytic continuation in √ t is not necessary in principle since all quantities of interest can be calculated directly but at the expense of an increasing numerical cost for Re √ t → 0. In any case, for the scope of this exploratory study it is clear that the contour-deformation method is well suited to compute light-front distributions in a quantitatively reliable manner.

V. GENERALIZATIONS
In this section we explore two extensions of the contour-deformation method within the scalar model to bridge the gap towards possible future applications in QCD: One is the generalization to unequal masses in the BSE and the other is the implementation of complex conjugate propagator singularities.

A. Unequal masses
A straightforward generalization is the case of unequal masses m 1 = m 2 of the constituents in the two-body BSE. To this end we write such that In this way the mass parameter m drops out from the BSE like before, which in turn depends on ε.
Alternatively, one may start from Eq. (9) with the original momenta q and P (instead of k and P ) with an arbitrary momentum partitioning parameter ε ∈ [−1, 1]. For unequal masses, the choice ε = (m 1 − m 2 )/(m 1 + m 2 ) maximizes the domain in t where the BSE can be solved without contour deformations, namely Im √ t < 1. When introducing the momentum k through Eqs. (11)(12), the momentum partioning parameter ε drops out from all subsequent equations in favor of α which takes its place, and setting q + = (α − ε)P + /2 in the LFWF leads to the same result as before, Eq. (50). Thus, ε only appears in the propagator product G 0 (q, P ) = 1 q 2 1 + m 2 through the masses m 1 and m 2 and assumes the role of the physical mass difference as defined above.  Fig. 2(a) does not return the correct limit δ → 0 (left), in contrast to the one according to Fig. 2(b) (right).
Writing u ± = 1 ± ε, Eqs. (52)(53) generalize to with A ± = ∓(1 ± α) √ t as before. This does not require any modifications in the contour deformation: a path connecting the origin with the point (1 + |α|) √ t and turning back to the real axis by increasing its real part and absolute value is sufficient to avoid the cuts for any u ± ∈ R, so it is equally applicable in this case.
In the unequal-mass case, the BSE eigenvalues are still independent of α but they change with ε, which is a physical parameter in the system. The BS amplitude is still invariant under the combined flip α → −α, ω → −ω and ε → −ε. Fig. 20 shows the resulting LFWF for different values of ε and selected values of x and t. In the unequal-mass case, the LWFW and corresponding lightfront distributions are no longer centered at α = 0 but tilted towards the heavier particle.

B. Complex propagator singularities
Another generalization concerns complex singularities in the propagators, which is the typical situation for QCD propagators obtained from functional methods within truncations, see e.g. [76][77][78][79][80]. For a single complex pole pair the propagator becomes D(q 2 ) = 1 2 1 q 2 + m 2 (1 + iδ) which for δ = 0 reduces to a real mass pole. Writing u 2 = 1 + iλδ withλ = ±1, Eqs. (52)(53) generalize to again with A ± = ∓(1 ± α) √ t. This can be generalized to the case discussed in Appendix B, cf. Fig. 24: For any singularity at q 2 /m 2 = −u 2 ∈ C in the upper half plane, iu lies in the right upper quadrant. A contour connecting the origin with the point (1 + |α|) √ t and turning back to the real axis by increasing its real and absolute value is sufficient to avoid the cuts as long as arg( √ t) < arg(iu) like in Fig. 10. If this condition is not satisfied, the cuts twirl in the opposite direction and different cuts may overlap so that no contour can be found, i.e., the region R no longer connects the origin with infinity. Therefore, the appearance of complex singularities does not impede the contour deformation method but merely restricts the calculable domain in the complex √ t plane. Some care needs to be taken, however, when working with complex conjugate poles in Minkowski space using residue calculus like in Sec. II B (see also Ref. [81] for a discussion). In that case the literal i prescription, where the i factors appear in the denominators and one integrates over the real w axis, is no longer meaningful. The w ± poles in the complex w plane are now separated from the real axis by a finite distance proportional to δ, and the infinitesimal term does not change that. An integration over the real axis as shown in the left of Fig. 21 therefore does not give the correct limit for δ → 0. By contrast, the integration path in the right panel corresponding to the i prescription in Eq. (31) yields a proper analytic continuation and returns the correct limit δ → 0. This is, of course, equivalent to the Euclidean integration path as long as one stays in the region R .   The results are very similar, which confirms that the LFWF and the quantities derived from it are insensitive to whether the propagator poles are real or complex. In conclusion, one can tackle complex singularities just like real singularities using contour deformations.

VI. SUMMARY AND OUTLOOK
In this work we explored a new method to compute light-front quantities based on contour deformations and analytic continuations. We applied the method to calculate the light-front wave functions and distributions for a scalar model, whose Bethe-Salpeter equation we solved dynamically, and found excellent agreement with the well-established Nakanishi method.
The method is quite efficient and has several advantages, as it is not restricted to bound-state masses below the threshold and it can also handle complex singularities in the integrands. Although we exemplified the technique for a situation where the propagators are known explicitly, in principle one does not even need to know the singularity locations as long as they are restricted to a certain kinematical region. Since the method is independent of the type of correlation functions, in principle one can extend it to the calculation of parton distributions such as PDFs, TMDs and GPDs in the future.

Appendix B: Cuts
In this appendix we analyze the branch cuts in Eq. (53) in more detail. Their general form is where A = a + ib and C = c + id are complex numbers (a, b, c, d ∈ R) and ω ∈ R is varied over the real axis. This defines a curve in the complex plane whose shape depends on A and C.
Because √ x ± is the solution of the equation solving for ω and setting Im ω = 0 gives an alternative parametrization of the cuts which separate the two regions in √ x: Because only C 2 enters in the formulas, we restrict C to the upper right quadrant, i.e., c > 0 and d > 0. For C in the lower right quadrant the whole structure is mirrored along the line √ x = λA with λ ∈ R.
In the six panels of Fig. 23, |C| is fixed and the phase of C varies from 0 to π/2 (in the figure we set A = i). To facilitate the construction, the open points indicate the values √ x = ±AC, and we draw the circles with radius |AC| and lines passing through the points AC 2 (both in orange). The shape of √ x ± is then as follows: The line connecting the origin with √ x = A (in Near the origin, the two branches follow the line in the direction √ x = AC 2 . For C ∈ R + , this is the direction of A (imaginary axis in Fig. 23, top left panel).
When we rotate C into the complex plane with arg C > 0, the line also rotates and drags the curves with it.
As C is rotated further, then for Im (AC 2 ) < 0 the curve √ x ± eventually crosses the real axis (filled squares in Fig. 23). For A = i this means Re C 2 = c 2 − d 2 < 0, i.e., arg C > π/4. For general A the crossing happens at with r = a/b. In a situation where one needs to integrate √ x from zero to infinity, a contour deformation is thus necessary. The angular region between the lines λA and −λAC 2 (λ > 0) is guaranteed to be free of any cuts and thus an integration path leading away from the origin in that direction is safe.
Finally, if C becomes imaginary (C = id), the orange line has rotated by π and is again the direction of A (bottom right panel in Fig. 23). For |ω| < d, the curves √ x ± lie on the circle with radius |AC|; for |ω| > d they follow the line in the direction of A.
The cuts in Eq. (53) correspond to C 2 = −1 − u 2 /A 2 , where u 2 = 1 for the propagator cuts and u 2 = γ for the monopole cut, and A ∝ √ t up to real constants. Eq. (B1) then turns into and Eqs. (B2-B3) become For general values of u ∈ C, the limiting cases in Fig. 23 correspond to where λ ∈ R. The location of u then divides the complex plane into four quadrants delimited by the lines λiu and λu, where the cuts generated by the points A inside these quadrants take intermediate shapes like in Fig. 23. The situation for iu in the upper right quadrant is illustrated in Fig. 24: depending on whether arg (A) < arg (iu) or arg (A) > arg (iu), the cuts twirl in one or the other direction.
Writing iu = g + ih with g, h ∈ R, C 2 = −1 − u 2 /A 2 entails so that Eq. (B4) for the zero crossing results in As long as |A| < h 2 − g 2 + 2ghr, the cuts do not cross the real axis. Also relevant is the condition for the cuts to cross the real axis if ω is restricted to the interval ω ∈ [−1, 1]. In that case Eq. (B8) must be satisfied for |ω| < 1, which implies b 2 > h 2 ; or in other words, the cuts do not cross the real axis as long as |Im A| < |Im iu|.
Another question concerns the general singularity locations in Eq. (52), q 2 1,2 = −m 2 u 2 for the propagators and q 2 = −m 2 u 2 for the BS amplitude, with u ∈ C. This is relevant for the propagators entering in G 0 (q, P ), because in general they may not be known in the whole complex plane, and for the BS amplitude Γ(q, P ) which may dynamically generate singularities in the course of the BSE solution. Both cases translate to Eq. (B5) and Fig. 24 for iu k (k = 1, 2, 3, . . . ) in the upper right quadrant. Any such singularity generates a branch cut in the complex √ x plane, where the region R in Fig. 6 arises from the intersection of all cuts for a given value A ∝ √ t. For example, if the amplitude generates only real singularities (then the respective iu k are imaginary) but the propagators produce complex singularities (so that iu k ∈ C), then as long as arg(A) < min {arg(iu k )} (B9) a contour deformation is always possible because there is always a path in √ x ∈ R that connects the origin with infinity. Vice versa, for arg(A) > min {arg(iu k )} like in the right panel of Fig. 24, the different cuts may intersect such that no such path can be found. In the extreme case where some singularity iu moves to the positive real axis, the equations can only be solved for √ t ∈ R + , whereas in the opposite extreme case where all singularities iu k are confined to the imaginary axis like in Fig. 6, a contour deformation is always possible and the equations can be solved for any √ t ∈ C.