Large-$N$ $\mathbb{CP}^{N-1}$ sigma model on a finite interval: general Dirichlet boundary conditions

This is the third of the series of articles on the large-$N$ two-dimensional $\mathbb{CP}^{N-1}$ sigma model, defined on a finite space interval $L$ with Dirichlet boundary conditions. Here the cases of the general Dirichlet boundary conditions are studied, where the relative $\mathbb{CP}^{N-1}$ orientations at the two boundaries are generic, and numerical solutions are presented. Distinctive features of the $\mathbb{CP}^{N-1}$ sigma model, as compared e.g., to an $O(N)$ model, which were not entirely evident in the basic properties studied in the first two articles in the large $N$ limit, manifest themselves here. It is found that the total energy is minimized when the fields are aligned in the same direction at the two boundaries.


Introduction
The two dimensional CP N −1 sigma model has enjoyed a special and constant attention of theoretical physicists since the pioneering work by D'Adda et. al. [1] and by Witten [2]. The model is interesting as an analogous model for nonperturbative dynamics of QCD, possessing asymptotic freedom and confinement; it can also be related to some phenomena in condensed matter physics such as quantum Hall effects [3,4,5,6,7].
A context in which this model emerges as an effective action is the quantum excitations of the monopole-vortex composite solitons [10,11,12,13]. Such systems occur under hierarchically broken gauge symmetry, e.g., SU (N + 1) → SU (N ) × U (1) → 1 in a color-flavor locked SU (N ) symmetric vacuum. The 2D CP N −1 model emerges as an effective theory describing the quantum fluctuations of the collective orientational modes of the nonAbelian vortex string [14,15,16]. The vortex boundaries are the massive magnetic monopoles which are generated in the higher-mass-scale gauge-symmetry breaking, and carrying the same orientational CP N −1 moduli. In other words, magnetic monopoles are "confined" by the nonAbelian vortex string. The nonAbelian monopoles arising this way are not plagued by the well-known difficulties (the topological obstruction and non-normalizable gauge zeromodes 1 ), albeit in a confinement phase.
The presence of the boundaries implies that one is dealing with a CP N −1 model on a finite-width worldsheet, with either Dirichlet, Neumann or mixed boundary conditions, depending on the details, such as the mass hierarchy ratios. These were part of the motivations for our previous work [8,9].
More generally this type of models are interesting on their own as a prototype of a quantum system of mixed dimensions. Since it possesses a dimensionless parameter LΛ consisting of a length of the string L and the mass gap Λ for an infinite length string, it interpolates between the known CP N −1 model in 2D in the LΛ → ∞ limit and a quantum mechanical (classical in the case of Dirichlet model) system in the LΛ → 0 limit.
In Refs. [8,9], the generalized gap equations have been studied analytically and solved numerically, for a wide range of values of LΛ. The energy density has then been studied carefully, by subtracting a quadratic divergence (an analogue of the QCD vacuum energy) consistently. The results show how the dynamical mass generation (well known in the 2D CP N −1 sigma model) and the classical L → 0 limit are consistently described by our solutions. In particular, it was found that the system possessed a unique phase for any L, which is smoothly connected to the "confining" phase of the 2D CP N −1 model in the L → ∞ limit. Also, the approach to the L → ∞ limit has been studied and shown to be purely exponential, with no Lüscher-like power-behaved terms present. Finally, the Casimir force has been studied and the presence of different regimes (repulsive force at moderate values of L ∼ O(Λ −1 ) and attractive force at large L, corresponding to a constant string tension) was predicted.
In the present work, we turn our attention to the cases in which the relative orientation, characterized by an angle α, of the CP N −1 fields at the two Dirichlet boundaries is generic (α = 0 in general). The richer structure of the gap equations is illustrated and the possible dependence of the solution on α is discussed. These equations are then solved numerically for various values of α and for different lengths L of the string, and the dependence of the total energy on the relative orientation is investigated, both analytically and numerically. We find that the total energy is minimized when the two orientations are parallel. Distinctive features of the CP N −1 sigma model manifest themselves here, in contrast to the first two articles, where in the large N limit the distinction from other sigma models such as the O(N ) model, was not always evident, as far as the static properties of the system were concerned.
The paper is organized as follows. In Sec. 2 we generalize the case of generic Dirichlet boundary conditions and show that the parameter space is governed by an angle α ∈ [0, π/2]. In Sec. 3 we present the new numerical method that we have employed and the numerical solutions found. Sec. 4 addresses the α-dependence of the total energy and Sec. 5 concludes the paper with a discussion. Some details of our calculations have been delegated to the Appendices A-C.

General Dirichlet boundary conditions
The generalized gap equation discussed in Refs. [8,9] has the form and r 0 ≡ 4π/g 2 stands for the bare coupling constant. The CP N −1 fields are split into σ(x) ≡ n 1 (x) and n i (x) for i > 1 and the latter are integrated out giving rise to the modes {f n }, and positive-definite energies (eigenvalues) {ω 2 n |ω n ∈ R >0 }. The shift of N 2π in Eq. (2.2) arises from an anomalous λ variation associated with the divergences in the sum over modes [9].
The equations (2.1) correspond to the Dirichlet boundary condition D-D : namely the orientation in the CP N −1 space was taken to be the same at the two boundaries. The fields are thus defined on the finite-width worldstrip x ∈ [− L 2 , L 2 ]; the length of the string is L. The other parameter in the model is the internal scale Λ which also sets the energy units. The only physical parameter is thus LΛ.
In the present paper, the cases in which the fields are orientated differently in the CP N −1 space at the two boundaries will be studied. Due to the global SU (N ) symmetry of the CP N −1 model and by the definition of the CP N −1 coordinates the most general D-D type boundary condition can be taken in the form, Accordingly, the generalized gap equation becomes together with the boundary conditions (2.5)-(2.7) and we have defined the two "classical" fields: σ 1 ≡ n 1 and σ 2 ≡ n 2 . {f n (x), ω 2 n } are the solutions of the Schrödinger equation Actually it is useful to start with the observation that the differential equation for any given λ(x) has two linearly independent solutions. A particular convenient choice turns out to be the solution σ R (x) and σ L (x) with the following characteristics. Near the left boundary, x = − L 2 , the two solutions behave as [8] σ The normalization of the divergent behavior (the first term in Eqs. (2.12) and (2.13)) is fixed by the gap equation, W is a constant determined later, and L 0 is a certain Length parameter. 2 For even λ(x), given a solution σ(x), the parity transformed function σ(−x) is also a solution. It follows that Exactly at the boundaries, the gap equation tells us that and similarly The Wronskian of the two solutions are defined by By evaluating the Wronskian (which is constant) near the two boundaries, the constant in Eqs. (2.12) and (2.13) is seen to be precisely the Wronskian itself. Let us define also the even and odd solutions: 2 Strictly speaking, we need to introduce another UV cutoff parameter b to impose condition (2.5) and condition (2.6) at a distance b from one of the boundaries. Using Eq. (2.2), we find the two UV cutoff parameters are related to each other as (2.14) That is, if we fix a ratio b / in the limit of → 0, L 0 is a given parameter and thus, independent of L and α.
In terms of σ R (x) and σ L (x), the most general solution can be written as which corresponds to the pair of boundary conditions In other words, the vectors q R a and q L a , describe the orientations of the classical field at the right and left boundaries, respectively. By using the global SU (N ) invariance of CP N −1 , one may choose 23) and this corresponds to the general boundary condition, anticipated in Eqs. (2.5)-(2.7).
The general boundary condition may be put in a form which looks more symmetric with respect to the exchange of the two boundaries. This can be done by rotating the CP N −1 frame by the angle α/2. By writing only the first two components of Eqs. (2.5)-(2.7), where .
(2.25) σ 1 andσ 2 are the components of the classical field in the new CP N −1 frame. By using the orientation vectors, this means that

The gauge field
In our first two papers [8,9], the A 0 = 0 gauge was chosen, which is appropriate for studying static vacuum configurations. Also the gauge field A x was assumed to be absent in the vacuum (the functional stationary point). Note that the equation of motion for A x is where n and σ are column vectors and hence n † n is the inner product. A x then satisfies and this vanishes if the solution for σ is real. More precisely, as equation (2.11) is real (λ(x) is real), the field σ can always be chosen to be real by the local U (1) gauge transformation. Eq. (2.29) shows that A x vanishes in the vacuum in such a gauge; a fact used in Refs. [8,9].
Similarly, the configuration with the most general boundary condition (2.20) gives the following (classical part of the) U (1) current After making the U (1) gauge transformation so that A x = 0, the right hand side of the above must vanish. Our choice with real q R a , q L a , in Eq. (2.23), thus yields J U (1) x = 0, corresponding to the choice of gauge A x = 0.

Range of α
Exchanging the boundaries x = L 2 and x = − L 2 and using an appropriate rotation, a configuration with α is seen to be equivalent to the one with −α.
Also, two solutions with α and with π − α can be regarded as the same boundary conditions. For consider the left boundary where we have used a global SU (2) ⊂ SU (N ) e.g., in the 2 − 3 plane first (which does not affect q R a ), and a U (1) equivalence relation q L a ∼ e iβ q L a at the end. Note that in the CP N −1 sigma model the U (1) symmetry (2.4) is implemented as a local gauge symmetry. Therefore the solution with the boundary condition can be regarded as a gauge transformation of the solution (2.23). In order to have the same physics (e.g. the same energy density, etc.), however, one must appropriately introduce the gauge field A x . Note that even though the σ 1 field is chosen to take real values at the two boundaries, it is necessary that it goes through a phase rotation, meaning that such a solution necessarily generates a current J x (x), hence a nonvanishing gauge field A x (see the previous Subsection). Once they are appropriately taken into account, the π − α solution (2.33) is simply a gauge transformation of the α solution, (2.23).
Accordingly, the range of the parameter α can be taken as without loss of generality. The result of Section 4 indeed indicates that the solutions with α in this range are the stable ones.
Nevertheless, the solution with α and the one with π − α are distinct if A x ≡ 0 and if σ 1 and σ 2 are kept real. In section 3, we will discuss the numerical solutions of the gap equation (2.8) and (2.9) for the whole range of α, i.e. [0, π]. Numerically, this is advantageous because we can avoid introducing the gauge field A x and we can keep the σ fields real in the calculations.

2.3
The "solutions" with α = 0 and α = π A particular case of interest about the gauge (non-)equivalence of the solutions with α and π − α, concerns the solutions with α = 0 and α = π. The classical field σ a (x) as a function of x ∈ [− L 2 , L 2 ] can be regarded as a path from a point in CP N −1 to another, through the "interior" of it: in the space CP N −1 × R. The solution α = 0 corresponds to the movement from the left to the right boundary, as this is the case studied in Refs. [8,9]. It describes a closed path in CP N −1 × R. Other solutions 0 < α < π involve two components σ 1 and σ 2 nontrivially, and in general do not describe a closed path.
It is interesting to consider the solution with another particular boundary condition α = π: it looks like (by using the first line of Eq. (2.24)) It may appear that the two solutions (2.36) and (2.37) correspond to topologically inequivalent classes of paths and that both solutions might hence be stable. Actually, (1, 0, . . .) and (−1, 0, . . .) represent the same point of CP N −1 . 3 In other words the α = π solution also describes a closed path with the same starting and end points as in the α = 0 solution. Since the space CP N −1 × R >0 is simply connected, these two closed loops are homotopically equivalent. Therefore only one of the solutions can be stable. The study of the α dependence of the energy in Section 4 indicates that the solution with α = 0 is indeed the stable one.

The solutions of the gap equation
With the most general D-D boundary condition, the equations to be solved are together with the boundary conditions (2.21) and (2.23). In Ref. [9], we used a randomwalk method for solving the first of the above equations. In this paper, however, we use a relaxation method akin to the procedure of solving the time evolution in a heat-like equation. In order to derive the equation that we will solve numerically, let us start from the energy functional of the form [9] from which a variation with respect to λ yields indeed this is how the generalized gap equation was derived in Ref. [9]. Although the above equation is formally convergent due to the regulator, it is not practical for numerical calculations as it still includes an infinite number of modes. Therefore, we will switch regularization to a finite number of modes, n ≤ n max . For consistency, we need to modify the coupling r accordingly was defined in Ref. [8].
Now, instead of setting Eq. (3.5) equal to zero (which it should be), we will introduce a fictitious time dependence and flow said fictitious time evolution where in the last equality, we have introduced a discretized first-order time derivative, h τ is the time step, and t is an index of the current time slice. Rearranging, we can finally write the formal evolution equation where h t should be chosen appropriately; in particular, it should be small enough to ensure convergence of the algorithm.
We are now almost ready to perform the numerical calculations. First we discretize the string into a one-dimensional lattice. We need to calculate the Schrödinger modes and energies, {f n , ω 2 n } from Eq. (2.10), and we will do that simply by diagonalizing the discretized version of the operator (−∂ 2 x + λ) subject to the boundary conditions f n (± L 2 ) = 0. Next, we will solve the equations of motion for σ 1,2 , (3.2), by back-solving the discretized version of the operator (−∂ 2 x + λ) on their respective boundary conditions (see Ref. [9] for details); in particular, we will use the conditions of Eqs. (2.21) and (2.23). We start the algorithm by using a guess for λ. The final step is to calculate the new λ from Eq. (3.8). This is the end of the cycle; now we start over by calculating the modes and energies and so on. The cycle continues until the integral of the absolute value of Eq. (3.5) is smaller than an appropriate small number, which we shall call ε numerical , see Appendix B.
The numerical solutions to σ 1 (x) and σ 2 (x) for α = nπ 16 with n = 0, 1, 2, . . . , 16 are given in Fig. 2. In order to distinguish the solutions for different alpha, we introduced a color scheme used throughout the paper, which is defined as 2α being mapped to the color circle 4 . In order to distinguish α = 0 and α = π, we changed the color for the latter to black. For a legend with all the colors, see  .23). An asymmetric appearance of σ 1 and σ 2 is due to the choice of the CP N −1 coordinates (the first line of (2.24)). The combination, σ 2 1 + σ 2 2 , which appears in the gap equation (3.1) is shown in Fig. 3. As σ 2 1 + σ 2 2 is invariant under the rotations in the 1 − 2 plane in CP N −1 , the picture is indeed symmetric under the exchange of the two boundaries. As can be seen clearly, it shows a very little dependence on α, except in the middle region of the string. In particular the value σ 1 = σ 2 = 0 is reached only for α = π and only at the center of the string, x = 0. This corresponds to σ 2 ≡ 0 on the entire string and σ 1 crossing through zero at the midpoint. 4 2α is mapped to the color angle called the hue. In our convention: α = 0 is red; α = π 3 is green; α = 2π 3 is blue and the anti-colors are in between: α = π 6 is yellow; α = π 2 is cyan and finally, α = 5π 6 is magenta. The gap function λ(x) is plotted in Fig. 4. Consistently with the previous Fig. 3, the gap function depends on α significantly only in the central region of the string. In order to better see the symmetric nature of the two boundaries, one can use the form of the boundary condition (2.26), (2.27) (see Fig. 1). The parametric plot of the solutions in terms ofσ 1,2 is shown in Fig. 5.
These results contain several interesting features which are not always manifest. To reveal some of them, we introduce the CP N −1 variables rather than the homogeneous coordinates n i . As only two components n 1,2 (σ 1,2 ) are involved in the solution, it is sufficient to use the CP 1 = S 2 variables related to them by where τ i are the Pauli matrices. Also, as in the solutions considered n i are real, one may restrict oneself to s 3 = σ 2 1 − σ 2 2 , s 1 = 2σ 1 σ 2 , s 2 ≡ 0 . Finally, we define some alternative CP N −1 variables by using the more the "symmetric" boundary conditions (2.26) and (2.27); thus by using theσ i 's given in Eq. (2.25). Using these variables, the L dependence of the solutions are illustrated in Fig. 7, for fixed values of α, α = π 4 and α = 3π 4 , as an example of two solutions with α and π − α, see Subsection 2.2.
It is seen from the figure that only for small values of L, the parametric solution bends inwards (for α ≤ π 2 ) or outwards (for α ≥ π 2 ). As L tends to infinity (for fixed and finite Λ), the parametric solution tends to exactly the diagonals. That is, in these coordinates, the solution will come in with an angle π 2 − α with respect to the y axis and go out with an angle α − π 2 . Indeed, as the two points α and π − α are the same point in CP N −1 , only one of the two solutions can be stable and the other will be metastable. The figure alludes to the claim that the solution with α ≤ π 2 is the stable one; in this case, the solution with α = π 4 , whereas the other is a metastable solution.
In order to make evidence for our claim, we will study the α dependence of the total string energy in the next section.
4 Dependence of the energy on the relative orientation α Let us now analyze how the total energy:  Figure 7: A parametric plot of the "symmetric" CP N −1 variables for α = π 4 (above the two diagonal straight lines emanating from the origin) and α = 3π 4 (below the same two diagonal lines), for Λ = 1 and L = 1, 2, 3, . . . , 12. The color scheme is shown in the legend and is unrelated to that defined for the angles α.
where σ a is taken to be real and repeated indices are summed over: a = 1, 2. Note that as λ(x) and σ a satisfy the equation of motion, the contribution from the interior of the string, (− L 2 , L 2 ), vanishes, and only the surface term remains. We will now use the general solution, the first line of (2.24), to write By using parity one can write this as (4.6) one gets the estimates The evaluation of terms involving ∂σ L ∂α requires a more careful consideration: although the leading behavior of σ L in (4.6) is independent of α, the α dependent effects coming from subleading terms may not be negligible. A study in Appendix A, however, shows that As for the first term of Eq. (4.5), it can be evaluated straightforwardly: where again the definition of the Wronskian (2.18) has been used. As one gets the net result The energy increases monotonically with α.
Even though at large L the Wronskian W is believed to be exponentially small 6 , O(e −LΛ ), so the α dependence of E is very small, this is not so at smaller L ∼ O(1/Λ). One expects a significant dependence of the energy on α. See Fig. 8  found W for various values of α, as a function of L. One sees that the exponential L dependence at large L is universal, i.e., independent of α, which is quite understandable, as the effect of mis-alignment at the CP N −1 variables at the far boundaries should be unimportant at large L.
Eq. (4.12) is one of the main results of the present work: it proves that, in the U (1) gauge in which the σ 1,2 (x) fields are real throughout [− L 2 , L 2 ] and A x ≡ 0, the solutions in the range 0 ≤ α ≤ π 2 (Eq. (2.35)) are the stable ones, as claimed.

Numerical checks of ∂E ∂α
The total energy of the system can be expressed in various ways: where we have integrated f 2 n by parts, used the fact that the boundary term vanishes [9], the equation of motion for the modes f n and the completeness relation (2.10). The latter expression was the starting point in Section 3 and defining it can be written neatly as Starting again from the first expression (4.13); instead of integrating by parts, we can collect the terms that are multiplied by the gap function λ: can be shown to be a constant (independent of x), by using the gap equation [9]. Also the equation of motion of λ N 2 n f 2 n ω n e − ωn + σ 2 1 + σ 2 2 = r , r = r 0 + N 2π , (4.19) has been used in the last step.
As the equations of motion for σ a and λ, Eq. (2.8) and Eq. (2.9) have been used in moving among the lines, these expressions are equivalent on-shell, i.e., when evaluated at the minimum of the action. For the purpose of numerically verifying the α dependence of the energy E found analytically above, however, use of different expressions will provide us with nontrivial, independent checks.
A particularly subtle issue in our discussion concerns the divergences. The constant part of the energy density, E 0 , is quadratically divergent, even after the logarithmic divergence in the sum over modes is eliminated by the standard coupling constant renormalization; the subtraction prescription [9] is implicit in Eqs. (4.13)-(4.17). Moreover, due to the behavior of the gap function near the boundaries, integration of the gap function in Eq. (4.17) diverges linearly, which should be canceled by the bare "monopole mass" terms [9]. Both the quadratic divergence of the energy density and the linear divergence of the λ integration are local effects, the former around a generic point in the string, x, and the latter at the boundaries; therefore, they should be independent of the relative CP N −1 orientation α at the two boundaries. The derivation of Eq. (4.12) relies on this tacit assumption: it is a highly nontrivial check whether this result is reproduced by the evaluation of some of Eqs. For the numerical check, a possibility is to use the formula which follows from differentiating Eq. (4.14) with respect to α. The dependence on α through the function λ(x), proportional to the gap equation (2.8), has been dropped (i.e., made use of), hence Eq. (4.22) expresses ∂E/∂α through the functional dependence on σ a only. If the equations of motion for σ a (Eq. (2.9)) were also used, the only thing that remains would be the surface terms -which have been evaluated analytically in Section 4.
Another possibility is to use Eq. (4.17) instead: where the constancy of E 0 can be exploited to evaluate the first term ∂E 0 ∂α , e.g., at the midpoint of the string, x = 0, where the numerical precision is best.
The result of such a comparison is shown in Fig. 9, for L = 1 and L = 4 (and Λ = 1). Shown in the Figure are the numerical evaluations of Eqs. (4.22) and (4.23) and the analytical expression ∂E/∂α = +2 sin α W , where the Wronskian is calculated numerically (Fig. 8). Apart from small numerical errors, the overall agreement between the direct numerical evaluation and the analytic formula is quite satisfactory. They clearly confirm our conclusion that the solutions in the range α ∈ [0, π 2 ] are the stable ones.

Discussion
In this work we have further examined the quantum vacuum configuration 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. Building on the results of our preceding work [8,9], the systems with two generic CP N −1 orientations at the two boundaries are studied here, for various values of the string length L and for different relative orientation angle α. The total Eq. energy of the system E at fixed L is found to increase with α monotonically. Taking into account of the defining properties of the CP N −1 sigma model, this means that the relative angle between the classical field orientations at the boundaries can be limited to [0, π 2 ], without loss of generality. The system has the lowest energy when the CP N −1 orientation is the same at the two boundaries.
The classical field component, in fact, traces a path from a point to another in CP N −1 × R >0 , in going from the left boundary to the right boundary. It is on a point in CP N −1 at the left boundary; it goes into its "interior", before emerging at another point of the CP N −1 surface at the other boundary (see Fig. 1). Due to the fact that the space CP N −1 × R >0 is simply connected, the solution with α and the one with π − α correspond to two paths which are homotopic to each other, as discussed in Subsection 2.2. Thus only one of them can be stable. The fact that the energy E monotonically increases with α in α ∈ [0, π] shows that the stable solutions correspond to those with α ∈ [0, π 2 ]. In this discussion the specific property of the CP N −1 model, not shared by other sigma models such as the O(N ) model, turned out to be crucial.
The instability of the solutions with α ∈ [ π 2 , π] cannot be seen perturbatively. We have indeed verified that our "solutions" with α > π 2 do not suffer from any zero or negative modes {f n , ω n } even if the potential λ(x) becomes negative in the central region of the string (the left of Fig. 14), see Fig. 10. The instability is a nonperturbative phenomenon. It would be an interesting problem to understand better the nature of such instabilities.

Note added
After this work is finished, we were informed by Muneto Nitta of their new paper [24] which deals with a similar system. As far as we can see our approach and the results obtained are different from theirs.
Let us consider the following quantity with a set of coefficients {b R,L n } given by There are relations between coefficients as b R n = (−1) n b L n due to the parity f n (−x) = (−1) n f n (x). It can be seen that ∆ α σ R,L (x) satisfies and it must vanishes at both boundaries. In other words, ∆ α σ R,L (x) is a normalizable, zero-energy solution of a Schrödinger equation. With a positive definite potential λ(x), the only solution is It follows that ∂σ R,L /∂α can be expanded in {f n (x)} as, Here the summation with respect to the infinite modes is assumed to converge even around the boundary since the α dependence does not affect UV modes. Therefore we find that their behavior for x ∼ −L/2 is, with constants B R,L obtained by, where {a L n } and {a R n } are defined as Using these properties, we find that

B Numerical accuracy
In this section, we will try to quantify the numerical accuracy of the solutions presented in this paper. The handle we have from the relaxation method is given by the quantity which is the numerically evaluated gap equation. In the numerical calculations presented in this paper, we have used ε numerical = 10 −4 . If we here denote by x the numerical accuracy, it is easy to see that The numerical lattice we used has LEN = 10 4 lattice points; thus locally on average, we have Notice, however, the combination σ 2 1 + σ 2 2 used in the gap equation has the local error, on average, which is very small.
Using now the fictitious time-flow equation (3.8), we can estimate the error of λ as where we have defined ∆λ ≡ λ t+1 − λ t and used that h τ < 1. Thus, we have also for λ, that the numerical precision on average, locally, is Using error propagation, we can estimate the precision of the Wronskian as Indeed, in Fig. 8 it is not possible to see any numerical error at the 10 −6 level.
Let us now estimate the precision of the α derivative of the energy using Eq. (4.22); starting with the first parenthesis, we have It is a bit harder to estimate the precision of the alternative formula for the α derivative of the energy, Eq. (4.23). Let us consider the error estimate of each term of the energy density in turn. The first term (in the brackets) in Eq. (4.18) can be estimated as the error of n ω n due to the normalization of the eigenmodes, which we can write as √ n max ε numerical LEN (B.10) Ignoring the spatial derivatives (as they are quite precise, see above), all the other terms in Eq. (4.18) as well as λ itself have an error level of ε numerical . Since, for n max = 350 and LEN = 10 4 , the largest error is ε numerical , we get the net precision of the α derivative using Eq. (4.23) is the same as using Eq. (4.22). In practice, however, the numerical error of using Eq. (4.23) is a bit larger than of Eq. (4.22); the difference is an order-one factor, which we did not evaluate.