Spectral Methods for Coupled Channels with a Mass Gap

We develop a method to compute the vacuum polarization energy for coupled scalar fields with different masses scattering off a background potential in one space dimension. As an example we consider the vacuum polarization energy of a kink-like soliton built from two real scalar fields with different mass parameters.

relativistic dispersion relations, the momenta are therefore related as We adopt the convention m 1 ≤ m 2 so that k 1 is real for the physical scattering process, and treat k 1 as the independent momentum variable. This choice induces a square root dependence of k 2 on k 1 which may lead to additional branch cuts in the complex k 1 plane. Within the gap, k 2 1 < m 2 2 − m 2 1 , the second momentum variable becomes complex and, at least in the lower complex half plane, analytic properties are most likely lost. Within that gap it is also important to ensure that the wave-function of the heavier particle decreases exponentially towards spatial infinity. In view of these obvious obstacles it is the main objective of the present investigation to establish a comprehensive scattering formalism for the VPE in this case, where we generalize Eq. (1) to complex momenta such that the above required properties of the Jost function F (k 1 ) are maintained.
There are additional technical advantages of analytically continuing to complex momenta. Rotating the momentum integral over the phase shift onto the imaginary axis automatically includes the bound state contribution to the VPE and they need not be explicitly constructed. Nevertheless we will do so in certain cases to verify the analytic properties of F (k 1 ). Additionally, at threshold, k 2 = 0, the phase shift typically exhibits cusps that are difficult to handle numerically. This problem is avoided on the imaginary axis, as is the one of numerically constructing a phase shift as a continuous function of k 1 from the phase of the determinant of the scattering matrix without 2π jumps. Furthermore, as we will discuss below, for m 1 = m 2 a consistent formulation of the no-tadpole renormalization condition is problematic when the VPE is formulated as an integral over real momenta.
The paper is organized as follows. In Section II we explain the problem and relate k 1 and k 2 such that the scattering problem is well-defined along the real k 1 axis and an analytic continuation to the upper half complex k 1 plane is possible. In Section III, we then define the Jost determinant for a spatially symmetric potential and demonstrate the necessary analytic properties of this determinant by numerical simulations. In Section IV we generalize this approach to a potential with mixed reflection symmetry similar to parity in the Dirac theory. In Section V we formulate the VPE in terms of an integral of the Jost determinant along the positive imaginary momentum axis. In Section VI we apply this spectral method to the soliton model with two real scalar fields proposed in Ref. [21] in Section VI to investigate whether the classical degeneracy of soliton configurations is broken by quantum corrections. Our results provide the exact one-loop VPEs, extending previous calculations that used the fluctuation spectrum to confirm stability of these solutions [22][23][24][25][26]. A short summary is given in Sec. VII. Technical issues are discussed in two appendixes. In Appendix A we verify the analytical structure of the Jost determinant by numerically simulating contour integrals and in Appendix B we investigate the VPE for the particular case of two un-coupled particles with a mass gap.

II. MULTI-THRESHOLD SCATTERING
We consider a field theory in one space dimension with two scalar fields Φ 1 and Φ 2 . We write the Lagrangian in terms of the field potential U (Φ 1 , Φ 2 ) and assume that the classical field equations produce a localized static solution Φ sol (x) = (Φ sol,1 (x), Φ sol,2 (x)). Fluctuating fields with frequency ω are then introduced as Φ( ). For simplicity, the frequency dependence of the fluctuations φ i (x) is not made explicit. To harmonic order these fluctuations are subject to the wave equation where the diagonal momentum matrix arises from the dispersion ω 2 = k 2 1 + m 2 1 = k 2 2 + m 2 2 and V (x) is the 2 × 2 potential matrix that couples the two fluctuating fields via the soliton background At this point we assume V (x) = V (−x), but we will also discuss a different scenario in Section IV. Note that for investigating the analytic properties of scattering data it is not necessary that the potential arises from a soliton model. As mentioned in the introduction, we take k = k 1 to be the independent momentum variable. Then the dependent momentum variable k 2 = k 2 (k) must have (at least) three important properties: 1. For any real k within the gap, k 2 < m 2 2 −m 2 1 , the dependent momentum k 2 is imaginary with a positive imaginary part, so that ik 2 ≤ 0 leads to a localized wave-function in the closed channel.
2. For any real k outside the gap, k 2 ≥ m 2 2 − m 2 1 , the dependent momentum k 2 is also real and k → −k must imply k 2 → −k 2 . This will ensure F (−k) = F * (k).
3. If Im(k) ≥ 0 then we must also have Im(k 2 ) ≥ 0 such that both momenta are in their respective upper half complex planes.
Properties 1. and 2. seem contradictory because the first one does not allow for a sign change of k 2 while the second one requires it. We will now show that indeed possesses these properties. Property 2. is obvious because outside the gap the real part of the radical is positive and the iǫ prescription that moves the pole into the lower half plane can be ignored. Within the gap, k 2 < m 2 2 − m 2 1 , we expand in ǫ with careful attention to the sign of its coefficient, yielding satisfying Property 1. above. Property 3. is established by introducing k = s + it with real s and t in the relation between k and k 2 , generalizing Eq. (6) We define X and Y such that The square root halves the phase of any complex number, i.e. it maps phases in [−π, π] to [−π/2, π/2], implying X ≥ 0. Also the sign of the imaginary part does not change when taking the square root, so that Y ≥ 0 as well. From this it follows that Im(k 2 ) = |s|sign(t)Y + tX ≥ 0 when t ≥ 0, which is the third required property.
With these momentum variables, we are now prepared to parameterize the Jost solution as a matrix-valued function. Its columns contain the two functions φ 1 and φ 2 subject to specific boundary conditions: Outside the gap, the rows refer to the two independent boundary conditions of out-going plane waves of the individual particles at positive spatial infinity. Inside the gap, the second boundary condition parameterizes a localized wave-function in the closed channel. To this end we apply the differential operator from Eq. (3) to and obtain the second order differential equation for the coefficient matrix function where primes denote derivatives with respect to the spatial coordinate while denote space-independent diagonal matrices. The boundary conditions then translate into lim x→∞ Z(k, x) = 1 and lim x→∞ Z ′ (k, x) = 0.
Since Eq. (3) is real, the physical scattering solutions are linear combinations of F (k, x) and F * (k, x) = F (−k, x). With the assumption that V (−x) = V (x) these solutions decouple into symmetric and anti-symmetric channels where S S and S A are the scattering matrices that are obtained from the boundary conditions Φ ′ S (k, 0) = 0 and Φ A (k, 0) = 0 in the respective channels. For the VPE we merely require the sum of the eigenphase shifts given by the (logarithm of the) determinant of the scattering matrices. Hence we write 2 with the Jost matrices In the symmetric channel, alternative definitions of the Jost matrix F S (k) with different factors of D(k) lead to the same determinant of S S , because D * (k) = D(−k). This ambiguity is removed by demanding that F S (k) be the unit matrix for V ≡ 0, which is equivalent to requiring lim k→∞ F S (k) = 1. The factor of D(k) as in Eq. (14) adds multiples of π to to the phase shift when k 2 ≤ m 2 2 − m 2 1 . In Appendix B we show that these multiples of π produce the correct VPE for the specific case that the two particles are uncoupled, i.e. when the potential matrix V is diagonal.

III. NUMERICAL EXPERIMENTS FOR SYMMETRIC POTENTIAL MATRIX
In this Section use a minimal example to show that the determinants of the matrices defined in Eq. (14), det [F S (k)] and det [F A (k)], have the following properties: • for real k, their real and imaginary parts are even and odd functions, respectively. This implies that their phases are odd (more precisely, the momentum derivatives of their phases are even), • they have no singularities or branch cuts for Im(k) ≥ 0, • for a particular form of Levinson's theorem the phases at k = 0 can be used to count the number of bound states, and • they have zeros along the positive imaginary axis k = iκ i corresponding to bound state energies ω 2 i = m 2 1 − κ 2 i . (For strong enough attraction, ω 2 i may be negative, indicating instability. For the models to be investigated in Section VI, this is not the case.) For decoupled channels, V 12 = 0, bound states of the heavier particle yield zeros at imaginary k 2 . If these bound states have energy eigenvalues m 2 1 < ω 2 i < m 2 2 , these zeros are on the real k-axis. When V 12 = 0, the heavier particle can decay into the lighter one and these bound states become Feshbach resonances [27].
While the first two items serve as tests for the properties established in the previous Section, the latter two provide important information about the zeros of the Jost determinant. They are crucial when it comes to compute integrals involving the scattering phase shift, i.e. the logarithm of the Jost determinant, by means of Cauchy integral rules.
We will confirm these properties by numerical experiments using a specific, but generic potential matrix, which we can easily tune to have an arbitrary number of bound states by changing the width w and the entries of the constant symmetric matrix V 0 . In the numerical simulations we vary the parameter ǫ from Eq.
The second equality gives the standard definition of the scattering phase shift [19] and is valid only if detF * (k) = detF (−k). The proper branch of the logarithm is determined by eliminating 2π jumps between neighboring points and imposing lim k→∞ δ(k) = 0.  The real and imaginary parts of the Jost determinants are clearly even and odd functions of the single momentum variable k, respectively. Discontinuities or even singularities may occur at threshold k 2 = m 2 2 − m 2 1 and at k = 0, which combine with jumps in the phases when bound states exist as required by Levinson's theorem. For repulsive potentials, this theorem implies [28] lim k→0 + δ S (k) = − π 2 . Hence the only phase without a jump is for a repulsive potential in the anti-symmetric channel. Yet, that phase still exhibits threshold cusps. For the attractive potential we also observe rapid changes slightly below threshold. They reflect Feshbach resonances, also called bound states in the continuum [19]. In total, the Jost determinants defined in Eq. (14) combined with Eq. (5) have exactly the expected features describing physical scattering for real k. We have checked these results for other values of V 0 and w as well.
A particularly interesting case is that of a system with interactions only in the heavy particle channel. The phase shift is treated as a function of the lighter particle's momentum and therefore must vanish inside the gap, or at least be a (piecewise constant) multiple of π. If it did not, the VPE would not be the simple sum of the two particles' VPEs when the potential matrix is diagonal. This property is verified in Figure 3. This result relies crucially on including the kinematic factors exactly as in Eq. (14). Also, the jumps by π within the gap are significant for the VPE, as will be discussed in Appendix B.
For the parameters of Figure 3, a close inspection of the symmetric channel reveals that the phase shift jumps from 0 to π at k = 1.270 and from π to π 2 at threshold. These jumps hamper the continuous construction of phase shifts. In the anti-symmetric channel a jump occurs at k = 1.615, while this phase shift is continuous at threshold. Even though the imaginary parts of the Jost determinants are zero within the gap, the jumps at k = 1.270 and k = 1.615 result from sign changes of the real parts as can also be observed from Figure 3. When the off-diagonal elements of V 0 are switched on, Im detF = 0 and detF has no zeros along the real axis. At the same time, the bound states from the uncoupled situation turn into (sharp) Feshbach resonances corresponding to zeros in the lower half of the complex momentum plane. We will comment further on this point at the end of Appendix A.
This discussion also reveals another interesting property of the phase shift related to Levinson's theorem. For single-channel scattering, this theorem states that the phase shift at zero momentum is π times the number of bound states in the anti-symmetric channel, However, for the symmetric channel, the ratio between the phase shift at zero momentum and π equals the number of bound states minus 1 2 [28]. For a two-channel problem one might therefore expect the phase shift at zero momentum to be π times the number of bound states minus 1. The right panels of Figures 1, 2 and 3 show that this is not the case. Rather, the discontinuities that result from Levinson's theorem appear at the respective thresholds. This is evident for the uncoupled system, and persists when the off-diagonal elements of the potential matrix are switched on, though the jump is slightly smoothed out in that case.
To identify genuine bound states as well as Feshbach resonances independently from the Jost function calculation we diagonalize the Hamiltonian associated with the field Equation (3) in a discretized basis. The 2N basis wave-functions in the symmetric channel are with x ∈ [0, L], which is sufficient since the potential is symmetric. The discretized momenta p n are determined such that Φ (0) (L) = 0. Replacing cos(px) → sin(qx) gives the basis wave-functions in the anti-symmetric channel, where we must also redefine the discretized momenta p n → q n in order to maintain Φ (0) (L) = 0. In general the eigenvalues ω 2 of H depend on L. However, this is not the case for bound states when L is sufficiently large, because their wave-functions vanish at spatial infinity. Similarly the Feshbach resonances are associated with those eigenvalues in the continuum that do not vary when L is altered. For the parameters of Figure 2 we observe bound state eigenvalues at ω 2 = −2.128 and ω 2 = 0.603 in the symmetric channel. They correspond to κ = 1.768 and κ = 0.630. This agrees with Levinson's theorem, since the phase shift at zero momentum equals 3π 2 . In the anti-symmetric channel only a single eigenvalue at ω 2 = −0.513 with κ = 1.230 occurs. Accordingly, the phase shift equals π at zero momentum. In the symmetric channel, a Feshbach resonance appears with the eigenvalue ω 2 = 3.435, i.e. k = 1.560, causing the corresponding phase shift to be π 2 at threshold. The Feshbach resonance in the anti-symmetric channel is just below threshold. It has the eigenvalue ω 2 = 3.968, i.e. k = 1.723 and the phase shift jumps by π at threshold. Increasing the potential strength slowly moves this eigenvalue deeper into the gap, e.g. v 22 = 1.2 yields the eigenvalue ω 2 = 3.922, while the resonance disappears for v 22 = 0.6 and the phase shift becomes continuous at threshold, though the typical cusps remains. We conclude that the Jost determinant as defined in Eq. (14) together with the definition of the dependent momentum variable, Eq. (5), possess all required properties for real momenta.
Finally, it is interesting to compute the ω 2 eigenvalues for the parameters of Figure 3, i.e. the case when the lighter particle does not interact. We have observed bound states below m 2 2 at ω 2 = 2.613 and ω 2 = 3.608 in the symmetric and anti-symmetric channels, respectively. These bound state energies correspond to momenta k 2 on the positive imaginary axis that relate to real k when ω 2 > m 2 1 . Hence at these values of k the phase shifts jump by π, as shown in the right panel of Figure 3.
When the coupling to the lighter particle is switched on again, these bound states turn into Feshbach resonances because the heavier particle may decay into the lighter one. Then a key additional property of the Jost determinant is that its only zeros (simple zeros in single channel scattering) are at purely imaginary bound state momenta k = k 1 = iκ i . Therefore we consider the Jost determinant for real t defined via k = it, so that the differential equation, Eq. (3), is purely real (the iǫ prescription can be omitted for t > 0) and so is the Jost determinant. In Figure 4 we show the corresponding numerical results. The positions of the zeros of detF clearly agree with the above listed values for κ i . In the symmetric channel the determinant diverges as t → 0 because of D −1 in Eq. (14). This is actually the origin for the π 2 in Levinson's theorem. This divergence is not a serious problem here because the VPE is only sensitive to t ≥ m 1 . At large t the determinant approaches unity, though this asymptotic value is assumed only slowly. Without any bound states the Jost determinant is positive along the positive imaginary axis.
In appendix A we present further numerical analysis of the Jost determinant using Cauchy integrals.

IV. SKEWED PARITY
So far we have studied the case of a completely symmetric potential matrix. For the model of Ref. [21], we will consider potential matrices symmetric under skewed parity, That is, diagonal and off-diagonal elements of the potential matrix are even and odd functions of the coordinate, respectively. Similar to the Dirac theory we introduce a parity operator that has eigenvalues λ p = ±1. It is then convenient to define the projectors so that φ 1 (x) = P + Φ(x) and φ 2 (x) = P − Φ(x). In this notation the channel with λ p = 1 has while the λ p = −1 channel has The Jost matrices for the two parity channels can also be obtained from the solution of Eq. (10). Suitable projection yields The factor matrices D + (k) = ik 0 0 1 and D − (k) = 1 0 0 ik 2 , with k 2 being again related to k via Eq. (5), ensure that lim |k|→∞ F ± (k) = 1. Note that any right multiplication of F (k, x) in Eq. (9) by a constant matrix is a solution to the field Equation (10) and translates into a right multiplication of Z(k, x) when this constant matrix is diagonal. 3 In the positive parity channel the bound state energies, and eventually the position of the Feshbach resonances, are obtained by diagonalizing the Hamiltonian, Eq. (16) with respect to the 2N basis states Again the discretized momenta p n and q n are obtained from Φ (0) (L) = 0. The basis states in the negative parity channel are constructed by swapping upper and lower components. Using the trial potential matrix we have repeated the numerical experiments from Section III and appendix A and established that the determinants defined by the matrices of Eq. (21) have the standard Jost function properties discussed in Section II.

V. FORMULATION OF VACUUM POLARIZATION ENERGY
Formally the VPE is the sum of the shifts of zero-point energies due to the interaction with the potential V (x). It can be decomposed as [29] E vac = 1 2 where the discrete sum is over the bound state energies (excluding Feshbach resonances) and k labels the continuum scattering states with energy ω(k) such that k = 0 denotes the scattering threshold. According to the Krein formula [30] the derivative of the total phase shift, in channel s gives the change in the density of scattering states. These channels can either be the symmetric and anti-symmetric channels discussed in Sections II and III, or the parity modes from the previous Section. We note that Eq. (23) can also be derived from the energy-momentum tensor in field theory [8]. Finally the subscript "ren." indicates that, as it stands, the integral diverges and requires regularization and renormalization. Working in one space dimension, the Born approximation is typically subtracted under the integral to impose the no-tadpole renormalization condition, which removes all terms linear in V (x) from the VPE. While the pole at k = 0 is canceled by the factor ω(k) − m 1 → k 2 2m1 as k → 0, the singularity at threshold k = m 2 2 − m 2 1 precludes a direct application of this prescription when different mass parameters are involved. Since has the same large k behavior as δ (1) (k), it serves well as a helper function such that is finite. In contrast to δ (1) (k) which has a square root branch cut, ∆ (1) (k) is analytic up the pole at k = 0 which is removed by the factor ω(k) − m 1 . Although E vac is finite, it does not obey the no-tadpole renormalization condition. By analytic continuation we will next obtain an integral expression for the VPE that allows us to replace ∆ (1) by δ (1) and also avoids integrating over the various cusps (or even jumps) of the phase shift within the gap. In a first step we use Eq. (24), replace δ s (k) by the odd function ln detF s (k) and express the VPE as a single integral along the whole real axis, As numerically verified in Appendix A, d dk ln detF s (k) has first order poles with unit residue exactly at those imaginary momenta that represent the bound state energies ω i . When closing the contour in the upper half plane, these poles compensate the explicit contributions from the bound states, i.e. the first term in the above equation. The only contribution then comes from enclosing the branch cut along the imaginary axis (t ∈ R)

Finally we integrate by parts and obtain
denotes the summed exponentials of the Jost determinants. Without difficulty we can now restore the original notadpole scheme, where we have replaced ∆ (1) (it) with It is obvious that this subtraction fully removes the O(V ) contribution to the VPE. Up to a factor of i it is the continuation of the Born approximation, Eq. (25), to the imaginary axis. Formally Eqs. (29)-(30) are standard [4], except for the second denominator in the Born approximation. It clearly reflects the infra-red singularity at k 2 = m 2 2 − m 2 1 above. Fortunately, that is not an issue when we integrate along the imaginary axis. Without off-diagonal interactions ν and ν (1) are the sums of the Jost functions for the two particles with the imaginary momentum variables t and t 2 + m 2 2 − m 2 1 , respectively. Thus it is obvious that Eq. (29) is just the sum of the individual contributions, as it should be. This is also the case for the real momentum formulation, Eq. (27). However, the additivity is less obvious on the real axis because it requires careful consideration of the phase shift within the gap, as we will discuss in Appendix B.

VI. BAZEIA MODEL
The Bazeia model [21] generalizes the φ 4 kink model in one space and one time dimensions by introducing a second scalar field χ. Since the model admits various soliton solutions that are classically degenerate it is a prime candidate to  apply the above developed method to compute the VPE and investigate whether this degeneracy is lifted by quantum corrections. In this Section, we adopt the notation for the fields from the original paper [21], The mass terms are included in the field potential, U (φ, χ), in contrast to Eq. (2). After appropriate rescaling of the fields and the coordinates the Lagrangian is characterized by the single parameter µ. The vacuum configuration is φ vac = ±1 and χ vac = 0 so that m χ = µ and m φ = 2. Subtracting the mass terms from ∂ 2 U yields the fluctuation potential matrix Later we will be interested solely in the case µ ≤ 2. Then χ is the lighter particle and hence plays the role of Φ 1 . There is always the pure kink soliton φ (1) = tanh(x) and In addition, for µ ≤ 1 a second solution is φ (2) = tanh(µx) and The construction of these solitons is described in Ref. [21]. For µ = 1 the two solutions are identical. They are actually degenerate in classical mass, E cl = 4 3 since the model is defined such that E cl can be uniquely determined from the asymptotic values of the fields. The specific cases of the fluctuation potential Eq. (32) for the two soliton solutions (33) and (34) are Since φ is the heavier particle, the translational zero mode wave-function is odd in the upper component and even in the lower one. Hence this zero mode lives in the negative parity channel. However, we also observe a zero mode in the positive parity channel. This second zero mode can be used to construct a family of soliton solutions parameterized by |χ(0)| ≤ 2 µ . This family has already been found in the context of supersymmetric domain wall models [31]. Except for the particular cases µ = 1 2 and µ = 2, which lead to the solitons discussed in Eqs. (37) and (38) below, these solitons can only be formulated numerically; we will provide a detailed discussion of their VPEs elsewhere [32].
For µ = 1 the particles decouple and the heavier one is subject to the potential v 22 (x) = −6/cosh 2 (x) which is known to have a zero mode and a bound state, the so-called shape mode, at ω 2 = 3 (for m 2 = 2) [1]. This corresponds to the imaginary dependent momentum k 2 = i and to the real independent momentum k = √ 2 since m 1 = 1 and m 2 = 2. This is a soliton example of a Feshbach resonance. Hence this state is not explicitly counted in the discrete sum in Eq. (23). Indeed for this real momentum the Jost determinant of the positive parity channel is zero and the phase shift jumps by π.
In Table I Table I is obtained by adding the σ = 2 and σ = µ entries of Table 4 in Ref. [11]. To our knowledge, there are no previous studies on the off-diagonal case, Eq. (36), in the literature.
Note that the µ = 2 entry corresponds to 2 × 1 2 √ 3 − 3 π ∼ −1.3326, i.e. twice the φ 4 -kink vacuum polarization energy for quantum fluctuations with mass m 2 = 2. The case µ = 1 is even more interesting not only because the two solitons are identical, but also because it is an application of our general method to the particular case of zero off-diagonal potential matrix elements, since χ ≡ 0, but with an actual gap present (see also Appendix B). While v 11 (x) = −2/cosh 2 (x) is the potential for fluctuations off the sine-Gordon soliton with mass m 1 = 1, v 22 (x) = −6/cosh 2 (x) is the potential for fluctuations off the φ 4 kink with mass m 2 = 2. Indeed our numerical result equals the sum of these well-established VPEs [1,4], 9846, providing verification of the threshold treatment in particular, since the shape mode of the φ 4 part has disappeared as a manifest bound state. The coupled channel formalism treats the two uncoupled particles in such a way that the shape mode of the heavier particle is no longer a zero of the Jost determinant on the imaginary k-axis. As discussed earlier, this is a kinematical feature induced by relating the momenta between the heavier and lighter particles.
As mentioned, the two solitons are classically degenerate. Obviously including one loop quantum corrections favors the second solution. The present conventions yield E cl = 4 3 suggesting that for small enough µ the soliton, Eq. (34) would unavoidably be destabilized by quantum corrections. These convenient conventions included scaling the fields by 1 √ λ , where λ is the fourth order coupling constant when the coefficient of the quadratic order does not contain λ. However, in quantum field theory the scale of the field cannot be chosen freely; rather, it is dictated by the equal-time commutation relations. Though this is irrelevant when comparing various VPEs in a given model, it must be taken into consideration when combining classical and quantum contributions. When (re)introducing physical parameters, E cl scales like 1 λ , but E vac does not depend on λ [1,4]. It is then clear that the occurrence of instability depends on the interaction strength.
Also note that for µ = 2 the VPE is twice that of the φ 4 model kink even though it is the same soliton configuration. This is due to the addition of χ type quantum fluctuations. The case µ = 2 is also interesting because there is an additional soliton solution and with an arbitrary parameter b > 1. For b = 1 this becomes the standard kink configuration. The classical energy is the same for all allowed values of b. As already observed in the heat kernel calculation of Ref. [12], we find that also the VPE does not vary with b. Actually the Jost determinant itself turns out not to depend on b. The reason is simple: for µ = 2 the transformation φ = η 1 + η 2 2 and χ = η 1 − η 2 2 decouples the potential, Eq. (31) into two φ 4 models for η 1 and η 2 [12]. The configuration of Eq. (37) corresponds to individual kinks for η i separated by arcosh(b). The solution η 1 = tanh(x) and η 2 (x) = 1 translates into φ(x) = [tanh(x) + 1] /2 and χ(x) = [tanh(x) − 1] /2. This case is interesting because it induces the symmetric potential matrix which has a zero eigenvalue. The vacuum polarization energy must then be computed from F S and F A defined in Eq. (14). It is reassuring that our numerical simulation yields E vac = −0.6625, consistent with the known VPE of a single kink. But then, the case µ = 2 is not a stringent litmus test for our approach because it does not exhibit a mass gap. As mentioned, for µ = 1 2 yet another soliton solution is known analytically [31] and with |a| < 1. The case a = 0 reproduces the pure kink soliton. Though the case |a| = 1 solves the field equations, it is not a soliton [21] because it is not localized. Consequently, the corresponding potential matrix, Eq. (32) does not vanish at spatial infinity. Comparing the numerical results Tables I and II, we see that the limiting case a → 1 has an even lower VPE than the second soliton, Eq. (34), for µ = 1 2 .

VII. CONCLUSION
We have generalized the spectral methods for computing the vacuum polarization energy (VPE) of static solitons to models containing coupled fields with different mass parameters. This result is non-trivial endeavor because the spectral methods rely on the analytic properties of scattering data, most prominently the Jost function/determinant F , but different masses induce a non-analytic relation between the momenta that describe the asymptotic behavior of the fields. Taking k to be the momentum associated with the smaller mass, the essential ingredient is to establish a relation, Eq. (5), between the two momenta that fulfills F * (k) = F (−k) for all real k and that preserves the analytic properties of F (k) in the upper half complex momentum plane, Im(k) ≥ 0. We have checked this relation by numerically verifying those properties for different potentials. The formulation of the VPE is then similar to the case without a mass gap. However, it is indispensable to consider momenta off the real axis because the notadpole renormalization condition would otherwise induce a singularity for models in one space dimension. Once the generalized spectral method is established, the computation of the VPE is numerically straightforward and does not require any further approximation or expansion. We have then applied this formalism to a soliton model with two real scalar fields to (i) verify that the approach is consistent and (ii) to show that the VPE, i.e. the quantum correction to the soliton mass, lifts the classical degeneracy.
There are many other soliton models in one space dimension to which the formalism can be applied. For example, the solitons constructed in Refs. [33,34] have different classical masses and it would be interesting to see whether the respective VPEs obey the same inequality.
Obviously the formalism is not constrained to one space dimension but can be applied in higher dimensions with cylindrical or spherical symmetry. A prime candidate for future investigation is the VPE of the 't Hooft-Polyakov monopole [35,36]. multiples of 2πi. The contour of the middle figure encloses all bound states identified in Section III and thus we find 4πi and 2πi for the symmetric and anti-symmetric channels, respectively. The higher bound state in the symmetric channel is not enclosed by the contour of the right panel. Hence in that case we obtain 2πi for both channels. The zeros in the lower half complex momentum plane that correspond to Feshbach resonances cannot be identified by computing I log because the Jost determinant is not expected to be an analytic function in that regime. However, we can test whether the zeros of Figure 3 unexpectedly moved to the upper half plane. To this end we consider a contour defined by the triangle 4 k = 1.0 + 0.0001 i → 2.0 + 0.0001 i → 1.5 + i → 1.0 + 0.0001 i for different values of the off-diagonal elements of V 0 , v 12 = 0.01 and v 12 = 0.2. In the vicinity of the Feshbach resonances and the threshold the numerical integration requires tiny step sizes. In the end the numerical results for I log are indeed compatible with zero.