Large-N CP(N-1) sigma model on a finite interval and the renormalized string energy

We continue the analysis started in a recent paper of the large-N two-dimensional CP(N-1) sigma model, defined on a finite space interval L with Dirichlet (or Neumann) boundary conditions. Here we focus our attention on the problem of the renormalized energy density $\mathcal{E}(x,\Lambda,L)$ which is found to be a sum of two terms, a constant term coming from the sum over modes, and a term proportional to the mass gap. The approach to $\mathcal{E}(x,\Lambda,L)\to\tfrac{N}{4\pi}\Lambda^2$ at large $L\Lambda$ is shown, both analytically and numerically, to be exponential: no power corrections are present and in particular no L\"uscher term appears. This is consistent with the earlier result which states that the system has a unique massive phase, which interpolates smoothly between the classical weakly-coupled limit for $L\Lambda\to 0$ and the"confined"phase of the standard CP(N-1) model in two dimensions for $L\Lambda\to\infty$.


Introduction
Recently we embarked on the investigation of the bosonic CP N −1 model [1,2], defined on finite space interval L, i.e., on a finite width worldstrip, in the large N approximation [3]. Such a system could provide a useful model for various physical situations. For instance, it appears as the low-energy effective theory describing the quantum excitations of the monopole-vortex soliton complex [4,5,6,7] in hierarchically broken gauge symmetries, such as SU (N + 1) → SU (N ) × U (1) → 1 in a color-flavor locked SU (N ) symmetric vacuum. The CP N −1 model describes the nonAbelian orientational zeromodes of the nonAbelian vortex (string) [8,9,10], whereas its boundaries represent the monopoles arising from a higher-scale gauge-symmetry breaking, carrying the same orientational CP N −1 moduli. NonAbelian monopoles, not plagued by the well-known difficulties, could emerge in such a context. The fate of the nonAbelian monopoles as a quantum mechanical entity is then linked to the phase of the low-energy CP N −1 effective action attached to it.
In [3] it was found that the quantum saddle-point equations describing the CP N −1 model with Dirichlet or Neumann boundary conditions has a unique solution under certain conditions. In the large-L limit, this solution approaches smoothly the well-known confining phase of the standard 2D CP N −1 system. A phase transition between a Higgs-like phase and the confining phase for a shorter L, which was claimed to be present in the literature [11], was shown not to exist in the system. 1 The model is interesting also from a formal point of view, as it provides a prototype model of a quantum system of varying dimensions in the presence of dynamical mass generation: it interpolates between a 2D QFT (in the L → ∞ limit) with all well-known phenomena such as asymptotic freedom and confinement and a 1D system in the L → 0 limit -quantum mechanics. For shorter strings of length L ≤ 1/Λ, quantum fluctuations of the CP N −1 fields n i , (i = 1, 2, . . . , N ) remain weakly coupled, as they lack sufficient 2D spacetime "volume" in which the fields fluctuate. With the Dirichlet condition, the system reduces effectively to a classical system in the L 1/Λ limit.
In this paper, we delve in more detail into the properties of the large-N CP N −1 model on a finite-width worldsheet. First, with a more refined numerical method we improve the precision of the solution to the generalized gap equation. This enables us to explore a larger region of the parameter space and in particular the limit of large L. The second problem is to understand the energy density of the string itself as a function of x, computed at the functional saddle point, completing the analysis presented in [3]. The third problem is to clarify the approach to the QFT (L → ∞) limit of our system; this involves the question of certain consistency with the known field-theory limit, as well as of figuring out interesting L-dependent effects. It will be seen that power-behaved corrections such as the Lüscher term are absent. This is consistent as all fields acquire dynamically generated mass; at the same time no spontaneous breakdown of the global SU (N ) symmetry takes place.
The paper is organized as follows. In Section 2 we review the CP N −1 model on a finite strip and also present new numerically improved results which allows to reach higher values of L than before. In Section 3 the generalized gap equation is re-derived, paying special attention to the anomalous term that arises in the functional variation, which is analogous to the axial anomaly. In Section 4 we consider the energy density in detail, its various contributions, and its L → ∞ limit. In Section 5 we study an analytical Ansatz that describes the large but finite L 1/Λ and the approach to the L → ∞ limit. The numerical results for the renormalized energy density are presented in Section 6. In Section 7 we discuss the Casimir force. Our conclusion is in Section 8. Some details of our analysis are given in Appendices A ∼ F.

Review of the CP N −1 model on a finite width worldsheet
The classical action for the CP N −1 sigma model is defined by where n i with i = 1, . . . , N are N complex scalar fields and the covariant derivative is given by D µ = ∂ µ − igA µ . Configurations related by U (1) gauge transformations n i → e iα n i are not only gauge-equivalent, but are equivalent because the U (1) gauge field A µ does not have a kinetic term in the classical action. λ is a Lagrange multiplier field that enforces the classical gap equation where r = 4π g 2 is related to the gauge coupling and can be thought of as the "size" of the CP N −1 manifold.
For the CP N −1 theory on a finite interval of length L, x ∈ [− L 2 , L 2 ], 2 the boundary conditions must be specified. One possibility is the Dirichlet-Dirichlet boundary condition which -up to a U (N ) transformation -is D-D : For the moment we take the boundary conditions for the n i fields in the same direction in the CP N −1 space at the two boundaries. Another possibility is the Neumann-Neumann boundary condition 3 N-N : In this paper, we will focus on the Dirichlet-Dirichlet boundary condition. Thus with this condition the N fields can naturally be separated into a classical component σ ≡ n 1 and the rest, n i (i = 2, . . . , N ). Integrating out the n i fields yields the effective action: Because the effective action only depends on |σ| 2 and |∂ µ σ| 2 and the boundary conditions take real positive values on both sides, one can take σ to be a real field and set the gauge field to zero. Finally, we will consider the leading contribution at large N only.
The generalized gap equations following from Eq. (2.5) (see also Section 3 below) have been solved numerically in [3], by a Hartree-like self-consistent method. The renormalized, finite functions of x, Λ, L, λ(x) and σ(x) have been obtained numerically for various values of L and Λ. The calculations of [3] have been extended to larger values of L, with a considerably improved method. The weak point of the Hartree-like method -from a numerical point of view -is the need to determine λ(x) from the second equation in Eq. (2.6), where σ tends to zero in the middle of the string.
It is sometimes convenient to rewrite the first equation in Eq. (2.6) as in terms of the two-point function The latter satisfies, for the D-D boundary condition, an equation (2.10) Note that the infinite number of mirror poles are required to satisfy the D-D boundary condition. See Appendix A.2 for more details.
Under the assumption the near-the-boundary behavior of the fields turns out to be [3] (see Section 2.2 below):

Numerical method and solutions
The new method is based on a random-walk algorithm and is reversed in some sense with respect to the old method. A guess can be made for the function λ(x), but the precise starting point is not important. The algorithm has two assumptions built in; basically just for saving computational costs; i.e. λ(x) is a symmetric function in x; the second is that λ(x) is a monotonically increasing function from the midpoint of the string to the boundary (both assumptions are indeed consistent with the results of [3]). Now the algorithm makes a random change to a part of the function λ(x) (viz. on an interval that is a randomly chosen subset of the full string interval) yieldingλ(x). Now the newλ(x) function is tested in the following way. σ(x) is calculated from its equation of motion (second equation of Eq. (2.6)) with the appropriate Dirichlet boundary conditions as well as from the generalized gap equation (the first equation of Eq. (2.6)); these two are compared. If the newλ(x) makes the two σs move closer to each other, then the newλ(x) function is accepted as the new improved λ(x), otherwise it is rejected. Then the cycle repeats until the precision is good enough (for the solutions we found, that is Some examples are shown in Figure 1 where the approach to the confined phase at large L is evident. The values of fields in the middle of the interval σ 2 (0) and λ 2 (0) − Λ 2 are shown in a logarithmic plot in Figure 2. The approach of σ 2 (0) to zero is clearly exponential and consistent with its mass. In Figure 3 we show the solutions in the interval (− L 2 , 0) by keeping one boundary fixed at − L 2 . This clearly shows the convergence at large L to the half-line solution (−∞, 0]. As already noted in [3], one can clearly see that the asymptotic (L → ∞) regime (with  λ(x) ∼ Λ 2 and σ(x) ∼ 0 except near the boundaries) has already set in at L 4, which is quite reasonable (LΛ = 4 1). The effect of the boundaries is seen to propagate only for ∆x ∼ 1/Λ from the latter: the system effectively reduces to the standard 2D CP N −1 model in an infinite spacetime, as one moves away from the boundaries by ∼ 1/Λ or more, as expected on general grounds. 2.2 L dependence, near-the-boundary behavior of λ(x), σ(x), and the classical limit Our system has two parameters, the dynamically generated mass scale Λ and the interval length L. But as Λ fixes the physical unit of length, the model actually possesses only one parameter: what distinguishes two physically distinct systems is the product ΛL. The crucial point is that the UV divergences are short-distance effects around any fixed space point, and are universal. They do not depend on the presence or absence of the boundaries. This is what allows us to define unambiguously systems of "different space width L".
The near-the-boundary behavior (2.12) can be obtained [3] as follows. Under the assumption (2.11), the large n modes are given at any fixed x by The finite part of the sum over modes in the gap equation behaves then as (2.14) near the left boundary. This singularity can only be compensated by σ(x) 2 in the gap equation, hence Eq. (2.12). The numerical solutions found in [3] and here clearly exhibit this logarithmic behavior.
There is an alternative way of understanding the behavior of σ(x) and λ(x) near the boundaries. Consider the regularized but un-renormalized form of the gap equation (2.6). By keeping the UV regularization parameter fixed and by going to x = ∓ L 2 , one finds that as f n = 0. This is nothing but the classical CP N −1 model (2.2) with the Dirichlet boundary condition With x close to but not exactly at a boundary, in (2.15) is replaced by x ± L 2 and one finds Eq. (2.12). This statement requires an explanation. What really happens is that in the gap equation (2.6), log 1/ which is in σ(x) at exactly x = −L/2, is moved at small but nonvanishing |x + L/2| to the first term involving the sum over the modes. After log 1/ is eliminated by the bare coupling constant term r and the gap equation is renormalized and made finite, it produces − log 1/|x + L/2|, which can only be compensated by σ(x) 2 , as in (2.12 This discussion clearly shows that the origin of the singular behavior of the mass gap λ(x) and of the σ(x) field is the fact that the system reduces to its classical limit 4 near the boundaries, not having sufficient 2D spacetime volume for the n i fields to fluctuate.
The same reasoning explains 5 the behavior of the value of λ(x) and σ(x) at the midpoint of the string at small L 1/Λ, found in [3] (see Fig. 3 there), (2.17)

Anomalous functional variation and the generalized gap equation
We now re-derive the generalized gap equation. Our starting point is the energy density where has been introduced as a regulator of the UV divergences coming from higher modes, E uv is a subtraction constant and f n (x), ω n are the eigenmodes of the n i field equations By integrating Eq. (3.1) over x ∈ − L 2 , L 2 and by using (3.2), one has the expression for the integrated energy, For instance the total derivative terms vanish as [3] lim By varying (3.3) with respect to λ(x), and by using one gets the generalized gap equation whereas the variation with respect to σ gives Note the extra N 2π term in (3.7). It arises when the variation δ/δλ(x) acts on the regulator factor e − ωn : this term is superficially of the order of O( ): however the sum in the last expression diverges as 1 , so it gives a nonvanishing contribution 6 . As the divergence comes from large n it may be calculated as where use was made of an approximate form for the eigenmodes The same result can be found by using the propagator representation (2.9). See Appendix A.2.

Energy density
We now go back to the density itself, and rewrite (3.1) as by collecting terms proportional to λ(x) and by using Eqs. (3.6) and (3.7). Note that the anomaly (3.7) is crucial to give the term proportional to λ(x) in the energy density, after using the gap equation.
It turns out that E 0 (x) is a constant. The space derivative of E 0 (x) is: Noting that the expression in the square bracket above is the derivative of the gap equation (3.6), we conclude that We thus find that the energy density is where the only dependence on x is through the function λ(x).
Let us study the L → ∞ limit of this expression. The value of the constant part of the energy density, E 0 , in the L → ∞ limit may be calculated by noting that λ(x) → Λ 2 and the fact that the spectrum is exactly known in that limit. See Appendix B. The result is This shows that the energy density of the system, after the standard regularization and renormalization of the coupling constant has been made to render the gap equation finite, still contains a quadratic divergence. This is a little similar to the vacuum density in QCD: the theory can be renormalized and all physical quantities can be calculated order by order in perturbation theory, but the vacuum energy density (a contribution to the cosmological constant) is still divergent, and requires a further subtraction. The result (4.6) however suggests that we take the vacuum energy subtraction constant simply as and the constant part of the energy density is thus one finds that the total energy density approaches a constant at any fixed finite x. This gives the quantum corrections 8 due to the fluctuations of the n i fields to the classical "tension" of the vortex, ξ , where Λ 2 ξ . This result is in agreement with the one [12], found in a finite worldstrip CP N −1 model with periodic boundary conditions.
As we shall see in the next section, and as can be verified by a WKB analysis done in Appendix D, the divergence in the energy density remains purely quadratic: the only subtraction needed is (4.7), in the case of finite (L) string also. No linear or logarithmic divergences are present. This is reasonable as the divergences due to the fluctuations of the n i fields is a short-distance effect, local in x, and cannot depend on the presence of the boundaries.

Large but finite L: an Ansatz and analytic calculation
We now study the corrections to (4.9) and (4.10) (and σ(0) = 0) for large but finite L. In order to do that, it is clearly necessary to analyze the large-L behavior of (the solution of) the gap equation, (3.6), (3.8) itself. To start with, λ(x, Λ, L) can be approximately taken to be a constant for |x| L/2; we parametrize its asymptotic approach to Λ 2 as The factor a represents the nonlocal effect due to the boundary condition. To estimate this, consider the propagator, Eq. (2.9). For x, x ∼ 0 (i.e., far from the boundaries), it satisfies locally thus its solution can be assumed to have the form with an unknown constant A = A(ΛL). 9 The subleading terms in the second line come from the nearest mirror poles in Eq. (2.10). 10 It is expected that A = A(ΛL) ∼ O(1), but due to the effects of the boundaries where λ(x) is non constant and singular [3] it will not coincide with the exact value A = 1 (for the Dirichlet boundary condition at x = ±L/2) or A = −1 (for the Neumann boundary condition) (cfr. Eq. (A.1)) [3].
The gap equation (3.6) then yields near x = 0 By requiring the consistency of this with the initial Ansatz (5.1) and making use of the identity one finds that Thus one finds at large L 1 Λ that, around x = 0, With these results in hand, one can now calculate the asymptotic behavior of the energy density itself: It turns out under the same approximation (5.3) that the constant part of the energy density, E 0 (Λ, L), is given at finite large L by The derivation of this result is given in Appendix C. Note that the divergence in the energy density is just the purely quadratic one, N π 2 , the same as in the L → ∞ case discussed in the previous section. This is correct, as the divergences arise from the UV fluctuations of the n i fields which is a local effect, independent of the boundaries, or of the value of L.
Finally one finds, by adding the λ(0, Λ, L) term and by making the same subtraction as before, viz. (4.7), where another identity has been used.
To conclude, we find that the approach to the asymptotic value of the energy density is exponential: no pure power corrections in 1/L (i.e. the Lüscher term) are present. This is perfectly consistent with the general result found in [3] that our system has a unique phase, which smoothly matches -in the large L limit -the "confinement phase" of the standard 2D CP N −1 model. All n i (i = 1) fields gain a dynamically generated mass ∼ Λ; at the same time σ ∼ 0 except at the boundaries. In other words, no dynamical breaking of the isometry group SU (N ) takes pace. No Nambu-Goldstone modes associated with the internal, orientational modes are generated. The absence of a long-range correlation in the large-L corrections in Eqs. (5.9), (5.12) is a simple reflection of this fact.

Numerical results
Due to the quadratic divergence present in the sum (4.2) the numerical calculation turns out to present quite a bit of a challenge. Any tiny errors in the eigenmodes and in the energy levels will introduce linear or logarithmic divergences in the sum, and the finite answer for E(x, Λ, L) one gets (including its dependence on x, L and Λ) depends on how these fake divergences are appropriately subtracted, together with the genuine quadratic divergences. Because of this, even the best results so far do not have a precision comparable to the solution of the generalized gap equation discussed in Section 3.
The check of the constancy of E 0 , is shown in Fig. 4. Note that it is found indeed to be constant everywhere, including values of x very close to the boundaries, where the distances from the latter are much smaller than 1/Λ. Fig. 5 shows the value of E 0 as a function of L (Λ = 1). The total energy density, including the N λ(x)/2π term, calculated at the midpoint, x = 0, is shown in Fig. 6. The numerical results are nicely consistent with the exact result at L = ∞, and with the analytic behavior for large but finite L, found in the previous section, with A ∼ 1. E 0 and E(0, Λ, L) are plotted against Λ at fixed L, in Fig. 7 and in Fig. 8. In particular we see that as Λ → 0 the energy density converges, although quite slowly (i.e. logarithmically), to the free-field value −π/12L 2 .

Boundary divergence and the Casimir force
The energy density after renormalization is a finite function of x, Λ, L. When integrated over x ∈ − L 2 , L 2 , it gives the total energy of the string where we introduced the mass gap function defined on the interval [0, L], The Casimir force is defined as 11 Before analyzing the behavior of F for various values of L, let us note that the second term in Eq. (7.2) gives rise to a new divergence in the integrated energy due to the singular behavior near the boundaries, e.g., near x = 0, Asλ(x) quickly approaches the constant value Λ 2 beyond x 1/Λ, for sufficiently large L, the divergent part can be extracted by considering the finite integral where a UV cutoff (x = ) has been introduced. Similarly E 2 for the contribution from the right boundary. E 1 (E 2 ) is an energy concentrated at the left (right) boundary: it can be interpreted as the quantum corrections to the monopole (antimonopole) mass, due to the n i field fluctuations (the factor N ). E 1,2 can be subtracted (i.e., compensated with the bare mass terms) from the total energy, leaving finite, renormalized monopole masses. They do not affect the discussion on the L dependent Casimir effect below.
The Casimir force can be rewritten, by differentiating Eq. (7.2), as Note that this is a finite quantity since it is not sensitive to the leading divergence ofλ, viz. (7.5).
At very small L (i.e. LΛ 1) one expects the dominant effect to come from the second term of (7.7) (see Fig. 8): it is an attractive free-field Casimir force.
At intermediate values, L ∼ O(1/Λ), instead, we find that the force is dominated by the third term in (7.7). The strong decrease (∼ 1/L 2 log(1/L)) of the mass gapλ(x) with L for all x (see e.g., Eq. (2.17)) cannot be compensated by a linear effect of integration, therefore there is an effective repulsive force at work. The total energy of the system is lowered when the space interval L gets larger.
At sufficiently large L, L 1/Λ, where the 2D regime sets in, the leading contribution comes from the first term in (7.7): corresponding to an approximately constant string tension (Eq. (4.10)). An external observer who attempts to pull the boundaries further apart will experience an attractive, constant force countering her/him.
A precise numerical verification of this nontrivial behavior of the force turned out to be exceedingly difficult because of the singular behavior of λ(x) near the boundaries. Our preliminary result (not shown) however clearly confirms the change from a repulsive regime at L = O(1/Λ) to an attractive force at L 1/Λ.

Conclusion
In this paper we have examined the energy density function E(x, Λ, L) of the large-N CP N −1 sigma model on a finite string, defined with the Dirichlet boundary conditions. We find that it is a sum of two terms, the first expressed as a sum over fluctuation modes, which turns out to be constant in x, and the second term proportional to the mass gap λ(x). The only x dependence arises from the second. The first term is quadratically divergent, analogous to the vacuum energy in QCD. The L-dependence of E(x, Λ, L) at fixed Λ shows that the effect of the boundaries are limited to their vicinity of width ∼ 1/Λ: the system approaches quickly the standard 2D large-N CP N −1 sigma model, with dynamical generation of the mass gap, and with no dynamical breaking of the isometry group SU (N ).
In the small LΛ limit, the system approaches the classical weakly-coupled CP N −1 model, as appropriate for the Dirichlet boundary conditions.
The approach to the limit E(x, Λ, L = ∞) = N 4π Λ 2 is found to be purely exponential: no power corrections in 1/L such as the Lüscher term are present. This is perfectly consistent with the general result found in [3] that our system has a unique phase, which smoothly matches the "confinement phase" in the large-L limit of the standard 2D CP N −1 model. All n i (i = 1) fields gain a dynamically generated mass ∼ Λ; at the same time σ ∼ 0 except at the boundaries. In other words, no dynamical breaking of the isometry group SU (N ) takes pace, and no associated Nambu-Goldstone modes are generated. The absence of the long-range correlation reflects this fact.
Recently a paper appeared [16] in which some analytical large-N CP N −1 sigma model solutions for inhomogeneous condensates are presented by a mapping to the Gross-Neveu model [15]. These solutions correspond to periodic boundary conditions and, as far as we can see, none of the solutions proposed there correspond to our system defined with the Dirichlet boundary conditions. It would certainly be very interesting if our type of solution could be found analytically with developments of these techniques in the future.
Comments on Ref. [17] In Ref. [17], submitted to the ArXiv on the same day as ours, the same system is analyzed in a different approach, and the authors there claim to find analytic solutions for mass gap function λ(x) and for σ(x), both in confinement and in Higgs phases.
By imposing Dirichlet boundary conditions and solving the generalized gap equations, we instead find the solutions (e.g., illustrated in Fig. 1) in a unique phase with mass gap for all values of L, which smoothly approaches the well-known solution in the infinite L limit (the standard 2D CP N −1 model). Our solutions are moreover consistent with the classical CP N −1 model in the L 1/Λ limit as discussed in Sec. 2.2.
It is possible that, if a Higgs-like solution (as in Fig. 1(b) of [17]) would exist, it represents an unstable solution, whereas our procedure necessarily picks up the stable solution, and that actually the confinement-type solution is always the stable one.
However, we find it difficult to make a proper comparison, as the renormalization of the gap equation and the generation of the mass scale Λ are not explained in [17].
The boundary behavior of the mass gap function and the field σ given in [17] is powerlike, whereas the logarithmic behavior found by us reflects the situation characteristic of a finite-space-width system. The system must compromise between the 2D physics at L ≥ 1/Λ -the divergences of the n i field fluctuations and the generation of the mass scale Λ -and the classical limit to which the model must reduce correctly in the L 1/Λ region, as explained in Sec. 2.2.
As the physical values of L (the space width) are not given in reference to 1/Λ in [17], in contrast to what is done in the present paper, it is not clear to us which physical values of L their solutions in Fig. 1(b) or Fig. 1(a) refer to, for instance.
As a consequence, it is unclear how and when (at which value of L) the Higgs phase vacuum disappears, as L is increased. Or, vice versa, at which L, if L decreases toward zero, the Higgs vacuum takes over, if it does at all. The authors of [17] do not give the criteria to decide which solutions should be chosen at any given L. As far as we can see, the analysis of the vacuum energy density, as made in the present paper, has not been done yet there.
A Propagator D(x, τ ; x , τ ) A.1 Exact forms with λ(x) = m 2 When λ(x) = m 2 , Eq. (2.10) can be easily solved by for the Dirichlet-Dirichlet boundary condition (2.3). In particular, note that (A.2) In the case of the Neumann-Neumann boundary condition (2.4), the sign of the last terms of the r.h.s. is flipped.

A.2 Alternative derivation of the anomalous functional variation
In terms of the propagator (2.9), the extra factor in Eq. Since the UV divergence comes from a short-distance effect, to extract divergent terms of D(x, , x, 0), it is sufficient to consider the contribution of the nearest poles, δ(τ −τ )δ(x−x ) (n = 0 term) in Eq. (2.10). Furthermore, at short distances λ(x)(|x − x | 2 + 2 ) 1, the potential λ(x) can be omitted and thus the propagator behaves as one for a massless field in two dimensional space, With a general potential λ(x), therefore, the divergent part of D(x, ; x, 0) (for − L 2 < x < L 2 ) is universal as In the simplest case with λ(x) = m 2 , one can easily check this property using Eq. (A.1) and K 0 (m ) ∼ log 1 . We find which gives the extra constant term − N 2π in the gap equation. Similarly, in a region where the dominant behavior of the propagator is At large L and at finite x, where λ(x) ∼ Λ 2 , one can make an approximation valid at all levels n (simply assume λ = m 2 = Λ 2 ). Then As E 0 has been shown to be a constant, it can be calculated at any fixed x, for example at the midpoint x = 0, where σ = 0: In the L → ∞ limit, the sum may be replaced by an integral, by πn L → z. Also, let us make a replacement e − ωn → e − πn/L , i.e., in the exponential damping factor. One finds In going from (B.1) to (B.5), however, we made a replacement (B.4) in the exponential damping factor. The correction due to this approximation must be taken into account. The effect of this replacement can be studied by writing so it gives a finite contribution. It is (B.14) This can be easily calculated to give This must be added to (B.7) obtained under the approximation (B.4): the final answer is To compute E 0 (Λ, L) at large but finite L we observe that the constant part (evaluated at x = 0) of the energy density (4.2) can be written as (see Eq. (2.9)) By using (5.2) this can be rewritten as where we set σ (0) = 0 by symmetry. We now use (5.3): and (5.1), to get where the identity has been used. Finally at small z, and we made the replacement in the last line. In the L → ∞ limit (C.8) approaches the function E 0 (Λ, L = ∞), calculated in Appendix B, exponentially fast.

D WKB analysis
Assume that for a given value of L, λ(x) has been found. We adopt the WKB approximation to the Schrödinger equation in order to study the nature of the divergences in i.e., the high-n behavior of the summand. As E 0 is constant in x, we shall set x = 0, where σ (0) = 0.
The WKB quantization condition is given by 12 where a ∼ − L 2 , b ∼ L 2 for large n. The wave function and its derivative are given by to first order in (implicit here). Near the boundaries λ(x) behaves as Let us check the large-L limit first. There 12 Due to the sharp rise of the potential λ(x) near the boundaries, the phase shift in the WKB wave function is 0 rather than π 4 (the Maslov index being 0 rather than 1). One has n instead of the familiar n + 1 2 on the right hand side of (D.3). The situation is analogous to the case of the rigid wall. We thank G. Paffuti for discussions on this point. L p(x) = L ω 2 n − Λ 2 = πn ; (D.10) This leads to the calculation for E 0 (Λ, L = ∞) described in Appendix B. To find the corrections, write The quantization condition is corrected to As |∆ n | is small compared to n, it may be calculated by inserting the zeroth-order WKB for ω n (D.11) in (D.13). Thus (D.17) A straightforward calculation leads to The last ingredient needed is the large n behavior of ∆ n . It is easy to estimate at large n. It follows from (D.18), (D.19) that at large n, where C 1 and C −2 are constants of order of unity. No n 0 and n −1 terms appear.
Thus the divergences in E 0 is purely quadratic and is equal to N π 2 , the same as in the L → ∞ system.

E Random walk algorithm
In this appendix, we will describe the algorithm we have used for the numerical calculations in more detail using pseudo-code.
For the numerical calculations, we can only include a finite number of modes in the sum in the left-hand equation in Eq. (2.6), henceforth we shall denote this number as nmax.
Next we have to discretize all the numerical functions on the interval on a lattice with LEN lattice points, which we for convenience will take to be an odd integer. As explained in the text, we will use the symmetry of the problem to make λ manifestly symmetric (with respect to x → −x) in the calculation.
Take a λ which is guessed or just λ = 1 (we started indeed with this) Now we need a function to calculate the error of using the current λ as compared to the true solution. We will define the function where Delta is the discretized second-order differential operator. The '\' notation is an implemented operator in MATLAB and Octave for a linear-algebra operation sometimes called back solving. Formally it is equivalent to multiplying by the matrix inverse of Delta from the left. Numerically, however, that is much more computationally expensive and hence one should instead use a back solving algorithm. In Mathematica it is implemented as a function called LinearSolve. Then we calculate σ again using the gap equation sigma2sq = r [V , D ] = eigensystem ( Delta ) for i = 1: nmax fn = V ( i )/( sqrt ( hx * sum ( V ( i )* V ( i )))) sigma2sq = sigma2sq -fn^2/(2* sqrt ( D ( i ))) end for Most computational packages have a built-in function for finding the eigenvectors and eigenvalues of a given matrix, here we will call it eigensystem and denote by V the eigenvectors and by D the eigenvalues. Other programming languages have libraries for linear algebra manipulations that include such a function, e.g. LAPACK for Fortran90 or CLAPACK for C. Now calculate the error as err = hx * sum ( abs ( sigma1^2 -sigma2sq )) end function Start the algorithm err = 1 errtol = 1e -5 while ( err > errtol ) do Randomly select an interval that should be changed istart = round ( LENHALF * random ()) istop = round ( LENHALF * random ()) Decompose λ into a difference vector diffvec = lambda ( LENHALF -1: end -1) -lambda ( LENHALF : end ) Act on the selected range with a random multiplication factor and a random addition diffvec = [ diffvec (1: istart -1) , scalefactor * random ( istop -istart +1) .*( diffvec ( istart : istop ) + additionfactor * random ( istop -istart +1)) , diffvec ( istop +1: end )] where .* denotes an inner product on the vector space. The addition factor is necessary in the beginning if one chooses to start with λ = 1. At the end of the convergence, it should be small or turned off. 13  If the discrepancy between the σ calculated from the equation of motion and the σ calculated from the gap equation has decreased, then store the new λ and continue; on the other hand, if the error has increased, then discard the new step and try again. The cycle continues until the error is small enough (set by errtol).
Various small tweaks can be implemented in the algorithm depending on the part of parameter space one is interested in. Those tweaks, however, just make the algorithm converge faster, but to the same solutions.
We should mention that if one suspects that the guess will converge to a local vacuum and not to the true vacuum of the functional space, then the metropolis algorithm can be used to accept increases in the error at an initial stage of the random walk. When the error decreases or when the running time increases, this allowance of "going in the wrong direction" should then be decreased. On the grounds of knowing the solutions from Ref. [3], we have not used this possibility in most of the calculations.

F Algorithm test
In this appendix, we will test the algorithm by choosing a poor initial condition, i.e. λ ini (x) = 0, and check which solution the algorithm will find. Most solutions presented in the text were found by starting with a much better guess for λ.
Slightly more advanced than what is described in App. E, we will run the algorithm on a computing cluster and only the best improvement of each cycle will be accepted.
Since the algorithm prefers the largest decrease in the numerical error (err) at all times, the first thing it wants to do is to bring down σ(0) towards zero. This happens very quickly by randomly adding arbitrary values to λ near the border, see Fig. 9. Recall that the algorithm is programmed to make λ a monotonically increasing function on the interval [0, L/2]. The algorithm randomly chooses where and how much to increase the function and uses the gap equation to accept or discard the random steps. Unfortunately, the randomly chosen (by the algorithm) values of λ near the boundary yield too "sharp" a solution for σ; to mitigate this, the algorithm sees the numerical error can be reduced by "pushing out" the corners of λ and adjusting the midpoint, λ(0), which after enough cycles yields a solution for σ to the gap equation and hence a solution for λ, see Fig. 10. The algorithm terminates when the error is below a given acceptable threshold (errtol). The solution is shown as a black line in Fig. 10; i.e. this solution has been accepted with an error tolerance of errtol = 6 × 10 −5 .
Finally, in Fig. 11 we display the midpoint values of λ and σ 2 as functions of acceptance numbers (which roughly corresponds to running time of the numerical calculation).