1 Introduction↩︎

The Green’s function of a system of coupled radial Schrödinger equations is a fundamental object in quantum scattering theory. It governs the propagation of a particle or system through a set of coupled channels and appears naturally whenever one reduces a many-channel problem to an effective single-channel description. In nuclear, atomic, and molecular physics, coupled-channel Green’s functions arise in the construction of optical potentials, the evaluation of transition amplitudes, and the treatment of resonance phenomena. Despite their central role, the construction of the full matrix Green’s function for \(N\) coupled channels, and the proof that the standard construction is the unique solution satisfying the correct boundary conditions, has received surprisingly little attention in the literature.

For a single uncoupled channel, the radial Green’s function is a textbook result: it is built from two linearly independent solutions (regular and irregular) of the homogeneous equation, joined at the source point with the appropriate derivative discontinuity, and normalized by the Wronskian of the two solutions. The extension to \(N\) coupled channels requires \(2N\) linearly independent solution vectors organized into two \(N \times N\) fundamental solution matrices, and the scalar Wronskian is replaced by a Wronskian matrix. This multi-channel construction was developed in the context of atomic and molecular physics by Levine and Soven [1] and further discussed by Levine [2]. Those references constructed a coupled-channel Green’s function by first building the regular and irregular solutions through eigenchannel (\(K\)-matrix) decomposition, and then verified that the resulting bilinear formula satisfies the correct Wronskian relations and boundary conditions. The subsequent Brief Report by Levine [2] addressed the important practical problem of numerical stabilization during integration. The Levine–Soven analysis establishes that the bilinear formula is a Green’s function. It does not address whether it is the only one, and the structural properties of the Wronskian matrix and the reciprocity of the Green’s matrix are imposed and verified rather than derived from the defining equation itself. In the present work we close this gap: starting from the defining equation and the regular/outgoing boundary conditions, we derive the bilinear formula directly and prove uniqueness through the symplectic structure of the \(2N\)-dimensional phase space of the coupled equations. The constancy and diagonal form of the Wronskian matrix, and the reciprocity of the Green’s matrix under the exchange \((R,\gamma)\leftrightarrow(R',\gamma')\), then emerge as necessary consequences of the construction rather than as properties to be checked case by case. An analogous connection between transfer matrices and matrix Green’s functions for systems of coupled second-order differential equations has been developed within the surface Green function matching (SGFM) framework by Rodriguez-Coppola et al. [3] and treated systematically in the monograph by Pérez-Álvarez and García-Moliner [4], where the symplectic character of the transfer matrix and its relationship to the matrix Green function are central themes. That framework, however, addresses boundary-value matching problems in slab geometries rather than the radial scattering problem with regular and outgoing-wave boundary conditions relevant for coupled-channel reactions.

The coupled-channel Green’s function finds important applications across many areas of physics. Wherever the Feshbach projection formalism [5], [6] is used to reduce a multichannel problem to an effective single-channel description, the Green’s function of the eliminated channels serves as the propagator that generates the nonlocal, energy-dependent polarization potential [7]–[9]. A particularly demanding application is the scattering of weakly bound nuclei [10], [11], where the continuum-discretized coupled-channels (CDCC) method [12]–[15] yields large systems of coupled equations with strong continuum-continuum couplings. The resulting dynamical polarization potential (DPP) has been shown to produce long-range attraction and absorption for halo projectiles [16], [17], and its microscopic construction has acquired renewed importance at rare-isotope beam facilities [18]. Most previous studies have adopted a weak-coupling approximation that neglects the continuum-continuum couplings and renders the Green’s matrix diagonal [19]–[21], but recent work has shown that the off-diagonal elements generate coherent interference contributions to the total absorption [22]. We emphasize that this simplification is specific to the Green’s-matrix construction of the DPP and is not a feature of loosely-bound-nucleus reaction theory in general. The continuum-discretized coupled-channels method itself retains the full continuum-continuum coupling, as do the adiabatic approximation for deuteron-induced reactions [23] and its halo-nucleus extension [24] through closure on the projectile internal Hamiltonian, the Weinberg state expansion that systematically improves the latter [25]–[27], and the Faddeev treatment of three-body reaction dynamics [28]–[30]. The weak-coupling simplification does not enter the solution of the three-body dynamics; it is introduced only when that dynamics is reduced to an effective two-body interaction through a diagonal Green’s matrix. The present work removes this approximation by constructing the exact coupled-channel Green’s matrix for the finite CDCC model space, retaining all continuum-continuum couplings among the included bins together with the exact channel-energy structure, from which the effective two-body interaction follows; this full coupled-channel Green’s function is precisely the object whose construction and uniqueness we establish in this work. Analogous coupled-channel Green’s functions appear in electron-atom and electron-molecule collisions [1], multichannel ultracold atomic scattering with multiple hyperfine channels, and coupled vibrational-rotational dynamics.

The central results of this paper are as follows. Starting from the defining equation and boundary conditions, we derive the coupled-channel Green’s matrix \(g_{\gamma\gamma'}(R,R')\) and prove that the bilinear construction in terms of regular and outgoing solution matrices is the unique solution. We show that the associated Wronskian matrix is necessarily diagonal with elements \(W_n = -k_n\) and independent of the radial coordinate, provided that the coupling potential matrix is symmetric, and that the symmetry of the Green’s matrix under channel and coordinate exchange is a consequence of the symplectic structure of the \(2N\)-dimensional phase space. The vanishing of the self-Wronskians of the regular and irregular solution matrices plays a central role in this proof by eliminating cross-terms in the matching conditions. The construction and uniqueness proof apply to any system of \(N\) coupled radial Schrödinger equations with symmetric coupling potentials and open channels, regardless of the physical origin of the coupling. As an illustrative application, we show how the Green’s matrix enters the nonlocal DPP within the CDCC framework, retaining all continuum-continuum couplings. We also discuss the numerical challenges of inward propagation of irregular solutions and the stabilization techniques required to maintain accuracy.

The numerical implementation of this formalism and its application to deuteron-induced reactions on \(^{58}\)Ni are presented in a companion paper [31], where the full-coupling effective potential is shown to reproduce CDCC observables exactly, while the weak-coupling approximation exhibits significant deviations.

The remainder of this paper is organized as follows. Section II presents the construction of the coupled-channel Green’s function from general fundamental solution matrices, proves the Wronskian properties, and establishes uniqueness through the symplectic structure of the \(2N\)-dimensional phase space. Section III illustrates the formalism by deriving the nonlocal dynamical polarization potential within the CDCC framework. Section IV provides a summary and outlook.

2 Formalism↩︎

We consider a general system of \(N\) coupled radial Schrödinger equations with symmetric coupling potentials. The Green’s matrix elements \(g_{\gamma\gamma'}(R,R')\) for this system satisfy the defining equation \[\sum_{\beta}\mathcal{L}_{\gamma\beta}(R)\,g_{\beta\gamma'}(R,R')=\delta_{\gamma\gamma'}\,\delta(R-R'), \label{eq46Gdef}\tag{1}\] where the differential operator is \[\mathcal{L}_{\gamma\beta}(R) = \left(E_\gamma - T_\gamma\right)\delta_{\gamma\beta} - V_{\gamma\beta}(R), \label{eq46Lop}\tag{2}\] with \[T_\gamma = -\frac{\hbar^2}{2\mu}\left(\frac{d^2}{dR^2} - \frac{L_\gamma(L_\gamma+1)}{R^2}\right)\] being the radial kinetic energy operator for channel \(\gamma\), \(\mu\) the reduced mass, \(E_\gamma\) the channel energy, \(L_\gamma\) the orbital angular momentum quantum number, and \(V_{\gamma\beta}(R)\) the coupling potential matrix element connecting channels \(\gamma\) and \(\beta\). The coupled radial equations for the homogeneous problem read \[(E_\gamma - T_\gamma) f_\gamma(R) - \sum_{\gamma'} V_{\gamma\gamma'}(R)\, f_{\gamma'}(R) = 0. \label{eq46cc}\tag{3}\] For a system with \(N\) coupled channels, this equation admits \(2N\) linearly independent solutions, which can be organized into two fundamental sets: \(N\) regular solutions and \(N\) irregular (outgoing) solutions.

The \(N\) regular solutions are defined by their behavior at small radial distance. Each regular solution vector \(u^n\) (\(n = 1, \ldots, N\)) satisfies the boundary condition \(u_\gamma^n(R \to 0) \propto R^{L_\gamma + 1}\,\delta_{\gamma n}\), meaning that at the origin, only the \(n\)-th channel component is nonzero, and it vanishes as \(R^{L_\gamma+1}\) in accordance with the centrifugal barrier. The \(N\) irregular solutions are defined by their asymptotic behavior at large distance. Each irregular solution vector \(h^n\) satisfies \(h_\gamma^n(R \to \infty) \to \delta_{\gamma n}\,H^{(+)}_{L_\gamma}(k_\gamma R)\) in open channels, where \(H^{(+)}_{L_\gamma}\) is the outgoing Coulomb–Hankel function (or Riccati–Hankel function in the absence of the Coulomb interaction) and \(k_\gamma = \sqrt{2\mu E_\gamma}/\hbar\) is the channel wave number. These irregular solutions are obtained by inward integration from a large matching radius, which poses specific numerical challenges that will be discussed below. These \(2N\) solution vectors can be organized into two \(N \times N\) solution matrices \[\begin{align} \mathbf{U} &= \begin{bmatrix} u^1 & \cdots & u^N \end{bmatrix} = \begin{bmatrix} u_1^1 & \cdots & u_1^N \\ \vdots & \ddots & \vdots \\ u_N^1 & \cdots & u_N^N \end{bmatrix}, \\ \mathbf{H} &= \begin{bmatrix} h^1 & \cdots & h^N \end{bmatrix} = \begin{bmatrix} h_1^1 & \cdots & h_1^N \\ \vdots & \ddots & \vdots \\ h_N^1 & \cdots & h_N^N \end{bmatrix}, \end{align} \label{eq46UH}\tag{4}\] where the rows are labeled by the channel index \(\gamma\) and the columns by the solution index \(n\).

We will show below that the coupled-channel Green’s function is uniquely determined by the defining equation and boundary conditions, and that it takes the form (cf.Refs. [1], [2]) \[G(R,R') = \frac{2\mu}{\hbar^2} \left\{\begin{align} \mathbf{U}(R)\, \mathbf{W}^{-1}\, \mathbf{H}^T(R') & \quad R < R', \\ \mathbf{H}(R)\, \mathbf{W}^{-1}\, \mathbf{U}^T(R') & \quad R > R', \end{align}\right. \label{eq46g95full}\tag{5}\] where \(\mathbf{W}\) is the Wronskian matrix defined by \[\mathbf{W}(R) = \mathbf{U}^T(R)\, \frac{d}{dR}\mathbf{H}(R) - \left[\frac{d}{dR}\mathbf{U}^T(R)\right] \mathbf{H}(R). \label{eq46Wdef}\tag{6}\] The structure of this expression merits careful explanation. For \(R < R'\), the Green’s function is built from the regular solutions at the smaller coordinate and the irregular solutions at the larger coordinate, ensuring that \(G(R,R')\) satisfies the correct boundary conditions as \(R \to 0\) and \(R' \to \infty\). The Wronskian matrix \(\mathbf{W}\) encodes the linear independence of the two fundamental sets and plays the role of a normalization factor. Its properties, constancy and diagonal structure, are crucial for the practical utility of Eq. 5 and are established in what follows.

We now demonstrate that, when the coupling potential matrix is symmetric (\(V_{\gamma\gamma'} = V_{\gamma'\gamma}\)), the Wronskian matrix \(\mathbf{W}\) is independent of the radial coordinate \(R\) and is, moreover, a diagonal matrix. The element \(W_{nm}\) of the Wronskian matrix is given explicitly by \[W_{nm} = \sum_{\gamma=1}^N \left( u_\gamma^n\, {h_\gamma^{m}}' - {u_\gamma^{n}}'\, h_\gamma^m \right), \label{eq46Wnm}\tag{7}\] where primes denote derivatives with respect to \(R\). Taking the derivative of \(W_{nm}\) with respect to \(R\) yields \[\frac{d}{dR} W_{nm} = \sum_{\gamma=1}^N \left( u_\gamma^n\, {h_\gamma^{m}}'' - {u_\gamma^{n}}''\, h_\gamma^m \right). \label{eq46dWnm}\tag{8}\] The second derivatives appearing here can be expressed through the coupled-channel equations 3 . After multiplying by \(2\mu/\hbar^2\) and rearranging, any solution \(f_\gamma\) of the coupled equations satisfies \[\begin{align} {f_\gamma}'' &= \sum_{\gamma'}\mathcal{M}_{\gamma\gamma'}(R)\,f_{\gamma'},\\ \mathcal{M}_{\gamma\gamma'}(R) &= \left[\frac{L_\gamma(L_\gamma+1)}{R^2}-\frac{2\mu}{\hbar^2}E_\gamma\right]\delta_{\gamma\gamma'} + \frac{2\mu}{\hbar^2}V_{\gamma\gamma'}(R). \end{align} \label{eq46M}\tag{9}\] Substituting the second-derivative expressions for \(u_\gamma^n\) and \(h_\gamma^m\) into Eq. 8 , we obtain \[\frac{d}{dR} W_{nm} = \sum_{\gamma=1}^N\sum_{\gamma'=1}^N \left(u_\gamma^n\,\mathcal{M}_{\gamma\gamma'}\,h_{\gamma'}^m - h_\gamma^m\,\mathcal{M}_{\gamma\gamma'}\,u_{\gamma'}^n\right). \label{eq46dWnm2}\tag{10}\] The symmetry of the coupling potential, \(V_{\gamma\gamma'} = V_{\gamma'\gamma}\), implies that \(\mathcal{M}_{\gamma\gamma'} = \mathcal{M}_{\gamma'\gamma}\). Consequently, the two terms in the double sum are related by the interchange \(\gamma \leftrightarrow \gamma'\), and the entire expression vanishes identically: \(dW_{nm}/dR = 0\). This establishes that each element of the Wronskian matrix is a constant, independent of the radial coordinate \(R\), and its value is therefore determined once and for all by the boundary conditions.

The constant values of \(W_{nm}\) are most conveniently evaluated from the asymptotic forms of the solutions. In the limit \(R \to \infty\), the regular and irregular solutions approach \[\begin{align} u_\gamma^n &\rightarrow \frac{i}{2} \left( H^{(-)}_{L_\gamma}(k_\gamma R)\, \delta_{n\gamma} - S_{n\gamma}\, H^{(+)}_{L_\gamma}(k_\gamma R) \right), \\ h_\gamma^n &\rightarrow \delta_{n\gamma}\, H^{(+)}_{L_\gamma}(k_\gamma R), \end{align} \label{eq46asympt}\tag{11}\] where \(H^{(\pm)}_{L}\) are the incoming and outgoing Coulomb–Hankel functions and \(S_{n\gamma}\) are the elements of the scattering \(S\)-matrix. For the off-diagonal elements with \(n \neq m\), the asymptotic boundary condition \(h_\gamma^m \to \delta_{\gamma m}\,H^{(+)}_{L_\gamma}(k_\gamma R)\) ensures that only the \(\gamma = m\) term survives in the sum over channels. However, this surviving term involves the Wronskian of \(H^{(+)}\) with itself, \(H^{(+)}{H^{(+)}}' - {H^{(+)}}'H^{(+)} = 0\), which vanishes identically. Therefore \(W_{nm} = 0\) for all \(n \neq m\), and the Wronskian matrix is diagonal.

For the diagonal elements (\(n = m\)), only the \(\gamma = n\) term contributes, again because \(h_\gamma^n\) vanishes asymptotically for \(\gamma \neq n\). Evaluating the remaining term using Eq. 11 gives \[\begin{align} W_n &= u_n^n\, {h_n^{n}}' - {u_n^{n}}'\, h_n^n \\ &= \frac{i}{2} \Big[ \left( H^{(-)}(k_n R) - S_{nn}\, H^{(+)}(k_n R) \right) {H^{(+)}}'(k_n R) \\ &\quad - \left( {H^{(-)}}'(k_n R) - S_{nn}\, {H^{(+)}}'(k_n R) \right) H^{(+)}(k_n R) \Big] \\ &= \frac{i}{2} \left[ H^{(-)}(k_n R)\, {H^{(+)}}'(k_n R) - {H^{(-)}}'(k_n R)\, H^{(+)}(k_n R) \right] \\ &= -k_n, \end{align} \label{eq46Wn}\tag{12}\] where the terms proportional to \(S_{nn}\) cancel, and the final equality follows from the standard Wronskian relation of the Coulomb–Hankel functions under the Riccati normalization convention adopted throughout this work. Different conventions for the Hankel functions modify only the overall constant prefactor and sign of \(W_n\), without affecting the physical content of the Green’s function. Thus the Wronskian matrix takes the simple diagonal form \(\mathbf{W} = \text{diag}(-k_1, -k_2, \ldots, -k_N)\).

With the diagonal Wronskian established, the Green’s function of Eq. 5 can be written in a form that makes the role of each linearly independent solution explicit. Expanding the matrix products, the individual matrix elements become \[g_{\gamma\gamma'}(R, R') = \frac{2\mu}{\hbar^2} \left\{ \begin{align} \sum_n \frac{u_\gamma^n(R)\, h_{\gamma'}^n(R')}{W_n} & \quad R < R', \\ \sum_n \frac{h_\gamma^n(R)\, u_{\gamma'}^n(R')}{W_n} & \quad R > R', \end{align}\right. \label{eq46g95elements}\tag{13}\] where the sum runs over the \(N\) solution indices. Each term in this sum represents the contribution of one linearly independent solution pair to the propagation from channel \(\gamma'\) at \(R'\) to channel \(\gamma\) at \(R\). For symmetric couplings, the Green’s matrix also satisfies the reciprocity relation \(g_{\gamma\gamma'}(R,R') = g_{\gamma'\gamma}(R',R)\), which provides a useful numerical consistency check during implementation.

We now derive the Green’s matrix of Eq. 5 directly from the defining equation 1 and the boundary conditions, and prove that it is the unique solution. Unlike the original construction [1], [2], which arrived at the bilinear formula through an eigenchannel-based procedure and then verified its correctness, the derivation below obtains the formula as the only possible outcome of the boundary conditions and matching requirements. As a byproduct, the Wronskian constancy \(W_n = -k_n\) and the reciprocity of the Green’s matrix emerge as necessary consequences, rather than properties that need to be checked separately. The proof proceeds in three steps: establishing the self-Wronskian identities, demonstrating the uniqueness of the functional form, and solving the matching conditions at the source point.

For any two \(N \times N\) solution matrices \(\mathbf{F}\) and \(\mathbf{G}\) of the coupled radial equations 3 , the bilinear concomitant \(\mathcal{W}[\mathbf{F},\mathbf{G}] \equiv \mathbf{F}^T\mathbf{G}' - \mathbf{F}'^T\mathbf{G}\) is independent of \(R\) when \(V_{\gamma\gamma'} = V_{\gamma'\gamma}\), as demonstrated above. Applying this result to the self-Wronskians of the fundamental solution matrices, we find that for the regular solutions, the small-\(R\) boundary condition \(u_\gamma^n(R \to 0) \propto R^{L_\gamma+1}\delta_{\gamma n}\) implies that every element of \(\mathbf{U}^T\mathbf{U}' - \mathbf{U}'^T\mathbf{U}\) vanishes at the origin, since the product of any two regular solutions both vanish as \(R \to 0\). By constancy, this gives \[\mathcal{W}[\mathbf{U},\mathbf{U}] = 0. \label{eq46WUU}\tag{14}\] Similarly, the large-\(R\) boundary condition \(h_\gamma^n(R \to \infty) \to \delta_{\gamma n}H^{(+)}_{L_\gamma}(k_\gamma R)\), together with the vanishing of the single-function Wronskian \(H^{(+)}{H^{(+)}}' - {H^{(+)}}'H^{(+)} = 0\) and the column-wise orthogonality of the asymptotic conditions, yields \[\mathcal{W}[\mathbf{H},\mathbf{H}] = 0. \label{eq46WHH}\tag{15}\]

The Green’s matrix \(G(R,R')\) must be regular as \(R \to 0\) (to avoid unphysical singularities at the origin) and purely outgoing as \(R \to \infty\) (to enforce causality). For a fixed source point \(R'\), the most general matrix function with these properties can be expressed as \[G(R,R') = \begin{cases} \mathbf{U}(R)\,\mathbf{A}(R'), & R < R',\\ \mathbf{H}(R)\,\mathbf{B}(R'), & R > R', \end{cases} \label{eq46Gform}\tag{16}\] where \(\mathbf{A}\) and \(\mathbf{B}\) are \(N \times N\) coefficient matrices to be determined. This is the only form compatible with the boundary conditions, since any regular solution of the homogeneous equation is a linear combination of the columns of \(\mathbf{U}\), and any outgoing solution is a linear combination of the columns of \(\mathbf{H}\). No ansatz is involved in writing Eq. 16 ; it follows uniquely from the completeness of the fundamental solution sets.

The coefficient matrices are determined by the matching conditions at the source point \(R = R'\). Continuity of the Green’s function requires \[\mathbf{U}(R')\,\mathbf{A}(R') = \mathbf{H}(R')\,\mathbf{B}(R'), \label{eq46cont}\tag{17}\] while integrating the defining equation 1 across the delta-function source at \(R = R'\) yields the derivative jump condition \[\mathbf{H}'(R')\,\mathbf{B}(R') - \mathbf{U}'(R')\,\mathbf{A}(R') = \frac{2\mu}{\hbar^2}\,\mathbf{I}. \label{eq46jump}\tag{18}\] We now solve for \(\mathbf{B}\). Left-multiplying the continuity condition 17 by \(\mathbf{U}'^T\) gives \[\mathbf{U}'^T \mathbf{U}\,\mathbf{A} = \mathbf{U}'^T \mathbf{H}\,\mathbf{B}, \label{eq46step1}\tag{19}\] while left-multiplying the jump condition 18 by \(\mathbf{U}^T\) gives \[\mathbf{U}^T \mathbf{H}'\,\mathbf{B} - \mathbf{U}^T \mathbf{U}'\,\mathbf{A} = \frac{2\mu}{\hbar^2}\,\mathbf{U}^T. \label{eq46step2}\tag{20}\] Subtracting Eq. 19 from Eq. 20 , we collect the terms containing \(\mathbf{A}\) and \(\mathbf{B}\) separately: \[\begin{align} &\underbrace{(\mathbf{U}^T \mathbf{H}' - \mathbf{U}'^T \mathbf{H})}_{\displaystyle =\,\mathbf{W}}\,\mathbf{B}\\ &\quad -\;\underbrace{(\mathbf{U}^T \mathbf{U}' - \mathbf{U}'^T \mathbf{U})}_{\displaystyle =\,\mathcal{W}[\mathbf{U},\mathbf{U}]\,=\,0}\,\mathbf{A} = \frac{2\mu}{\hbar^2}\,\mathbf{U}^T. \end{align}\] The first bracket is precisely the Wronskian matrix \(\mathbf{W}\) as defined in Eq. 6 . The second bracket is the self-Wronskian \(\mathcal{W}[\mathbf{U},\mathbf{U}]\), which vanishes identically by Eq. 14 . This eliminates \(\mathbf{A}\) entirely, leaving \[\mathbf{W}\,\mathbf{B}(R') = \frac{2\mu}{\hbar^2}\,\mathbf{U}^T(R'),\] from which \[\mathbf{B}(R') = \frac{2\mu}{\hbar^2}\,\mathbf{W}^{-1}\,\mathbf{U}^T(R'). \label{eq46B}\tag{21}\]

The derivation of \(\mathbf{A}\) follows the same strategy with a different choice of multipliers. Left-multiplying Eq. 17 by \(\mathbf{H}'^T\) gives \[\mathbf{H}'^T \mathbf{U}\,\mathbf{A} = \mathbf{H}'^T \mathbf{H}\,\mathbf{B}, \label{eq46step3}\tag{22}\] and left-multiplying Eq. 18 by \(\mathbf{H}^T\) gives \[\mathbf{H}^T \mathbf{H}'\,\mathbf{B} - \mathbf{H}^T \mathbf{U}'\,\mathbf{A} = \frac{2\mu}{\hbar^2}\,\mathbf{H}^T. \label{eq46step4}\tag{23}\] Subtracting Eq. 22 from Eq. 23 : \[\begin{align} &\underbrace{(\mathbf{H}^T \mathbf{H}' - \mathbf{H}'^T \mathbf{H})}_{\displaystyle =\,\mathcal{W}[\mathbf{H},\mathbf{H}]\,=\,0}\,\mathbf{B}\\ &\quad -\;\underbrace{(\mathbf{H}^T \mathbf{U}' - \mathbf{H}'^T \mathbf{U})}_{\displaystyle =\,-\mathbf{W}^T}\,\mathbf{A} = \frac{2\mu}{\hbar^2}\,\mathbf{H}^T. \end{align}\] The first bracket is the self-Wronskian \(\mathcal{W}[\mathbf{H},\mathbf{H}] = 0\) by Eq. 15 , which eliminates \(\mathbf{B}\). For the second bracket, we note that \(\mathbf{H}^T\mathbf{U}' - \mathbf{H}'^T\mathbf{U}\) is the negative transpose of the Wronskian: since \(\mathbf{W} = \mathbf{U}^T\mathbf{H}' - \mathbf{U}'^T\mathbf{H}\), taking the transpose gives \(\mathbf{W}^T = \mathbf{H}'^T\mathbf{U} - \mathbf{H}^T\mathbf{U}'\), so that \(\mathbf{H}^T\mathbf{U}' - \mathbf{H}'^T\mathbf{U} = -\mathbf{W}^T\). The equation therefore reduces to \[\mathbf{W}^T\,\mathbf{A}(R') = \frac{2\mu}{\hbar^2}\,\mathbf{H}^T(R'),\] from which \[\mathbf{A}(R') = \frac{2\mu}{\hbar^2}\,\mathbf{W}^{-T}\,\mathbf{H}^T(R'). \label{eq46A}\tag{24}\] Since \(\mathbf{W}\) is diagonal for symmetric couplings, we have \(\mathbf{W}^{-T} = \mathbf{W}^{-1}\). Substituting the expressions 21 and 24 for \(\mathbf{B}\) and \(\mathbf{A}\) back into Eq. 16 reproduces exactly the construction of Eq. 5 .

It is instructive to verify that the resulting Green’s matrix is indeed self-consistent by checking that the continuity condition 17 and the jump condition 18 are both satisfied. Substituting Eqs. 21 and 24 into the continuity condition 17 , we need \(\mathbf{U}\,\mathbf{A} = \mathbf{H}\,\mathbf{B}\), which after inserting the explicit expressions becomes \[\mathbf{U}\,\mathbf{W}^{-1}\,\mathbf{H}^T = \mathbf{H}\,\mathbf{W}^{-1}\,\mathbf{U}^T, \label{eq46symmetry}\tag{25}\] where the argument \(R'\) and the common factor \(\frac{2\mu}{\hbar^2}\) have been suppressed for brevity. This identity states that the matrix \(\mathbf{C}\equiv\mathbf{U}\mathbf{W}^{-1}\mathbf{H}^T\) is equal to its own transpose \(\mathbf{C}^T=\mathbf{H}\mathbf{W}^{-1}\mathbf{U}^T\). To see why this is nontrivial, note that the element \((\gamma,\gamma')\) of \(\mathbf{C}\) and \(\mathbf{C}^T\) read \[C_{\gamma\gamma'}=\sum_n\frac{u_\gamma^n\,h_{\gamma'}^n}{W_n},\qquad C^T_{\gamma\gamma'}=C_{\gamma'\gamma}=\sum_n\frac{u_{\gamma'}^n\,h_\gamma^n}{W_n}.\] These two sums are manifestly different term by term: for each \(n\), the outer product \(h^n u^{nT}\) (with elements \(h_\gamma^n u_{\gamma'}^n\)) is not the same matrix as \(u^n h^{nT}\) (with elements \(u_\gamma^n h_{\gamma'}^n\)), yet the claim is that after summing over all \(n\) with weights \(1/W_n\), the result is symmetric. This is a completeness property of the full set of \(2N\) fundamental solutions, and we now prove it using the symplectic structure of the solution space.

We assemble the \(2N\) fundamental solutions and their derivatives into the \(2N\times 2N\) phase-space matrix \[\boldsymbol{\Phi}(R) = \begin{pmatrix} \mathbf{U}(R) & \mathbf{H}(R) \\ \mathbf{U}'(R) & \mathbf{H}'(R) \end{pmatrix}, \label{eq46Phi}\tag{26}\] whose upper-left and upper-right \(N\times N\) blocks contain the solution values, and whose lower blocks contain their derivatives. The key observation is that the coupled radial equations preserve a symplectic structure encoded in the \(2N\times 2N\) matrix \(\mathbf{J}=\bigl(\begin{smallmatrix}0 & \mathbf{I}_N \\ -\mathbf{I}_N & 0\end{smallmatrix}\bigr)\). The generalized Wronskian \(\boldsymbol{\Phi}^T\mathbf{J}\boldsymbol{\Phi}\) is constant in \(R\) (as a consequence of the symmetry \(V_{\gamma\gamma'}=V_{\gamma'\gamma}\), by the same argument that proved constancy of \(\mathbf{W}\)). Computing this product block by block, \[\begin{align} \boldsymbol{\Phi}^T \mathbf{J}\,\boldsymbol{\Phi} &= \begin{pmatrix} \mathbf{U}^T\mathbf{U}'-\mathbf{U}'^T\mathbf{U} & \mathbf{U}^T\mathbf{H}'-\mathbf{U}'^T\mathbf{H} \\ \mathbf{H}^T\mathbf{U}'-\mathbf{H}'^T\mathbf{U} & \mathbf{H}^T\mathbf{H}'-\mathbf{H}'^T\mathbf{H} \end{pmatrix}\\ &= \begin{pmatrix} \mathcal{W}[\mathbf{U},\mathbf{U}] & \mathbf{W} \\ -\mathbf{W}^T & \mathcal{W}[\mathbf{H},\mathbf{H}] \end{pmatrix}, \end{align} \label{eq46PhiJPhi}\tag{27}\] where we have identified the four blocks: the \((1,1)\) block is \(\mathcal{W}[\mathbf{U},\mathbf{U}]=0\) by Eq. 14 , the \((1,2)\) block is \(\mathbf{W}\) by its definition 6 , the \((2,2)\) block is \(\mathcal{W}[\mathbf{H},\mathbf{H}]=0\) by Eq. 15 , and the \((2,1)\) block is \(\mathbf{H}^T\mathbf{U}'-\mathbf{H}'^T\mathbf{U}=-\mathbf{W}^T=-\mathbf{W}\) (the last equality holds because \(\mathbf{W}\) is diagonal). Thus \[\boldsymbol{\Phi}^T \mathbf{J}\,\boldsymbol{\Phi} = \begin{pmatrix} 0 & \mathbf{W} \\ -\mathbf{W} & 0 \end{pmatrix} \equiv \mathbf{S}. \label{eq46S}\tag{28}\] Since \(W_n=-k_n\neq 0\) for all channels, \(\mathbf{W}\) is invertible, and so is \(\mathbf{S}\) with inverse \[\mathbf{S}^{-1} = \begin{pmatrix} 0 & -\mathbf{W}^{-1} \\ \mathbf{W}^{-1} & 0 \end{pmatrix}, \label{eq46Sinv}\tag{29}\] as can be verified directly from \(\mathbf{S}\,\mathbf{S}^{-1}=\mathbf{I}_{2N}\). Since \(\det\mathbf{S}\neq 0\) and \(\boldsymbol{\Phi}^T\mathbf{J}\boldsymbol{\Phi}=\mathbf{S}\), the matrix \(\boldsymbol{\Phi}\) itself must be invertible. From \(\boldsymbol{\Phi}^T\mathbf{J}\boldsymbol{\Phi}=\mathbf{S}\) we obtain \(\boldsymbol{\Phi}^T\mathbf{J}=\mathbf{S}\,\boldsymbol{\Phi}^{-1}\), and therefore \[\boldsymbol{\Phi}^{-1} = \mathbf{S}^{-1}\,\boldsymbol{\Phi}^T\,\mathbf{J}. \label{eq46Phiinv95formula}\tag{30}\] We now evaluate the right-hand side step by step. First, \[\boldsymbol{\Phi}^T\,\mathbf{J} = \begin{pmatrix} \mathbf{U}^T & \mathbf{U}'^T \\ \mathbf{H}^T & \mathbf{H}'^T \end{pmatrix} \begin{pmatrix} 0 & \mathbf{I} \\ -\mathbf{I} & 0 \end{pmatrix} = \begin{pmatrix} -\mathbf{U}'^T & \mathbf{U}^T \\ -\mathbf{H}'^T & \mathbf{H}^T \end{pmatrix}. \label{eq46PhiTJ}\tag{31}\] Then, left-multiplying by \(\mathbf{S}^{-1}\), \[\begin{align} \boldsymbol{\Phi}^{-1} &= \begin{pmatrix} 0 & -\mathbf{W}^{-1} \\ \mathbf{W}^{-1} & 0 \end{pmatrix} \begin{pmatrix} -\mathbf{U}'^T & \mathbf{U}^T \\ -\mathbf{H}'^T & \mathbf{H}^T \end{pmatrix}\\ &= \begin{pmatrix} \mathbf{W}^{-1}\mathbf{H}'^T & -\mathbf{W}^{-1}\mathbf{H}^T \\ -\mathbf{W}^{-1}\mathbf{U}'^T & \mathbf{W}^{-1}\mathbf{U}^T \end{pmatrix}. \end{align} \label{eq46Phiinv}\tag{32}\] We can now extract the desired identities from the completeness relation \(\boldsymbol{\Phi}\,\boldsymbol{\Phi}^{-1}=\mathbf{I}_{2N}\). Writing out the block matrix product explicitly, \[\begin{pmatrix} \mathbf{U} & \mathbf{H} \\ \mathbf{U}' & \mathbf{H}' \end{pmatrix} \begin{pmatrix} \mathbf{W}^{-1}\mathbf{H}'^T & -\mathbf{W}^{-1}\mathbf{H}^T \\ -\mathbf{W}^{-1}\mathbf{U}'^T & \mathbf{W}^{-1}\mathbf{U}^T \end{pmatrix} =\begin{pmatrix} \mathbf{I} & 0 \\ 0 & \mathbf{I} \end{pmatrix}. \label{eq46completeness}\tag{33}\] Each of the four \(N\times N\) blocks in this equation yields an independent identity. The upper-right \((1,2)\) block gives \[\begin{align} &\mathbf{U}\bigl(-\mathbf{W}^{-1}\mathbf{H}^T\bigr) +\mathbf{H}\bigl(\mathbf{W}^{-1}\mathbf{U}^T\bigr)\\ &\quad= -\mathbf{U}\,\mathbf{W}^{-1}\,\mathbf{H}^T + \mathbf{H}\,\mathbf{W}^{-1}\,\mathbf{U}^T = 0, \end{align} \label{eq46cont95check}\tag{34}\] which is precisely the continuity condition 25 . This identity is the completeness relation that guarantees \(\mathbf{C}=\mathbf{C}^T\), despite the individual outer products being non-symmetric. The lower-right \((2,2)\) block gives \[\begin{align} &\mathbf{U}'\bigl(-\mathbf{W}^{-1}\mathbf{H}^T\bigr) +\mathbf{H}'\bigl(\mathbf{W}^{-1}\mathbf{U}^T\bigr)\\ &\quad= \mathbf{H}'\,\mathbf{W}^{-1}\,\mathbf{U}^T - \mathbf{U}'\,\mathbf{W}^{-1}\,\mathbf{H}^T = \mathbf{I}. \end{align} \label{eq46jump95check}\tag{35}\] Substituting the explicit forms of \(\mathbf{A}\) and \(\mathbf{B}\) into the jump condition 18 , the left-hand side reads \[\begin{align} \mathbf{H}'\mathbf{B}-\mathbf{U}'\mathbf{A} &=\frac{2\mu}{\hbar^2}\left(\mathbf{H}'\mathbf{W}^{-1}\mathbf{U}^T-\mathbf{U}'\mathbf{W}^{-1}\mathbf{H}^T\right)\\ &=\frac{2\mu}{\hbar^2}\,\mathbf{I}, \end{align}\] where the second equality uses Eq. 35 . This matches the right-hand side of Eq. 18 exactly. The remaining two blocks yield analogous identities involving derivatives: the upper-left \((1,1)\) block gives \[\mathbf{U}\,\mathbf{W}^{-1}\,\mathbf{H}'^T - \mathbf{H}\,\mathbf{W}^{-1}\,\mathbf{U}'^T = \mathbf{I},\] and the lower-left \((2,1)\) block gives \[\mathbf{U}'\,\mathbf{W}^{-1}\,\mathbf{H}'^T - \mathbf{H}'\,\mathbf{W}^{-1}\,\mathbf{U}'^T = 0.\] All four blocks of the completeness relation are satisfied simultaneously, confirming that Eq. 5 is the unique Green’s matrix satisfying the correct boundary conditions at the origin and at infinity, continuity at the source point, and the prescribed derivative discontinuity.

We emphasize that the construction and uniqueness proof presented above rely on three conditions: (i) the coupling potential matrix is symmetric, \(V_{\gamma\gamma'} = V_{\gamma'\gamma}\); (ii) all \(N\) channels are open, so that \(k_\gamma > 0\) and the Wronskian matrix \(\mathbf{W}\) is invertible; and (iii) the regular and outgoing boundary conditions are well defined. No assumption has been made about the physical origin of the coupled channels. The formalism therefore applies to any system of \(N\) coupled radial Schrödinger equations satisfying these conditions, including coupled-channels descriptions of inelastic nuclear scattering, electron-atom and electron-molecule collisions, ultracold atomic scattering with multiple hyperfine channels, and coupled vibrational-rotational dynamics. When closed channels are present (\(k_\gamma^2 < 0\)), the outgoing boundary condition is replaced by exponential decay, and the corresponding Wronskian elements acquire a modified form; the extension of the present proof to mixed open-closed systems requires a separate treatment of the closed-channel boundary conditions.

A practical issue that must be addressed in evaluating Eq. 13 is the numerical stability of the irregular solutions during inward propagation. As discussed by Levine and Soven [1] and Levine [2], the irregular solution matrix \(\mathbf{H}\) is obtained by integrating inward from the asymptotic region, where the column-wise boundary conditions are imposed. During this inward propagation, channels with different orbital angular momenta exhibit vastly different growth rates: in the approximately decoupled small-\(R\) limit, the irregular solution in channel \(\gamma\) behaves as \(h_\gamma^n(R) \sim R^{-L_\gamma}\), while the regular solution scales as \(u_\gamma^n(R) \sim R^{L_\gamma+1}\). As \(L_\gamma\) increases, the dynamic range among channels grows rapidly, and the faster-growing components can numerically overwhelm the slower ones, leading to a loss of linear independence among the columns of \(\mathbf{H}\). This contamination propagates into the off-diagonal elements of the Green’s function and manifests as a drift of the computed Wronskian from its exact constant value. Reference [2] points out that such deviations can grow approximately exponentially if stabilization is not enforced during the propagation.

To maintain numerical accuracy, practical implementations should periodically re-stabilize or re-orthogonalize the irregular solution matrix during the inward integration, using techniques such as modified Gram–Schmidt orthogonalization or related procedures. The constancy of the Wronskian matrix provides a sensitive diagnostic: any deviation of the computed \(W_n\) from \(-k_n\) signals a loss of numerical precision. Even with careful stabilization, the Green’s function construction encounters a fundamental limitation at sufficiently high partial waves, where the inward integration of irregular solutions becomes intrinsically ill-conditioned. In applications, one should verify convergence with respect to \(L_{\max}\) and, when necessary, treat very high-\(L\) contributions with asymptotic or perturbative approximations.

3 Application to breakup dynamical polarization potentials↩︎

As an illustrative application, we show how the coupled-channel Green’s function enters the construction of the dynamical polarization potential (DPP) within the continuum-discretized coupled-channels (CDCC) framework for the scattering of weakly bound composite projectiles. Weakly bound nuclei, such as \(^{6}\)He, \(^{11}\)Be, \(^{8}\)B, \(^{6}\)Li, \(^{7}\)Li, and \(^{9}\)Be, have separation energies of the order of one MeV or less, so that coupling between the ground state and the continuum profoundly modifies the reaction dynamics [10], [11]. In the CDCC method [12]–[15], the projectile continuum is discretized into \(N\) square-integrable bin states, yielding a system of \(N\) coupled radial equations of the form 3 for the relative motion between the projectile and target.

The Feshbach projection formalism [5], [6] partitions the model space into the elastic (\(P\)) and continuum (\(Q\)) subspaces and yields an effective single-channel Hamiltonian for the elastic channel: \[H_{\text{eff}} = PHP + PHQ \frac{1}{E + i\varepsilon - QHQ} QHP, \label{eq46Heff}\tag{36}\] where \(P = |\phi_{bx}^0\rangle\langle\phi_{bx}^0|\) projects onto the projectile ground state and \(Q = \sum_{i=1}^{N} |\phi_{bx}^{i}\rangle\langle\phi_{bx}^{i}|\) onto the discretized continuum. The first term contains the bare elastic-channel Hamiltonian, while the second term is the DPP, which is manifestly energy-dependent and nonlocal. The DPP can be expressed in terms of the coupled-channel Green’s matrix constructed in Sec. II as \[\Delta U(R,R') = \sum_{\gamma,\gamma'}^Q U_{0\gamma}(R)\,g_{\gamma\gamma'}(R,R')\,U_{\gamma'0}(R'), \label{eq46DPP}\tag{37}\] where \(U_{0\gamma}(R) = \langle\phi_{bx}^0| U_{bA} + U_{xA} |\phi_{bx}^n\rangle_{\boldsymbol{r}}\) are the coupling form factors between the ground state and the continuum bins, and \(g_{\gamma\gamma'}(R,R')\) are the matrix elements of the \(Q\)-space Green’s operator in the channel basis. The physical interpretation is transparent: the projectile is excited from the ground state into a continuum bin \(\gamma'\) at position \(R'\) through the coupling \(U_{\gamma'0}(R')\), propagates through the coupled continuum channels via \(g_{\gamma\gamma'}(R,R')\), and is de-excited back to the ground state at position \(R\) through \(U_{0\gamma}(R)\). The sum over all intermediate channels encodes all possible excitation and de-excitation pathways. Applying the optical theorem within the \(P\)-space, the imaginary part of the expectation value of \(\Delta U\) with respect to the elastic wave function provides a direct measure of the total reaction flux absorbed into all nonelastic channels.

Most previous studies [19]–[21] have adopted a weak-coupling approximation in which the couplings among the continuum channels are neglected, reducing the Green’s matrix to \(g_{\gamma\gamma'} = g_\gamma \delta_{\gamma\gamma'}\). The DPP then simplifies to a sum of separable terms, \(\Delta U(R,R') = \sum_\gamma U_{0\gamma}(R)\,g_\gamma(R,R')\,U_{\gamma 0}(R')\), which is computationally attractive but neglects the multistep propagation pathways through which continuum-continuum couplings redistribute flux among the breakup channels. Alternative strategies based on trivially equivalent local potentials [32], [33] suppress the intrinsic nonlocality and introduce partial-wave dependence that obscures the microscopic origin of the effective interaction. The full coupled-channel construction of Sec. II retains all such pathways, including the coherent interference contributions highlighted in Ref. [22]. The numerical implementation and application to deuteron scattering on \(^{58}\)Ni are presented in a companion paper [31], where the full-coupling effective potential is shown to reproduce CDCC elastic scattering cross sections exactly, while the weak-coupling and folding-model approaches show significant deviations.

4 Summary and Outlook↩︎

In this work we have presented a rigorous construction and uniqueness proof of the matrix Green’s function for \(N\) coupled radial Schrödinger equations with symmetric coupling potentials. The Green’s matrix was built from two fundamental sets of solutions, regular and outgoing, organized into \(N \times N\) solution matrices. While the bilinear form of the Green’s matrix was previously obtained by Levine and Soven [1], [2] through an eigenchannel construction and verified by substitution, the present work derives it directly from the defining equation and boundary conditions and proves that it is the unique solution. The proof relies on the symplectic structure of the \(2N\)-dimensional phase space. The self-Wronskian identities \(\mathcal{W}[\mathbf{U},\mathbf{U}] = 0\) and \(\mathcal{W}[\mathbf{H},\mathbf{H}] = 0\) play a central role by eliminating the cross-terms in the matching conditions, while the completeness relation \(\boldsymbol{\Phi}\,\boldsymbol{\Phi}^{-1} = \mathbf{I}_{2N}\) guarantees the continuity and symmetry of the Green’s matrix at the source point.

The construction and uniqueness proof hold for any system of coupled radial Schrödinger equations with symmetric coupling potentials and open channels, regardless of the physical context. The same framework is directly applicable to coupled-channels descriptions of inelastic nuclear scattering, effective optical potential constructions via the Feshbach formalism, electron-atom and electron-molecule collisions, and multichannel ultracold scattering problems. Wherever the Feshbach projection formalism is used to reduce a multichannel problem to an effective single-channel description, the coupled-channel Green’s function derived here provides the exact propagator for the eliminated channels.

As a specific application, we showed how the Green’s matrix enters the nonlocal dynamical polarization potential within the CDCC framework, retaining all continuum-continuum couplings. The resulting DPP kernel, expressed as a bilinear sum over the full Green’s matrix elements, naturally incorporates the off-diagonal propagation mechanisms through which continuum-continuum couplings generate additional nonlocality, energy dependence, and coherent interference in the elastic effective interaction.

We also discussed the practical challenges of implementing this construction, particularly the loss of numerical conditioning during inward propagation of irregular solutions at high partial waves. The constancy of the Wronskian provides a built-in diagnostic for monitoring numerical precision.

The numerical implementation of this framework, including its application to deuteron-induced reactions on \(^{58}\)Ni and a quantitative benchmarking of the weak-coupling approximation, is presented in a companion paper [31]. That work demonstrates that the full-coupling effective potential reproduces CDCC elastic scattering cross sections exactly across a range of incident energies, while the weak-coupling and folding-model approaches show significant deviations. Extensions to systematic studies of how breakup-induced nonlocality varies across different projectile-target systems and energy regimes are in progress.

References↩︎