On coupled-channel dynamics in the presence of anomalous thresholds

We explore a general framework how to treat coupled-channel systems in the presence of overlapping left and right-hand cuts as well as anomalous thresholds. Such systems are studied in terms of a generalized potential, where we exploit the known analytic structure of t- and u-channel forces as the exchange masses get smaller approaching their physical values. Given an approximate generalized potential the coupled-channel reaction amplitudes are defined in terms of non-linear systems of integral equations. For large exchange masses, where there are no anomalous thresholds present, conventional $N/D$ methods are applicable to derive numerical solutions to the latter. At a formal level a generalization to the anomalous case is readily formulated by use of suitable contour integrations with amplitudes to be evaluated at complex energies. However, it is a considerable challenge to find numerical solutions to anomalous systems set up on a set of complex contours. By a suitable deformations of left-hand and right-hand cut lines we managed to establish a framework of linear integral equations defined for real energies. Explicit expressions are derived for the driving terms that hold for an arbitrary number of channels. Our approach is illustrated in terms of a schematic 3-channel systems. It is demonstrated that despite the presence of anomalous thresholds the scattering amplitude can be represented in terms of 3 phase shifts and 3 in-elasticity parameters, as one would expect from the coupled-channel unitarity condition.

In particular coupled-channel systems involving the nonet of vector mesons with J P = 1 − or the baryon decuplet states with J P = 3 2 + can only be studied with significant results after such a framework has been developed. The latter play a crucial role in the hadrogenesis conjecture that expects the low-lying resonance spectrum of QCD with up, down and strange quarks only, to be generated by final state interactions of the lowest SU(3) flavor multiplets with J P = 0 − , 1 − and J P = 1 2 + , 3 2 + [13,14,16,17,[29][30][31][32][33][34][35]. The coupled-channel interaction of such degrees of freedom leads to a plethora of subtle effects, like numerous anomalous thresholds [36][37][38].
The physical relevance of anomalous threshold effects has been discussed recently in [39][40][41][42][43][44]. To the best knowledge of the authors there is no established approach available that can treat such phenomena in coupled-channel systems reliably. In the previous works which attempted to deal with such systems the strategy was to perform an analytic continuation of an N/D ansatz for the reactions amplitudes in the external mass parameters as to smoothly connect a normal system to an anomalous system. This was studied for two channel systems only [37,38]. Even there the first study of Ball, Frazer and Nauenberg [37] was rejected by the later work of Greben and Kok as being incorrect [38]. So far we have not been able to track any numerical implementation of either of the two schemes [37,38]. Following this strategy an extension to a truly multi-channel systems appears prohibitively cumbersome.
A powerful framework to study coupled-channel systems is based on the concept of a generalized potential. A partial-wave scattering amplitude T ab (s) with a channel index a and b for the final and the initial state respectively is decomposed into contributions from leftand right-hand cuts where all left-hand cut contributions are collected into the generalized potential U ab (s). For a given generalized potential the right-hand cuts are implied by the non-linear integral equation whereρ cd (s) is a channel dependent phase-space function and all integrals are fors on the real axis. By construction any solution of (1) does satisfy the coupled-channel s-channel unitarity condition for normal systems. Typically, the matching scale µ 2 in (1) is to be chosen in a kinematical region where the reaction amplitude can be computed in perturbation theory [13,17,20]. For normal systems the non-linear and coupled set of equations (1) can be solved by an N/D ansatz [20,21,[45][46][47][48][49][50][51][52]. While the general framework is known from the 60's of last century only recently this framework has been successfully integrated into an effective field theory approach based on the chiral Lagrangian [20][21][22][23][24][25][26]. As it stands (1) this approach breaks down once a coupled-channel system involves unstable particles or anomalous thresholds arise. In this work we will construct a suitable adaptation that overcomes this gap.
Our formal developments will be illustrated by a schematic three-channel model where explicit numerical results for key quantities will be presented along the way. In Section II and III we discuss the analytic structure of the coupled-channel reaction amplitudes and propose an efficient deformation of the left-and right-hand cut lines such that anomalous systems can be dealt with. In addition our schematic model is specified. In Section IV an ansatz for the solution of the non-linear integral equations on the complex contour lines is derived. In the next step in Section V we consider the limit where the complex contour lines are deformed back towards the real energy line. The main formal result of our work is derived in Section VI, where a set of linear integral equations for anomalous systems that is suitable for a numerical implementations is established. We close with a summary and outlook in Section VII.

II. ANALYTIC STRUCTURE OF PARTIAL-WAVE SCATTERING AMPLITUDES
In a first step we will recall a spectral representation for a generic t-channel and u-channel term as shown in Fig. 1. In our previous work [28] we established the following general form with properly constructed spectral weights (t) ±,ab (m 2 , m 2 t ) and (u) ±,ab (m 2 , m 2 u ). The contour functions c (t) ±,ab (m 2 ) and c (u) ±,ab (m 2 ) depend on the masses of initial and final particles of the given reaction ab for which we use the convenient notation While the derivation of the spectral weights (t) ±,ab (m 2 , t) and (u) ±,ab (m 2 , u) is quite cumbersome the identification of the contour functions c (t) ± (m 2 ) and c (u) ± (m 2 ) is straight forward. Owing to the Landau equations any possible branch point of a partial-wave amplitude must be associated with an endpoint singularity of the partial-wave projection integral that involves some Legendre polynomials in cos θ of the scattering angle θ. This leads to the well known result An anomalous system arises if either of the spectral weights  i (m 2 ) approaches any of the thresholds or pseudo thresholds of the given reaction ab. In this case the decomposition (1) breaks down: an anomalous threshold behavior is encountered. The latter is characterized by a particular branch point µ A ab of the partial-wave amplitude that is associated with the given amplitude ab. For simplicity of the presentation we consider in (5) an anomalous threshold behavior at a normal threshold point only. A similar phenomenon may occur at a pseudo threshold. Also the case µ A ab > Max{m a + M a , m b + M b } is not excluded in general. For both cases, our approach can be adapted in a straight forward manner. What is excluded, however, from general considerations, that an anomalous threshold arises for a diagonal reaction with a = b. In the previous works an analytic continuation in the external mass parameters as to smoothly connect a normal two-channel system to an anomalous two-channel system was attempted [37,38]. After all using this method Mandelstam gave a transparent presentation of the anomalous threshold phenomenon in a typical one-loop diagram [36]. However, we feel that an application of this method directly to multi-channel reaction amplitudes, as attempted in [37,38], appears futile due to the proliferation of branch points and cut structures that need to be properly deformed and followed up in various limits.
We argue that there is a significantly more transparent path to arrive at a framework that is capable to treat coupled-channel systems in the presence of anomalous thresholds. Given a multi-channel system various anomalous thresholds may appear in different reactions. Our starting point, is based directly on the general representation (2), which already established a suitable analytic continuation of the tree-level potential terms valid for arbitrary exchange masses. As was emphasized in our previous work [28] the reason why a one-loop diagram develops an anomalous threshold can be read off easily from its driving tree-level potential terms.
We need to be somewhat more specific on the generic form of the generalized potential.
The tree-level potential takes the generic form where the left-hand contour, L ab , is a union of contours required for the spectral representation of a general t-or u-channel exchange process as recalled with (2). In order to achieve our goal it is instrumental to separate that left-hand contour into a normal and an anomalous part where the part A ab starts at the anomalous threshold (µ A ab ) 2 slightly above the real axis, extends to the normal threshold point Min{(m a + M a ) 2 , (m b + M b ) 2 } and returns at (µ R ab ) 2 to the anomalous one slightly below the real axis. This is illustrated in Fig. 2. Note that according to the general representation (2) the anomalous contour line would be on real axis only. With Fig. 2 such a contour is deformed into the complex plane in a manner that leaves its contribution to the generalized potential unchanged for s outside the area encircled by the dashed line. While the values for µ A ab are determined by the specifics of the coupeld-channels interactions, there is some freedom how to choose the location of the return points µ R ab . Given the analytic structure of the spectral weight ρ U ab (s) the generalized potential U tree−level ab (s) does not depend on the choice of the return point µ R ab at energies s outside the dashed line of Fig. 2.
We point at a subtle issue such a contour separation is based on. Here we tacitly assumed that the anomalous threshold µ A ab is real. However, from the general representation (2) there is the possibility of a pair of complex conjugate anomalous points µ A± ab . In this case the associated contour lines have to be followed till both reach the common pseudo-threshold The latter we take as the anomalous threshold value in our work. This can be justified since in this case the spectral weight ρ U ab (s) can be shown to be analytic an left to that pseudo threshold point ats = (µ A ab ) 2 − . It is useful to identify an anomalous threshold value µ a associated with a given channel a. For the clarity of the development we assume in the following a strict channel ordering according to the nominal threshold value, i.e. we insist on where for later convenience we permit the equal sign for channels with the same two-particle states in different spin configurations only. We can now identify the desired anomalous threshold value. It is the minimum of all accessible anomalous branch points where the value µ A ab gives the anomalous threshold value of the partial-wave amplitude ab. Note that it does not necessarily follow that the channel ordering (8) implies µ a+1 ≥ µ a . If such a channel crossing with µ a > µ b for any pair ab with a < b occurs, we will further move the point µ a with µ a → µ b − such that we ultimately arrive at a strict channel ordering where again the equal sign is permitted for channels with the same two-particle states in different spin configurations only. It is useful to introduce a similar streamline of the plethora of return points. We choose universal return pointsμ a of the particular form where we consider a and b with µ A ab < (m a + M a ) 2 only. We can now introduce our anomalous contour A a that starts at µ 2 a slightly above the real axis, passes (m a + M a ) 2 and returns atμ 2 a to µ 2 a below the real axis. In turn we can now write where we will exploit that the first term in (12) is uncritical to the extent that it is analytic in the vicinity of any normal or anomalous threshold point. The challenging term is the second one, which is associated with a contour that entangles the normal threshold point The dots in (12) remind us that we did not yet discuss the generic form of contributions to the generalized potential from loop effects.
The key starting point is an adaptation of the non-linear integral equation (1). The integral in (1) starts at the point where the phase-space matrix ρ cd (s) vanishes ats = 2 and remains on the real axis going to infinity. In the anomalous case the integral must also start ats = (m c + M c ) 2 = (m d + M d ) 2 but will leave the real axis following suitable paths on higher Riemann sheets. The integral on the real axis has to be replaced by a contour integral. Eventually, somewhere on its way the contour path will touch the anomalous threshold ats = µ 2 c . Note that the case c = d appears only for systems with non-vanishing spin.
The reason for this complication is a consequence of the anomalous threshold behavior of the generalized potential. If we wish to separate the left-from the right-hand cuts in the reaction amplitudes, it is necessary to apply suitable deformations of the cut-lines. In the normal case the separation of left-from right-hand cuts is trivially implied by (1). The second term on the right-hand-side has a cut on the real axis extending froms = (m 1 + M 1 ) 2 to infinity. This is the right-hand cut. All left-hand cuts sit in U (s). As long as the cut lines of the generalized potential do not cross any of the right-hand cut lines on the real axis the non-linear integral equation (1) is well defined and can be solved numerically in application of conventional methods. This is not the case if the generalized potential develops an anomalous threshold. Without an appropriate adaptation the associated left-hand cut lines would cross the right-hand cut-lines. An avoidance is possible only if the integral on the real lines in (1) is replaced by contour integrals in the complex plane.
In Fig. 3 we illustrate the analytic structure of the reaction amplitudes defined with respect to suitably deformed left-and right-hand cut lines. The dashed red lines show the location of the anomalous left-hand cuts. We suggest a deformation of the right-hand cut lines as illustrated by the solid lines in the figure. Now, none of the anomalous left-hand cut lines A a crosses any of the right-hand contour lines C a .
Before writing down the generalization of (1) we need some more notation. The reaction amplitudes T ab (s) need to be evaluated slightly below and above any of the many different contour lines. With where the points s + and s − are an distance above and below the contour C c ats respectively. The value of is chosen smaller than the minimal distance of any of the horizontal lines in Fig. 3.
With this notation we can write down the desired generalization where C c = C d holds by construction. We consider the phase-space matrix ρ cd (s) to be analytic on the real axis at s > (m c + M c ) 2 . This implies that at s < (m c + M c ) 2 the function ρ cd (s) has a cut line. Note that none of the contour paths C c cross any of the cut-lines of the phase-space functions. It is emphasized that in the absence of an anomalous threshold in the generalized potential our generalization (14) reproduces the conventional expression (1).
The representation (14) appears deceivingly simple and ad-hoc, however, we derived it by an analytic continuations in the t-and u-channel exchange mass parameters. The starting point for this is provided by the general representation (2). First we assume all exchange particles to have a very large mass, so that all left-hand branch cuts are well separated from the right-hand branch cuts. Already at this stage we can deform the right-hand cuts from their conventional choice to take the form as indicated in Fig. 3 by the solid lines.
This contour deformation cannot change the value of any of the coupled-channel reaction amplitudes T ab (s) at Im s > 0, simply because along the cut deformations the reaction amplitudes are analytic by assumption. Here we exploit the fact that in the limit of very large exchange masses all left-hand branch cuts are moved outside the figure. In a second step we follow the left-hand cuts as they move right in Fig. 3 with the exchange masses getting smaller approaching their physical values [36-38, 48, 52]. At the end we claim that they can be presented by the dashed lines in the figure. In turn we take (14) as a transparent starting point of our presentation.
Our numerical examples will rest on a minimal model as it is implied by schematic treelevel interactions. Along the formal derivations in the next few sections we will illustrate the important auxiliary quantities in terms of a schematic three-channel model specified with In Fig. 4  In Fig. 5 we anticipate the usefulness of our formal developments and present the solution to the non-linear integral equation (14) as implied by (15). It is noted that our results do not depend on the particular choices for the return pointsμ a in (15). The reaction amplitudes show significant and non-trivial structures that are implied by the presence of the anomalous threshold effects in the generalized potential. While the amplitudes at subthreshold energies show various singular structures, they are smooth and well behaved in the physical region.
We emphasize that the subthreshold structures are not driven exclusively by the contribution from the generalized potential as shown in Fig. 4. There are in addition significant structures generated by the right-hand cut contributions. As to the best knowledge of the authors with

III. ANOMALOUS THRESHOLDS AND COUPLED-CHANNEL UNITARITY
Before we explain how to find numerical solutions to the non-linear integral equations (14) in the presence of anomalous threshold effects it is useful to pause and discuss a critical issue that we will be confronted with. Given a specific approximation for the generalized potential, a solution of (14) does not necessarily imply that the coupled-channel unitarity condition is fulfilled. From the latter we expect a representation of the reaction amplitudes T ab (s) in terms of a set of real quantities, channel dependent phase shift and inelasticity parameters φ a (s) and η a (s). It should hold where on general grounds one expects 0 ≤ η a ≤ 1. While any solution to (14) satisfies the time-reversal invariance condition T ab (s) = T ba (s), the Schwarz reflection principle, cannot be derived in general for right-hand cut lines off the real axis. It is not surprising that then the coupled-channel unitarity condition is not necessarily obtained.
Let us be specific and identify the generalized potential by tree-level t-and u-channel exchange processes, a typical strategy in hadron physics. Though we may solve the nonlinear system (14) in this case, the unitarity condition (14) will not be fulfilled once an anomalous threshold effect is encountered. A supposedly related problem was noted in the previous studies [37,38]. It was argued that the problem is caused by the neglect of second order contributions to the generalized potential. A minimal ansatz for a physical approximation to the generalized potential requires some additional terms to be constructed properly. However, no conclusive form of the latter was presented and illustrated in the literature so far. Conflicting suggestions were put forward in [37,38] for two-channel systems. Whether any of such forms lead to physical results remains an open challenge. So far there is no numerical implementation for a specific example worked out, at hand of which one may judge the significance of any of the two approaches. Unfortunately, for us both works [37,38] are rather difficult to follow. Nevertheless, we tend to agree with the conclusions of [38] that the approach advocated in [37] is incorrect.
An important achievement of the current work is the construction of a minimal second order term applicable for systems of arbitrary high dimensions. Our approach is based on the request that we arrive at the coupled-channel unitarity condition (16) and recover the Schwarz reflection principle (17). The key observation is that in the absence of the anomalous box term U box ab (s) the second order reaction amplitude T ab (s), as it would follow from (14), is at odds with (17). The extra term is unambiguously determined by the condition that its inclusion restores the Schwarz reflection principle (17). This leads to the following form where the spectral weight ρ A ab (s) characterizes the anomalous behavior of the tree-level potential. It is identified in analogy to the general representation (12). Unfortunately, a direct comparison of our result with the ansatz in [37,38] is not so easily possible. Nevertheless, we note that for the two-channel case our result (19) appears quite compatible with the ansatz discussed in [38].
We anticipate the phase shifts and in-elasticities as implied by our schematic model (15) as properly supplemented by the anomalous box term (18,19). We affirm that the reaction amplitudes as already shown in  It is instrumental to realize the quite different nature of the two contributions in (18).
Consider first the tree-level term in (18,19). The spectral weight ρ A ab (s) is real as s approaches the real line below s < Min{(m a + M a ) 2 , (m b + M b ) 2 }. In contrast it is purely imaginary for This follows from the general results presented in [28].
The corresponding generalized potential satisfies the Schwarz reflection principle with We turn to the second order term in (18,19). While for Re the two factors, ρ A ac (s) and ρ A cb (s), are real quantities as s approaches the real line, the phase-space factor ρ cd (s) turns purely imaginary in this case. This leads to the property and illustrates the particular feature of the anomalous box term our construction is based on. Given our schematic model (15) the anomalous box term U box ab (s) is illustrated with Fig.7. Like in the previous Fig. 6 the left-hand panel shows the non-vanishing potentials on the real axis, the right-hand panel the corresponding potentials slight below the real axis at It should be emphasized that there are further second order contributions to the generalized potential, however, they add to the 'normal' spectral weight ρ normal ab (s) in (12) only.
Since the latter terms do not jeopardize the coupled-channel unitarity condition, there is no stringent reason to consider such effects in an initial computation. Typically one may hope that the effect of the latter is suppressed in some suitable power-counting scheme. This should be so since higher loop effects are characterized by left-hand branch cuts that are further separated from the right-hand cuts. In turn such contributions to the generalized potential cannot show any significant variations at energies where the generalized potential is needed in (14).
We would speculate, that the ansatz (18) is quite generic, i.e. it should hold also for contributions including higher loop effects where we expect ρ anomalous ab,c (s) to be determined by ρ anomalous ac (s) and ρ anomalous cb (s) in analogy to (19).

IV. NON-LINEAR INTEGRAL EQUATION ON COMPLEX CONTOURS
The key issue is how to numerically solve that non-linear set of equation (14) and cross check its physical correctness. After all one may consider it merely as a definition of the generalized potential in the presence of anomalous thresholds. For a given approximated generalized potential U ab (s) we will device an appropriate N/D like ansatz that will eventually lead to a framework which is amenable to numerical simulations of (14).
We introduce a set of contour functions ς ab (s) defined initially on distinct contours C b , the choice of which depends on the channel index b where we apply the convenient ± notation introduced already in (13). Assuming the existence of such a set of functions ς ab (s) we seek to express the reaction amplitude T ab (s) in terms of them. This requires a few steps. Like in the previous sections we anticipate with It is emphasized that none of the function D ab (s + i ) depend on the particular choice of the return pointsμ ab in (15). The functions have a significant imaginary part starting at the anomalous threshold s ≥ µ 2 b . In the first step we study the analytic properties of the function where we applied the master equation (14). In Fig. 9 we provide such functions B ab (s + i ) as implied by our model (15) for s > (m b + M b ) 2 , a region which encompasses the physical domain probed by the phase shifts and in-elasticity parameters in (16). Here we again do not encounter any dependence on the return pointsμ a . However, the functions B a1 (s) show a strong variation close to the anomalous threshold point at s = µ 2 2 = 29 m 2 π . It is important to note that this is not propagated into the amplitude T 11 (s + i ) as is evident with Fig. 5 and Fig. 6. Since the contours C c partially overlap it is useful to decompose the contours with where the lower C − c contours are all on the straight line crossing the value s = −i . By construction the upper contours C + c do not overlap. Real and imaginary parts are shown with solid blue and dashed red lines respectively.
From this we conclude that a dispersion integral representation of the form can be assumed. In the next step we compute the discontinuity along the contours C + c with where we observe the cancellation of most terms in (27). This follows from the defining equations for ς ab (s) in (23) together with symmetry of the reaction amplitude It is left to compute the discontinuity of the B functions along the straight contour C A , i.e.
the second term in (26) is considered. For anys ∈ C A we derive where the cancellations in (29) follow from (23,28). We observe that owing to the identities (27,29) the two terms in (26) can be combined into a dispersion integral written in terms of the partially overlapping contours C c . It holds where the crucial identity in the last line of (30) is a consequence of the properly deformed contour lines as illustrated in Fig. 3.
With this we arrive at the anticipated representation of the scattering amplitude in terms of the spectral density ς ab (s). It holds The functions D ab (s) were already expressed in terms of ς ab (s) in (23).
We are one step before a more practical rewrite of the defining request (23) and in particular checking the consistency of the construction. After all we had to use both equations in the first line of (23) as to arrive at (31). Inserting our result (31) into (23) we will obtain two distinct equations, for which we have to show their equivalence. We derive where s is on the contour C b strictly. With s ± we introduce values of s slightly above and below the contour C b , i.e. it holds |s − s ± | < , where is chosen smaller than the minimal distance of any of the horizontal lines in Fig. 3. The two choices correspond to the two identities in the first line of (23). Since the numerator in (32) strictly vanishes ats = s both choices lead to identical results.
While with (32) we arrive at a mathematically well defined linear integral equation, it remains to construct a numerical solution to it. This is not quite straight forward and will require further developments. In the following we will analyze the linear system (32) in more detail and eventually establish a framework that can be used to numerically solve it on a computer. The key issue is to systematically perform the limit → 0 in the system of complex contours.

V. FROM COMPLEX CONTOURS TO REAL CONTOURS
We consider first the D function, for which its definition is recalled with in terms of the spectral weight ς ab (s). We consider the limit → 0, in which the complex contours C c all approach the real axis. If we are interested in values of s only that are below or above all right-hand cut lines, we may simplify the integral into Riemann sums on the real axis. This is achieved as follows.
We first note that we may consider ς ab (s) to be an analytic function in s, with various branch cuts. This follows from the integral representation (32). More precisely, if the linear system has a solution ς ab (s) with s ∈ C b then the equation (32) can be used to analytically continue ς ab (s) away from the contour line C b . The branch cuts are readily identified. First it carries the branch cuts of the phase-space function ρ cc (s) that is on the real axis strictly in our convention. Second, the cut-lines of the generalized potential U cb (s) for any c are inherited. The important observation is the absence of any right-hand cut-lines.
According to the cut lines summarized in Fig. 3 Fig. 3 and returns back. With this in mind we introduce where the integrals overs in (34) are on the real axis strictly. In (34) we apply the useful notation ς ± ab (s) = ς ab (s ± ) with Ims = 0 and s ± ∈ C b and Im (s + − s − ) > 0 and Re s ± =s .
which defines ς + ab (s) fors < (m b + M b ) 2 initially only, but is naturally extended up to the point ats =μ 2 b . Here we assume the availability of the analytic continuation of the function ς ab (s) from the nominal threshold value at s = (m b + M b ) 2 up to the return point of the left-hand cut at s =μ 2 b in Fig. 3. The particular location ofμ b > m b + M b is irrelevant. While the spectral weightsς ab (s) and ∆ς ab (s) depend on it, by construction the D function does not. Given (35) the limit → 0 in the right-and left-hand contours of Fig. 3 can be applied without changing the form of (34). In this limit the contributions of any vertical parts of the contour lines vanish. Such terms are already omitted altogether in (34). Note that modulo those vertical lines (34) is nothing but a regrouping of the various contour contributions in (33).
We illustrate the generic form of the spectral weightς ab (s) in our model (15). As shown in Fig. 10  There is a subtle point as to where to evaluate the D ac (s − − i ) function in the first line of (36). Since D ac (s) has a branch cut along C c ands − ∈ C c fors > µ 2 c it is necessary to specify whether we should evaluate the function below or above the cut. Our prescription follows unambiguously if we slightly deform the contours C c in (32). In Fig. 3 the solid line passing throughs − is deformed a bit towards, but still avoiding, the dashed lines above.
With this it is manifest from (32) that ς ab (s) is analytic along the horizontal line through s − . This should be so since ς ab (s) is analytic along all contours C c by construction. In turn our prescription D ac (s − − i ) is justified. Note that for the term D ac (s + ) no further specification is needed simply becauses + / ∈ C c . Here it always holds D ac (s + ) = D ac (s + i ).
We ask the reader to carefully discriminate the objects U c± ab (s) with s ∈ C c as introduced in (13) from the newly introduced object U ±c ab (s) withs defined on the real axis only in (36). In the following we will show that given the functionς ab (s) only the D ab (s) function can be computed unambiguously. Note that this requires the solution of a linear integral equation since ∆ς ab (s) requires the knowledge of the D ab (s) function. In order to solve this system it is useful to introduce some notation wheres is on the real axis. Note that due to the particular form of the anomalous box term U box ab (s) in (19) there is no contribution from the latter to ρ L ab (s). Our result (34, 36) is expressed as follows where all integrals overs are on the real axis. The bounds of the integrals are provided by the properties of ρ L cb (s) as summarized in (54). It remains to express D ab (s) in terms ofς ab (s) and ρ L ab (s). In the non-overlapping case this is readily achieved by an iteration in the index b. At b = max the second term in (38) does not contribute as a consequence of the second condition in (54). In turn we can compute D ab at b = max for all a. In the next step we study D ab (s) at b = max − 1, where now the second term in (38) turns relevant. However, here only the previously computed D ab (s) at b = max are needed. This process can be iterated down to the computation of D a1 (s). While this strategy is always leading to the correct result it is not very efficient for the following developments.
A more powerful framework can be readily established as follows. We introduce a Green's function L(x, y) via the condition where we suppress the coupled-channel matrix structure for notational clarity. All objects will be written in the correct order, so that the matrix structure can be reconstructed unambiguously for any identity presented below. The i prescription in the definition of (39) is inherited from the i prescription in (38). Given the Green's function D ab (s) can be expressed in terms ofς ab (s) with The representation (40) does not look very promising for numerical simulations since the Green's function is a highly singular object. However, a closed form can be derived in terms of six analytic matrix functions u L n (s) and U L n (s) with n = 1, 2, 3. The latter are determined by appropriate integrals involving the anomalous spectral weight ρ L (s) only. After some algebra we established the following form for the Green's function FIG. 11. The non-vanishing functions u L n,ab (s) of (43) for n = 1, 2 in our schematic model (15).
Real and imaginary parts are shown with solid blue and dashed red lines respectively.
We illustrate the general form of the complex objects u L n (s) and U L n (s) at hand of our schematic model (15). In Fig. 11 and Fig. 12 all elements are shown.
An important property of our result (40) is that the imaginary part of the spectral weight ς D ab (s) does not vanish in the presence of anomalous threshold effects. This follows from the results (40)(41)(42)(43). In turn the functions D ab (s) do not satisfy the Schwarz reflection principle even in the limit with → 0. We emphasize that this is unavoidable in our formulation.
Nevertheless we expect that our final reaction amplitudes T * ab (s) = T ab (s * ) will satisfy the Schwarz reflection principle, which plays a crucial role in the derivation of the coupledchannel unitarity condition. Note that (44) should not be taken as a surprise since we have already discussed that the anomalous box term U box ab (s) is at odds with (17). The relation (44) can be confirmed explicitly upon an expansion of the D-function in powers of the generalized potential. At second order there is a term that confirms (44).
We continue with the B function, for which its representation in terms of the spectral weight ς ac (s) and the generalized potential U ac (s) is recalled. Again we perform the limit → 0, in which the complex contours C c all approach the real axis.
Following the decomposition of the contour lines C c as introduced in our study of the D function in (34) we readily derive the corresponding rewrite where the integrals overs in (46) are on the real axis strictly. Again modulo contributions in (45) from vertical parts the contours C c both representations (45) and (46) are identical at any finite . In (46) we apply the convenient notation ς ± ac (s) and U ±c cb (s) introduced already in (35,36).
In the following we will express the spectral weightsβ c ab (s) and ∆β c ab (s) in terms ofς ab (s) and the ∆ς ab (s). For this we will have to consider the limit → 0 again, in which we find where the generalized potential U cb (s) is evaluated on the real axis strictly. It is important to realize that U cb (s) is needed in between the upper and lower anomalous cut lines (dashed lines in Fig. 2) even after the limit → 0 has been performed. For the anomalous spectral weight ∆β c ab (s) we derive three distinct contributions ∆β c ab (s) = ∆ς ac (s) U cb (s) which contribute depending on the various cases b > c, b < c or b = c. It is important to note that in (48) the contribution of the box term (19) requires particular attention. Its contribution to U cd (s) is to be taken slightly below the double cut structure of (19). In the following we will promote this to our convention. Whenever we use U cd (s) with s on the real axis the above prescripton is implied. We express our results (46,47,48) in the notation (54) as follows where all integrals overs are on the real axis. The bounds of the integrals are provided directly or by the properties of ρ L cb (s) as summarized in (54). We point at an important subtlety. While in the previous section we managed to express D ab (s) in terms ofς ab (s) at s > (m b + M b ) 2 , this is not possible for B ab (s). The term in (53) involving ς + ac (s) is required at µ 2 c < s <μ 2 c partially outside the domain whereς ac (s) was introduced in (34). It is instrumental to introduce a function ς • ac (s) defined on that interval, that coincides withς ac (s) in the region where their defining domains overlap. We do so as with µ 2 c < s <μ 2 c on the real axis strictly and Note that while the integral overs in (50) is on the complex contour C d it is on the real axis in (51) strictly.
It may be useful to point out that for a non-overlapping case systemμ n < µ n+1 for all n, it follows ς • ac (s) = ς + ac (s). We recall thatς ac (s) was introduced in (m c + M c ) 2 < s <μ 2 c to be the analytic continuation of ς + ac (s). In turn it holds: With this we arrive at the central result of this section in terms of ς • ac (s) and D c ad (s) as introduced in (50, 51) and It is useful to provide yet a further rewrite of the B-function that further streamlines the various prescriptions. We find with where we apply the particular notation U ±c ab (s) as introduced in (36). It should not come as a surprise that like the D-also the B-function does not satisfies the Schwarz reflection principle with even in the limit with → 0. This is readily verified. If expanded to second order in powers of the generalized potential there must be an anomalous contribution with (57) that cancels the effect of the anomalous box contribution (21). This is so since by construction the full reaction amplitude was constructed to satisfy (17) at least to second order in a perturbative expansion.

VI. LINEAR INTEGRAL EQUATION ON REAL CONTOURS
In the previous two sections we have expressed the D ab (s) and B ab (s) functions in terms of the spectral weightς ab (s). In this section we wish to establish a set of linear integral equations forς ab (s) given a generalized potential U ab (s). These equations will serve as an alternative formulation of (32), which is suitable for numerical simulations. After some necessary steps detailed in this section we will arrive at an integral equation of the form in terms of a set of analytic matrix functionsÛ (s),Û mn (s) and u L m (s), u R m (s). The latter will be expressed in terms of the generalized potential U (s). Note that in (58) we suppressed the coupled-channel indices. The terms are ordered properly so that the coupled-channel structure is correctly implied by standard matrix multiplication rules.
How to cast the contour integral equation (32) into an integral equation (58) where all integrals are on the real axis strictly withς ac (s) = 0 for s > (m c + M c ) 2 ? Several steps are required. The first task is to express the functions ς ± ab (s) in terms ofς ab (s) only. We begin with the consideration of ς − ab (s), which we evaluate according to the second identity in (23) with for which we consider the contour limit → 0 in Fig. 3. The reaction amplitude T b− cd (s) in (23) is expressed in terms of the B b− ab (s) function as evaluated in (53). Similarly for the required D b− ac (s) function we use (38). Then for → 0 and s >μ 2 b we find where we again use the particular notations ∆U ab (s) and U ±c ab (s) as introduced in (56, 36). We emphasize that s ands in (60) are strictly on the real axis as is implied by the limit → 0. It is important to realize that U cb (s) is evaluated in between the upper and lower anomalous cut lines, the dashed lines in Fig. 2. The double cut line of the box contribution is not probed here.
We continue with the more complicated case ς + ab (s) at s <μ 2 b . This time we start with the first identity in (23) with and again consider the contour limit → 0. For an evaluation of ς + ab (s) we will need the D ac (s) and B ab (s) functions evaluated slightly above the contour C b with s ∈ C b . The latter can not be deduced directly from the results of the previous two sections. This is so because sometimes the functions are required in between two right-hand cut lines, for which the results (38, 53) can not be applied. Some intermediate steps are required. Consider first the B b+ ab (s) term, the first contribution in (61). In the limit → 0 we can derive: where we consider s in (62) to be strictly real again. For the required kinematics with s <μ 2 b and → 0 it follows where we applied (62) in combination with (53). We point the reader to the different prescriptions s ± i in the various terms in (63) with The change in prescription is caused by the terms in the right-hand side of (62).
We proceed with the second term in (61). Here we need to evaluate D b+ ac (s) U b+ cb (s) in the limit → 0. Progress is based on the identity D b+ ac (s) = D ac (s + i ) which again follows in the limit → 0. From (65) it now follows where we assumed (m b + M b ) 2 < s <μ 2 b and → 0. Combining our results (63, 66) we arrive at the desired expression where we exploited the corresponding prescription changes in (63)  We turn to ς • ab (s) = ς + ab (s) as was introduced in (50) and required in the evaluation of ς + ab (s). In this case we consider the intermediate result where we now assumed µ 2 b < s < (m b + M b ) 2 and → 0. After some algebra we find The reader is pointed at the formal similarity of the two expressions (67, 69). This is in line with our construction with ς • ab (s) =ς ab (s) on (m b +M b ) 2 < s <μ 2 b . It remains to rewrite our result (67, 69) into a more practical form. This requires three steps. First we multiply (69) by the pseudo inverse of the phase-space matrix where the result is expressed in terms of the more convenient building blocks The integrals overs in (70) With (75) we derived a convenient basis for the derivation of a suitable integral equation to numerically solve forς ab (s). After a few more steps we will find with a well behaved integral kernelK db (s, s) and potential termÛ ab (s). It is emphasized that (75) is valid for s on the real axis strictly, in particular also for s < The auxiliary functionsÛ ab (s) are presented in Fig. 14 as derived from our model (15).
where we emphasize that each of the objects u R n (x) and U R n (x) has a coupled-channel matrix structure at any n = 1, 2, 3. The formal expressions for u R n (x) and U R n (x) can be extracted from our previous expression for u L n (x) and U L n (x) in (43). The close relations amongst the two, left and right, Green's function is a consequence of the identity which is reflected in (54, 71). This implies where we use U R n<4 (x) and u R n<4 (x) as previously introduced in (77). An analogous result holds for the action of the left Green's function L(x, y) on such a structure. This can then be applied to derive the integral kernel K(s, s) = U 4n (x) = dz π ρ L (z) ∆Ū eff n (z, x) z − x ,Ū 5n (x) =Ū eff n (x) , in terms of a set of analytic functionsŪ mn (s) and u L m (s), u R n (s). We note that U L m (x),Ū eff n (x) were already introduced in (43) where the object U eff + (x) was defined already in (74). Finally, the potential term takes the formÛ where with (71) we finally arrived at the anticipated result (58).
We illustrate the role of the auxiliary matricesÛ mn (s) as derived from our model (15).
Since there are quite a few functions, we focus on the particular combinations H L n (s) =  Fig. 13. Note that H R n,ab (s) are discontinuous only at the return points s =μ 2 b . The analogous property holds for the functions H L n,ab (s) at s > (m a + M a ) 2 . We should briefly summarize the general work flow how to derive the phase-shifts and in-elasticity parameters for a given model interaction. Given the driving termsÛ (s) and U nm (s) together with u L n (s) and u R n (s) as specified in (43,79, 80-84) we use the linear set of equations (58) to determine the functionς ab (s) for s > (m b + M b ) 2 . This is a numerical stable task since all driving terms in (58) are sufficiently regular in the needed domain.
In the next step we can compute the functions D ab (s + i ) in application of (40)(41)(42)(43). The functions B ab (s + i ) are evaluated from (53, 55), which, however, requires the knowledge of the functions ς + ab (s). This goes in two steps. First, given the functionsς ab (s) we can computê N ab (s) from (85) at subthreshold energies the desired object ς • ab (s) is available and we can finalize the derivation of the functions

VII. SUMMARY AND OUTLOOK
In this work we presented a novel framework how to deal with coupled-channel systems in the presence of anomalous threshold effects. The framework is formulated for an arbitrary number of channels and is suitable for numerical simulations. We list the main corner stones of our development.
• Given a generalized potential the coupled-channel reaction amplitudes are defined in terms of a set of non-linear integral equation formulated on contours in the complex plane.
• The analytic structure of the generalized potential in the presence of anomalous threshold effects was clarified. The key observation is the fact that the later must not satisfy the Schwarz reflection principle.
• We confirm previous studies that a physical approach must consider second order terms in the generalized potential. The minimal contributions are identified with the terms that are odd under a Schwarz reflection.
• The non-linear integral equation can be solved numerically by a suitable ansatz with Riemann integrals over real energies only. The specific form of the later was derived for the first time. Explicit expressions for the driving terms were presented for an arbitrary number of channels.
• A schematic 3-channel model was analyzed in the presence of anomalous thresholds.
Along our formal developments all key quantities were illustrated in this model. In particular the reaction amplitdues as well as the phase shifts and in-elasticity parameters were computed and discussed.
Given our framework it is now possible to investigate coupled-channel systems including J P = 1 − and J P = 3 2 + states using realistic interactions. Such systems are notoriously challenging since a plethora of anomalous threshold effects are present. We expect such studies to shed more light on the possible relevance of the hadrogenesis conjecture, which predicts such computations to generate a large part of the hadronic excitation spectrum in QCD.