Casimir squared correction to the standard rotator Hamiltonian for the O($n$) sigma-model in the delta-regime

In a previous paper we found that the isospin susceptibility of the O($n$) sigma-model calculated in the standard rotator approximation differs from the next-to-next to leading order chiral perturbation theory result in terms vanishing like $1/\ell\,,$ for $\ell=L_t/L\to\infty$ and further showed that this deviation could be described by a correction to the rotator spectrum proportional to the square of the quadratic Casimir invariant. Here we confront this expectation with analytic nonperturbative results on the spectrum in 2 dimensions, by Balog and Heged\"us for $n=3,4$ and by Gromov, Kazakov and Vieira for $n=4$. We also consider the case of 3 dimensions.


Introduction
In the pioneering paper [1] Leutwyler showed that to lowest order in chiral perturbation theory (χPT) the low energy dynamics of QCD in the δ−regime is described by a quantum rotator for the spatially constant Goldstone modes. We recall that for a system in a periodic spatial box of sides L the δ−regime is where the time extent L t L and m π L is small (i.e. small or zero quark mass) whereas F π L, (F π the pion decay constant) is large.
There are other important physical systems, in particular in condensed matter physics where anti-ferromagnetic layers are described by the O(3) sigma-model for d = 3 [2], where the order parameter of the spontaneous symmetry breaking is an O(3) vector. In the analogous perturbative regime these systems and also the non-linear sigma models in d = 2 are described by a quantum rotator to leading order.
In all such systems the lowest energy momentum zero states of isospin I have to leading order χPT energies of the form E I ∝ C n;I , (1.1) where C n;I = I(I + n − 2) , (1.2) is the eigenvalue of the quadratic Casimir for isospin I. At 1-loop level it turns out that the Casimir scaling (1.1) still holds, but it is of course expected that at some higher order the standard rotator spectrum will be modified. The standard rotator describes a system where the length of the total magnetization on a timeslice does not change in time. This is obviously not true in the full effective model given by χPT. To our knowledge, the actual deviation from the Casimir scaling was first observed by Balog and Hegedüs [3] in their computation of the spectrum of the d = 2 O(3) nonlinear sigma model in a small periodic box (circle) using the thermodynamic Bethe ansatz (TBA).
Of course, a deviation from the standard rotator spectrum could be established by explicit perturbative computations. However, in our previous paper [5] we pointed out that by comparing the already obtained NNLO results for the isospin susceptibility calculated in χPT at large ≡ L t /L with that computed using the standard rotator one can establish, under reasonable assumptions, that the leading correction to the rotator Hamiltonian occurs at 3-loops and is proportional to the square of the Casimir operator with a proportionality constant determined by the NNLO LEC's of χPT.
In Appendix A we illustrate that such corrections are expected by recalling the spectrum of the O(3) rotator in d = 1 with lattice regularization. Note that the appearance of terms proportional to C 2 I is not a lattice artifact. This simple example serves only to show that distorting the Lagrangian of the standard 1d rotator by an O(3) invariant perturbation leads naturally to a spectrum which (at a higher order) contains such term.
After reviewing some preliminary results in Sect. 2, in Sect. 3 we test our claim above in the d = 2 O(3) and O(4) non-linear sigma models in a periodic box. For O(3) we find excellent agreement of our prediction with the analytic computations of the lowest isospin 1 and 2 energies for O(3) by Balog and Hegedüs [3].
For O(4) Balog and Hegedüs [4] computed only the ground state energy; Gromov, Kazakov and Vieira [6] computed also higher state energies but the results at small volumes presented there are not sufficiently precise for our purposes. We thus generated more data; our methods used to solve the TBA equations are described in Appendix C. Again there is good agreement with our prediction.
The derivation of our result depends on the validity of a (plausible) assumption; but also there is, to our knowledge, no rigorous derivation of the TBA equations (or even the S-matrix) from first principles starting with the 2d O(n) QFT 1 . Hence the agreement of the results provides extra evidence for the validity of both scenarios.
Finally in Sect. 4 compute the effect for the d = 3 O(n) model, but we have not yet found data for which a comparison can be made.

The isospin susceptibility
In this section we recall some results on the isospin susceptibility. The Hamiltonian of the O(n) standard quantum rotator with a chemical potential coupled to the generatorL 12 of rotations in the 12-plane is where Θ is the moment of inertia. In d = 4 dimensions to lowest order χPT one has Θ F 2 L 3 . The isospin susceptibility is defined as the second derivative of the free energy wrt h: In ref. [5] we showed that for small u = L t /(2Θ) the isospin susceptibility computed from the standard rotator, which we call χ rot is given by On the other hand in a previous paper [7], we computed the isospin susceptibility in an asymmetric L t × L d−1 box (with periodic boundary conditions) and the mass gap, in this case the lowest energy in the isospin 1 channel, to NNLO (next-to next-to leading order) χPT. For the susceptibility we recall the results in eqs. (3.54)-(3.57) of [7] with dimensional regularization: and and R (b) 2 involves the 4-derivative couplings which we shall not include in this paper. The expressions above require some explanation. Firstly g 0 is the bare coupling for d = 2; g −2 0 = ρ s , the spin stiffness for d = 3; and g −2 0 = F 2 for d = 4 where F is the pion decay constant in the chiral limit. V D is the volume V D = L t L D−1ˆ q where with DR we have added q = D − d extra dimensions with extent Lˆ ; in this paper we will usually setˆ = 1 unless stated otherwise. (Note that the renormalized quantities do not depend onˆ .) I nm;D are 1-loop dimensionally regularized sums over momenta (cf eq. (3.44) of [7]) 2 . Finally W in (2.7) is an integral associated with a two-loop vacuum massless sunset diagram, which is discussed in Appendix B.
The lowest energy state above the vacuum (the mass gap) is given by E 1 = m 1 in eqs. (5.9)-(5.11) of [7]: Here R(z) is the propagator for an infinitely long strip without the slow modes and W is a 2-loop sunset integral discussed in Appendix B.
For the simple rotator (2.1) at zero isospin chemical potential, the energy gap is given by By inserting the expression for Θ using (2.8) into (2.3), we obtain the susceptibility for small u as a function of F, L, . This can then be compared to the direct χPT computation (2.4) in the -regime for 1. The comparison requires knowledge of the large -behavior of shape functions and the sunset integral appearing in the latter, which are discussed in Appendix B of ref. [5] for d = 4.
The two results in NNLO differ by ∝ 1/ terms for 1, (plus terms vanishing exponentially with ): for d = 4 and similarly for d = 2, 3. In ref. [5] we showed that the deviation above can be accounted for if the spectrum to the order we are considering is given by a modified rotator with eigenvalues of the form E I (L) = 1 2Θ(L) C n;I + Φ(L) L C 2 n;I . (2.13) Θ(L) appearing here contains a higher order correction with respect to the one in (2.11), denoted below by Θ old (L). (2.14) The leading order of the coefficient Φ(L) in (2.13) is directly related to ∆ 2 in (2.12). In [5] we considered only the case d = 4 ; in the next sections we consider the lower dimensions d = 2, 3.
It is interesting to note that the correction discussed here shows up at NNLO in the isospin susceptibility χ, and in the NNNLO (next order) in the spectrum of the effective rotator. The reason is that in the (standard) rotator approximation the Boltzmann weight is ∼ exp(−g 2 C I L t /L) = exp(−g 2 C I ) where g 2 = g 2 (1/L). The typical values of the isospin in the partition function are then given by C I ≈ 1/(g 2 ), hence the extra factor exp(−C 2 I Φ(L)L t /L) in the Boltzmann weight gives a correction δχ/χ ∼ −Φ(L)/(g 4 ) for the susceptibility. With Φ(L) = Φ 3 g 8 + O g 10 (cf. [5]), the leading ∝ C 2 I term in the effective Hamiltonian which is of NNNLO, results in an NNLO, ∝ g 4 / term in δχ/χ.

Delta regime in d = 2
The susceptibility computed in χPT is for d = 2 given by (cf eq. (3.67) of [7]) where g MS (q) is the minimal subtraction (MS) scheme running coupling at momentum scale q. The function r 2 ( ) appearing above is given by (cf eq. (3.68) of [7]) where w( ) appears in the expansion of the sunset diagram defined (see (B.5)).
The susceptibility computed from the simple rotator is given by (cf eq. (2.4) of [5]): where g FV is the LWW running coupling [8] defined through the finite volume mass gap: Its expansion in terms of the running coupling in the MS scheme of DR is given by (cf eq. (2.2) in [8]): The first two coefficients were computed in ref. [8] where Z ≡ ln(4π) − γ , (γ = −Γ (1)) . (3.11) The 3-loop coefficient c 3 can be obtained by combining Shin's 3-loop computation [9] of the finite volume mass gap using lattice regularization, with the computation of the 3-loop coefficient of the lattice beta-function by Caracciolo and Pelissetto [10]. 3 The result is given by 14) Next we insert the expansion (3.8) in (3.6) and compare the resulting expression with the perturbative result given in (3.1). First we see that the g 0 This is satisfied since (using (3.4)) which is a special case of the general relation Matching at order Using eq. (B.11) and the fact that κ 10 (∞) is finite consistency requires consistency is obtained if We arrive at the result On the other hand, from our considerations of the modified rotator (cf eq. (11.5) of [5]) we would expect where Φ 3 is the leading coefficient in the perturbative expansion of Φ(L), assuming the expansion starting at order g 8 MS : As a consequence, the low-lying spectrum to order g 8 MS is given by Hence we conclude, for example, where α MS = g 2 MS /(2π). In the next subsection we test this prediction for n = 3 and n = 4 for which independent computations of the spectrum of excited states exist. The first computations for O(3) were done by Balog and Hegedüs [3] using the TBA equations with the knowledge of the exact Zamolodchikov S-matrix [14]. Later Gromov, Kazakov and Vieira [6] computed the excited states for O(4) from Hirota dynamics; in earlier papers Balog and Hegedüs [4,15] computed the mass gap E 1 for O(2r) but not E 2 .

Running coupling functions
For the perturbative analysis of their data, Balog and Hegedüs [4] introduced a function g 2 J (L) of the box size L 5 through where b 0 , b 1 are the universal first perturbative coefficients of the β−function: and Λ FV is the Λ−parameter of the LWW finite volume coupling in (3.7). For Λ FV L small there are usually two solutions of (3.34) for g 2 J (L), one large and one small 6 . We chose the solution which is small for Λ FV L small, which has the property It is appropriate to call g J a "running coupling function" since it satisfies the equation i.e. like a perturbative running coupling with β−function Balog and Hegedüs consider g 2 J as a function of z = M L where M is the infinite volume mass gap: The ratio M/Λ FV is known from the work of Hasenfratz and Niedermayer [17] In particular, (3.42) 5 analogous constructions can be made for functions of lengths in infinite volume 6 The solution of the equation , first studied by Euler in 1783, and the index −1 refers to the real branch for which α ≈ −1/ ln(X) → +0 when X → +0. The LWW coupling has the following expansion in terms of g 2 J : with the coefficient j 3 given by In table 1 we reproduce the data of Balog and Hegedüs [3]; (f 3,est appearing in the last column is defined in (3.48)). Firstly we remark that by fitting the data for the mass gap E 1 one obtains a value of j 3 close to that given in (3.45) deduced from Shin's analysis.  Secondly, from the fifth column of table 1, we see that although the ratio E 2 /E 1 is close to the ratio of the Casimir eigenvalues, the data of Balog and Hegedüs establishes that the simple effective rotator model requires corrections.
We have for the I = 2 mass for n = 3 the expansion With our value of f 3 we get a very small 3-loop coefficient which is similar to the remark in [3] (before their eq. (4.7)) that they numerically established that the 3-loop coefficient was zero. Although Shin's numbers (which have not yet been checked by an independent computation) are consistent with the data, to indicate whether or not j 3 − 8f 3 is non-zero would require the values of the s i to a greater (at least a factor 10) accuracy than that obtained by Shin.  To see even more clearly the agreement of our analysis with the data of Balog and Hegedüs, in fig. 1 we plot estimates for f 3,est given by for the case n = 3 (and for n = 4 , see next subsection). Comparing the quadratic and cubic fits to the data in the range 0.05 < α J < 0.089 (the first 8 points) one gets f 3 = 0.151 (2), which is in excellent agreement with our prediction in (3.30).

The O(4) case
In [4] the 1-particle data for z ≥ 0.001 for the O(4) model is sufficiently precise for our purposes, however, the 2-particle results are missing. Gromov, Kazakov and Vieira [6] computed the 2-particle state energies, but for z < 0.1 we required higher precision than those Afterwards we have rewritten the code, discretizing the continuum formulation and using Fast Fourier Transform (FFT) to calculate the convolutions appearing. Due to FFT the calculation of an iterative step becomes much faster and yields more precise values. In addition, approximating the 1d continuum functions by a finite lattice of N points, and using several N values (typically N = 2 10 . . . 2 14 ), one can get the continuum limit by determining numerically the coefficients of the corresponding 1/N expansion. Some details are given in Appendix C.
With this technique we could lower z, the size of periodic box down to z = 0.0005 . For z < 0.0005 it becomes more time consuming to produce reliable estimates 9 , but the values of α J (z) (see table 3) are already quite small.
Our data for the TBA energies E I for I = 0, 1, 2 is presented in table 2. The corresponding values for the excitation energies for isospins 10 , I = 1, 2 together with α J are given in table 3.
In that table we see that Casimir scaling holds approximately for a large range of volumes. Also we show estimates for f 3,est given by (3.48); which are plotted in fig. 1 together with quadratic and cubic fits. In the range 0.0437 < α J < 0.075 (7 points) 8 We thank Nikolay Gromov for sending these data to us. 9 Note that with decreasing z the region of attraction shrinks, and it becomes more and more difficult to start the iterative method in the corresponding region. 10 In the sector with a given isospin I only states with particle number ν ≥ I contribute. Among these the ν = I is the lightest, and the TBA equations used describe just the ν = I = 1 and the ν = I = 2 states.
with for D = 3: and, (recall W = −L 4−2D Ψ), We require the large behaviors of the shape functions which we take from eqs. (B.20) and (B.27) of [5]: 3 ( ) α For the mass gap we have from eqs. (5.9)-(5.11) of [7] (settingˆ = 1) where we have set the coefficients of the four-derivative couplings l 1 = l 2 = 0. From eq. (3.16) of [18] for D ∼ 3: The expressions for χ and χ rot agree at 1-loop order up to terms vanishing exponentially with since (see eq. (5.10) of [18])  Suppose the spectrum is given by (2.13) and that Φ(L) has the expansion starting at order (ρ s L) −4 : then from our considerations of the modified rotator (cf eq. (11.5) of [5]) Nyfeler and Wiese [19] have investigated the histogram of the uniform magnetization for the spin 1/2 Heisenberg model on the honeycomb lattice in the cylindrical regime (see eqs. (3.9)-(3.12) and fig. 4.3). Having previously determined the low-energy constants ρ s and c for this model in the cubical regime, they took the NLO formula for Θ and the standard rotator spectrum, to compare the cylindrical regime. The agreement seemed to be very good, but their statistical errors are too large to see signs of our predicted Casimir 2 term (or even the NNLO term in Θ). We are not aware of any other measurements of the spectrum of the d = 3 O(n) models.

Acknowledgements
We would like to thank Janos Balog, Nikolay Gromov and Uwe-Jens Wiese for useful correspondence.

A Spectrum of O(3) rotator in d = 1 at finite lattice spacing
The aim here is to illustrate that the energy levels at a given order of the lattice spacing are polynomials of C I = I(I + 1), not only of I.
The eigenvalues of the transfer matrix are given by where P I (z) are the Legendre polynomials. These are polynomials in 1 − z with coefficients which are polynomials of C I of order I: In the continuum limit one should take → 0. Neglecting non-perturbative, exponentially small terms ∝ e −1/ , one obtains the expansion 11 Note that the summand vanishes for n > I hence the upper limit of the sum could be extended to infinity.
For the energies one obtains Here the coefficients of n are indeed polynomials in C I .

B The sunset diagram for D ∼ 2, 3, 4 for large
The sunset diagram for the susceptibility is (cf. [18] (4.1) and (4.36)) here we have reinstatedˆ , but since we are working withˆ = 1 we have V D = ˆ D−4 = .
For the analogous diagram in the infinite strip one had (cf. eqs. (5.11) and (5.61) in [18] ) (B.2) Further, for 1 one has up to exponentially small corrections (cf. eq. (5.20) in [18]) In Appendix (B.4) of ref. [5] we derived the relation For D ∼ 2 (settingˆ = 1) For the analogous diagram in the infinite strip Using similar methods to that used to derive (B.4) we obtain for D ∼ 2: (B.9) Comparing this with (B.5) and (B.8) requires with values given in eq. (4.45) of ref. [18] W( ) = 2.12506105522294 , for = 1 , For D ∼ 3 Ψ has a singularity: Again using similar methods to that used to obtain (B.4), we can derive the expansion of ∆Ψ for D ∼ 3 for large , and find that the expansion of Ψ( ) large is given by: .

(B.15)
For the purposes of this paper we do not need the value of the constant q 1 appearing in eqs. (B.14) and (B.15). However we can easily obtain its numerical value from the lattice computation of the perturbative coefficients of the moment of inertia given in eq. (7.7) of [20]: These coefficients are simply related to the coefficients in (4.7) through Using eq. (4.11) and the numerical values C The equations for the lowest energy state in a finite volume for the vacuum, 1-and 2-particle sectors We consider here the SU(2)×SU(2) principal chiral model which is equivalent to the O(4) non-linear sigma model. Al. and A. Zamolodchikov [14] proposed a self consistent exact S-matrix for the principal chiral models, which is built from the factors In the equations for I = 0, 1, 2 discussed below the following kernels occur.
In the relations above it is assumed that x is real. Unlike the others, K pp (x) has a pole term −i/(πx). Further, it satisfies the relation (C.6) 13 For d = 4 the analogous estimate for cw at = 2 was correct to 7 digits.
Below we shall consider separately the cases of the ground states in the sectors I = 0, 1, 2. The equations to be solved are clearly presented in ref. [4]; our discussion, however is based on the notations of ref. [6].

C.1 The iterative calculation of the vacuum energy
The finite-volume vacuum energy is defined through a single function A(x), which is obtained by an iterative solution.
The initial function is chosen for convenience as The iterative step with the input function A (k) (x) calculates the output function 14Â(k) (x) through the steps The auxiliary functions 15 are defined below. The reason for usingÂ (k) for the output, is that the input of the next iteration, A (k+1) (x), will be modified to improve the convergence of the iterations. After introducing and r c (x) = r(x) * , the function f (x) is given by Here ' * ' denotes convolution. Due to the singularity of K pp (x) at x = 0 a Principal Value integral appears in (C.9) instead of an ordinary convolution, which is more difficult to deal with numerically. One can circumvent this problem by introducing the function which besides of the PV integral includes also the extra −r c (x) in (C.9). One obtains then a simpler relation with two convolutions only, Finally, the output function of the iteration step iŝ The vacuum energy after k-th iteration is given by (C.13) 14 We denote the output of the iteration withÂ since it could be different from the input of the next iteration. 15 For intermediate quantities like r(x), f (x) we don't write the iteration number k, except for clarity when they appear in the final answer.
To improve convergence one can use the well known trick which takes the input for the next iteration a linear combination with a parameter 0 < α ≤ 1. Note that the choice of the parameter α does not influence the fixed point (provided that the iteration converges), however, it affects the properties of the iteration step, like the speed of convergence or whether the initial function A (0) (x) is in the domain of attraction of the iterative process. One expects that increasing α increases the speed of convergence to the fixed point but decreases the convergence region. The latter could be crucial for small L (i.e. small z) since decreasing L also shrinks the domain of attraction.
For later use, note that the FT of Φ(x) is quite simple, where Θ(p) is the Heaviside step function. The functionΦ(p) is continuous, but its first derivative has a finite jump at p = 0.

C.2 Using discrete FT
The discrete FT in Maple is defined as , (C.16) and its inverse Here x j = jdx , (j = 0, . . . , N − 1) , p k = k dp 2π , (k = 0, . . . , N − 1) For N → ∞ the resolutions 1/dx and 1/dp both in the original x-space (which is here related to rapidity) and in p-space are proportional to √ N . The (rapidity) cutoff is given by X = N dx/2 ∝ √ N . Therefore the N behavior given in (C.18) for constant Q obeys both requirements: for N → ∞ one approaches the continuum limit and the infinite cutoff limit. At this point Q is independent on N but otherwise arbitrary. Its value can be chosen to represent A(x) andÃ(p) by their corresponding discretizations equally well. The number of points in the interval Σ x , the width of the function |A(x)|, is Σ x /dx ≈ √ N Σ x /Q. Similarly, forÃ(p) the number of points in its relevant region in Fourier space is Σ p /dp ≈ √ N Σ p Q/(2π). The value of Q is expected to be optimal when these two numbers are roughly equal, i.e. Q 2 ≈ 2πΣ x /Σ p .
Further for the convolution one has Rewriting (C.11) by performing a discrete Fourier transformation one gets The factor Q appears only in the first term!)

C.3 N -dependence of the sums appearing in FT and the Euler-Bernoulli relation
Assuming that the function g(x) is smooth for x = 0, and it has a jump in the first derivative, according to the Euler-Bernoulli relation one has δg (5) (0)a 6 + . . .
Here a is the lattice spacing representing dx or dp. In our case a 2 ∝ 1/N ; hence under the assumption mentioned above, the corresponding sums have a 1/N expansion with integer powers. One can calculate a quantity for several N values and determine from these the expansion coefficients 16 . The corresponding fits produce a very stable continuum extrapolation, often to ∼ 12 digits.
It seems that it is better to use the naive non-improved sum, (appearing anyhow in the FT) whose 1/N behavior is known rather than an improved sum which has a more complicated (or even unknown) type of approach to the continuum limit. In [6] the authors took α = 1/2. 16 Allowing in the fit odd powers of 1/N 1/2 it turns out indeed that the coefficient of 1/N 3/2 is suppressed by a factor of ∼ 10 −5 .
C.4 The energy of the 1-particle sector For the 1-particle sector the input of an iteration step consists of a function A (k) (x) and a number θ (k) . Within an iteration one has the chain [ , (cf. [6]). The A (k) (x) → r(x) → f (x) part is given by the same expressions we had in (C.8) -(C.11).
In [6] the authors choose a general θ (0) as a starting point of iterations. However, one can verify that θ (k) → 0 for k → ∞. The physical reason is that θ is related to the rapidity of the particle [4]. 17 The state in the I = 1 sector with the smallest energy is the 1-particle state with zero momentum, corresponding to θ = 0. Starting the iteration with θ (0) = 0 all the subsequent θ (k) values remain automatically zero. For this reason, we set θ = 0, simplifying the iteration steps to The output function after one iteration (for θ = 0) is written in the form has the same form as for the vacuum, given by (C.8)-(C.12). (Note that due to the factor [S p (x)] 2 the agreement is only formal.) Introducing the averaging in this case, one obtains for the next input function The energy of the ground state in the 1-particle sector for θ = 0 after the k-th iteration is given by For the case of the 2-particle sector, the input of an iteration step consists of a function A (k) (x) and two numbers, θ (k) 1 and θ (k) 2 . In this case again it is enough to restrict the iteration to the zero total momentum, i.e. θ 1 + θ 2 = 0. (The iteration would drive the system anyhow to satisfy this condition at the FP.) In the rest frame one has θ 2 = −θ 1 = θ/2, where θ is related to the rapidity difference. Because of the symmetry θ → −θ we can restrict the discussion to θ > 0.
In the case of zero total momentum the input of the iteration is [A (k) (x), θ (k) ] and the output is [Â (k) (x),θ (k) ].
Consider the equation z sinh 1 2 πθ + Im ln S 2 0 (θ) + 2φ where we introduced the function (depending implicitly on the input variables) φ(x) = (K m * r)(x) . (C.26) The output value of θ is given by the positive solution of (C.25), The output functionÂ (k) (x) after the k-th iteration is given bŷ where the procedure to obtain f (k) (x) is formally the same as for the ν = 0, 1 cases. Here one can introduce two averaging parameters α A and α θ : and The presence of these two parameters plays an important role in the stability of the iteration and for the speed of convergence. For example, choosing α θ 1 the value of θ barely changes. In our O(4) simulations at small z we found that it is much easier to find a good starting point, (one in the domain of attraction of the fixed point) when one chooses α θ < α A . For our last point, z = LM = 0.0005 we took α A = 0.4, α θ = 0.14.
From the numerical analysis described above we found that θ(z) has a logarithmic dependence for small z, θ(z) ≈ − ln(z)/π + a, where a 1.0.
Using the running coupling g 2 J (z) (cf. (3.34), (3.40)) one obtains a good fit (with an error smaller than 10 −3 ), θ(z) ≈ 1 g 2 J (z) where A = −0.1754, B = −0.2385. As mentioned earlier, the domain of attraction seems to shrink with decreasing z, therefore it is important to choose a proper starting value θ (0) for the iteration. The value obtained by (C.32) can be used for θ (0) . Note, however, that we did not investigate systematically questions related to the domain of attraction.