Localised distributions and criteria for correctness in complex Langevin dynamics
aa r X i v : . [ h e p - l a t ] J un Localised distributions and criteria for correctnessin complex Langevin dynamics
Gert Aarts a ∗ , Pietro Giudice a,b ‡ and Erhard Seiler c § a Department of Physics, College of Science, Swansea UniversitySwansea, United Kingdom b Universit¨at M¨unster, Institut f¨ur Theoretische PhysikM¨unster, Germany ¶ c Max-Planck-Institut f¨ur Physik (Werner-Heisenberg-Institut)M¨unchen, Germany
June 13, 2013
Abstract
Complex Langevin dynamics can solve the sign problem appearing innumerical simulations of theories with a complex action. In order to justifythe procedure, it is important to understand the properties of the real andpositive distribution, which is effectively sampled during the stochasticprocess. In the context of a simple model, we study this distribution bysolving the Fokker-Planck equation as well as by brute force and relatethe results to the recently derived criteria for correctness. We demonstrateanalytically that it is possible that the distribution has support in a stripin the complexified configuration space only, in which case correct resultsare expected. ∗ email: [email protected] ‡ email: [email protected] § email: [email protected] ¶ Present address ontents A.1 Lowest-order solution . . . . . . . . . . . . . . . . . . . . . . . . . 26A.2 First-order correction . . . . . . . . . . . . . . . . . . . . . . . . . 27
Complex Langevin (CL) dynamics [1, 2] provides an approach to circumvent thesign problem in numerical simulations of lattice field theories with a complexBoltzmann weight, since it does not rely on importance sampling. In recentyears a number of stimulating results has been obtained in the context of nonzerochemical potential, in both lower and four-dimensional field theories with a severesign problem in the thermodynamic limit [3–8] (for two recent reviews, see e.g.Refs. [9,10]). However, as has been known since shortly after its inception, correctresults are not guaranteed [11–16]. This calls for an improved understanding,relying on the combination of analytical and numerical insight. In the recent past,the important role played by the properties of the real and positive probabilitydistribution in the complexified configuration space, which is effectively sampledduring the Langevin process, has been clarified [17, 18]. An important conclusionwas that this distribution should be sufficiently localised in order for CL to yieldvalid results. Importantly, this insight has recently also led to promising resultsin nonabelian gauge theories, with the implementation of SL( N, C ) gauge cooling[8, 10].The distribution in the complexified configuration space is a solution of theFokker-Planck equation (FPE) associated with the CL process. However, incontrast to the case of real Langevin dynamics, no generic solutions of this FPE2re known (see e.g. Ref. [19]). In fact, even in special cases only a few resultsare available [11, 17, 20, 21]. In Refs. [17, 18] this problem was addressed in aconstructive manner by deriving a set of criteria for correctness, which have tobe satisfied in order for CL to be reliable. These criteria reflect properties of thedistribution and, importantly, can easily be measured numerically during a CLsimulation, also in the case of multi-dimensional models and field theories [6].A widely used toy model to understand CL is the simple integral Z = Z ∞−∞ dx e − S , S = 12 σx + 14 λx , (1.1)where the parameters in the action are complex-valued. This model has beenstudied shortly after CL was introduced [11, 22, 23], but no complete solution wasgiven. As we will see below, its structure, with complex σ , is relevant for therelativistic Bose gas at nonzero chemical potential [4, 20]. Recently, a variant ofthis model (with σ = 0 and λ complex) was studied by Duncan and Niedermaier[21]: in particular they constructed the solution of the FPE, using an expansionin terms of Hermite functions. They considered the case of “complex noise”, inwhich both the real and imaginary parts of the complexified variables are subjectto stochastic kicks. Unfortunately, it has been shown in the past that genericallycomplex noise may not be a good idea, since it leads to broad distributions inthe imaginary direction and hence incorrect results [17, 18]. This was indeedconfirmed in Ref. [21].In this paper we aim to combine the insights that can be distilled from thecriteria for correctness discussed above with the explicit solution of the FPE,adapting the method employed in Ref. [21] to the model (1.1). The paper isorganised as follows. In Sec. 2 we discuss CL and the criteria for correctness.To keep the paper sufficiently accessible, we first briefly review how to arriveat the criteria for correctness and subsequently present numerical results, forboth real and complex noise. In Sec. 3 we study the probability distributionin the complexified configuration space, by solving the FPE directly as well asby a brute-force construction using the CL simulation, again for complex andreal noise (the latter was not considered in Ref. [21]). In Sec. 4 we combine ourfindings concerning the distribution and the criteria for correctness, and providea complete characterisation of the dynamics. Sec. 5 contains the conclusion.Finally, in order to see whether the structure found numerically can be understoodanalytically, a perturbative analysis of the FPE is given in Appendix A. We consider the partition function (1.1). We take λ real and positive, so thatthe integral exists, while σ is taken complex. Analytical results are available: a3irect evaluation of the integral yields Z = r ξσ e ξ K − ( ξ ) , (2.1)where ξ = σ / (8 λ ) and K p ( ξ ) is the modified Bessel function of the second kind.Moments h x n i can be obtained by differentiating with respect to σ . Odd momentsvanish.The aim is to evaluate expectation values numerically, by solving a CL process.We start from the Langevin equation,˙ z = − ∂ z S ( z ) + η, (2.2)where the dot denotes differentiating with respect to the Langevin time t and the(Gaussian) noise satisfies h η ( t ) η ( t ′ ) i = 2 δ ( t − t ′ ) . (2.3)After complexification, z = x + iy, η = η R + iη I , σ = A + iB, (2.4)the CL equations read˙ x = K x ( x, y ) + η R , ˙ y = K y ( x, y ) + η I , (2.5)with the drift terms K x ≡ − Re ∂ z S ( z ) = − Ax + By − λx (cid:0) x − y (cid:1) , (2.6) K y ≡ − Im ∂ z S ( z ) = − Ay − Bx − λy (cid:0) x − y (cid:1) . (2.7)The form of the drift terms is similar as in the Bose gas, after a reduction to asingle momentum mode [20].The normalisation of the real and imaginary noise components follows fromEq. (2.3) and is given by h η R ( t ) η R ( t ′ ) i = 2 N R δ ( t − t ′ ) , h η I ( t ) η I ( t ′ ) i = 2 N I δ ( t − t ′ ) , h η R ( t ) η I ( t ′ ) i = 0 , (2.8)with N R − N I = 1. Here N I ≥ N I , but inpractice they are not. Real noise amounts to N I = 0.Expectation values are obtained by averaging over the noise. After this aver-aging, holomorphic observables evolve according to h O i P ( t ) = Z dxdy P ( x, y ; t ) O ( x + iy ) , (2.9)4here the distribution P ( x, y ; t ) satisfies the FPE˙ P ( x, y ; t ) = L T P ( x, y ; t ) , (2.10)with the FP operator L T = ∂ x ( N R ∂ x − K x ) + ∂ y ( N I ∂ y − K y ) . (2.11)In order to justify the approach, we also consider expectation values with respectto a complex weight ρ ( x, t ), h O i ρ ( t ) = Z dx ρ ( x, t ) O ( x ) , (2.12)which satisfies its (complex) FPE˙ ρ ( x, t ) = L T ρ ( x, t ) , L T = ∂ x [ ∂ x + ( ∂ x S ( x ))] . (2.13)This equation has a simple stationary solution, ρ ( x ) ∼ e − S ( x ) , which is the desiredweight.The task is now to show that the two expectation values h O i P ( t ) and h O i ρ ( t ) are equal, h O i P ( t ) = h O i ρ ( t ) , (2.14)at least in the limit of large t , making use of the respective FPEs and the Cauchy-Riemann (CR) equations [17, 18]. Here it is essential that only holomorphicobservables are considered, which evolve according to ∂ t O ( z, t ) = ˜ LO ( z, t ) , (2.15)with the Langevin operator ˜ L = [ ∂ z − ( ∂ z S ( z ))] ∂ z . (2.16)We note that for holomorphic observables, ˜ L = L , where L is the transpose of L T introduced above. The equivalence (2.14) can indeed be shown, as discussedin detail in Refs. [17, 18], provided that integration by parts in y is allowed,without the presence of boundary terms at infinity. This construction involvesthe products P ( x, y ; t ) O ( x + iy ) for ‘all’ observables O ( x ), and hence it puts severeconstraints on the decay of the distribution at infinity. This will indeed be shownto be crucial below.From now on we consider only the equilibrium distribution P ( x, y ), assumingthat it exists, and hence drop the t dependence. In the large t limit, the equiva-lence (2.14) can then be expressed in terms of the criteria for correctness [17, 18] C O ≡ D ˜ LO ( z ) E = 0 , (2.17)5 n -0.3-0.1500.150.3 Re
0. Here we focus on σ = 1 + i and λ = 1. In Fig. 1CL results are shown for the real and imaginary parts of the observables n h z n i and for the criteria for correctness C n = n h ˜ Lz n i , for n = 2 , , ,
8. The figureshows the result for real noise, N I = 0: all expectation values agree with theexact result, denoted with the horizontal lines, and the criteria for correctnessare all consistent with 0, as it should be.In Fig. 2 we show how the observables and the criteria for correctness dependon the amount of complex noise. In the top figures we see that for small N I theobservables with n = 2 , N I they start to deviate. Perhaps surprisingly, the lowest-order criterium C is consistent with 0 for all N I shown. This implies that even though h z i and h z i have converged to the wrong result at larger N I , this occurs in such a way7hat the condition (2.22), i.e. the corresponding SD equation, is still satisfied.The possibility of multiple solutions to the SD equations when solving CL hasbeen observed earlier in Ref. [13] (see also Refs. [25, 26]).In order to detect problems, it is necessary to consider higher moments. InFig. 2 (below), we observe that for small N I the observables (with n ≥
6) and thecriteria (with n ≥
4) are only marginally consistent with the expected results,while for larger N I they suffer from large fluctuations and can no longer be sensiblydetermined. According to the analytical justification [17, 18], this implies thatthe results from CL cannot be trusted in the presence of complex noise. Belowwe give an interpretation of this in terms of the properties of the probabilitydistribution. For now we tentatively conclude that, if we assume that the largefluctuations reflect the slow decay of the distribution in the imaginary direction, P ( x, y ) should decay as 1 / | y | α , with 5 . α .
7, which will indeed be confirmedbelow.
A crucial role in the justification of the method is played by the equilibriumdistribution P ( x, y ) in the complexified space. In Refs. [17, 18] it was shown indetail that for CL to give correct results, it is necessary that the product of thedistribution and a suitable basis of observables drops off fast enough in the imag-inary direction. This condition can be translated into the criteria for correctness,as discussed above. Unfortunately the Fokker-Planck equation, satisfied by thedistribution, cannot be solved easily, except in the case of a noninteracting model( λ = 0), see Appendix A.In this section we study the distribution following two approaches. Firstly, itis possible to collect histograms of the (partially integrated) distribution duringthe CL evolution. Note that very long runs are required, in order to samplethe configuration space properly. Here we will in particular be interested in thepartially integrated distributions P x ( x ) = Z ∞−∞ dy P ( x, y ) , P y ( y ) = Z ∞−∞ dx P ( x, y ) . (3.1)We note that this approach can easily be extended to multi-dimensional integralsand field theories. We refer to this as the brute force method.Secondly, for the zero-dimensional model we consider here, it is possible toexpand the distribution in terms of a truncated set of basis functions and solvethe resulting matrix problem numerically, following Duncan and Niedermaier [21].We discuss this approach in the next subsection.8 .1 Solving the Fokker-Planck equation We consider the eigenvalue problem − L T P κ ( x, y ) = κP κ ( x, y ) , (3.2)where the FP operator L T was given in Eq. (2.11) and takes the explicit form L T = N R ∂ x + ( Ax − By ) ∂ x + N I ∂ y + ( Ay + Bx ) ∂ y + 2 A + λ (cid:0) x − xy (cid:1) ∂ x + λ (cid:0) x y − y (cid:1) ∂ y + 6 λ (cid:0) x − y (cid:1) . (3.3)We denote the eigenvalues of − L T with κ and the eigenfunctions with P κ ( x, y ).If there is a unique ground state P with eigenvalue κ = 0, and for all othereigenvalues Re κ >
0, the time-dependent distribution can be written as P ( x, y ; t ) = P ( x, y ) + X κ =0 e − κt P κ ( x, y ) , (3.4)and the equilibrium distribution is given by P ( x, y ). In the CL simulationswe observe convergence to well-defined expectation values (at least for the lowmoments, n = 2 ,
4) and hence we are certain that an equilibrium distributionexists.In order to solve the eigenvalue problem, we follow closely Ref. [21]. The FPoperator is invariant under x → − x, y → − y , which implies that eigenfunctionshave a definite parity, P κ ( x, y ) = ± P κ ( − x, − y ). The ground state is expected tosatisfy P ( x, y ) = P ( − x, − y ), such that observables of the type h ( x + iy ) n i , with n odd, vanish. If P κ is an eigenfunction of L T with eigenvalue κ , then so is P ∗ κ with eigenvalue κ ∗ . It is expected that P is real.In Ref. [21] P ( x, y ) was doubly expanded in a basis of Hermite functions, i.e. P ( x, y ) = N H − X k =0 N H − X l =0 c kl H k (cid:0) √ wx (cid:1) H l (cid:0) √ wy (cid:1) , (3.5)where ω is a variational parameter appearing in the harmonic oscillator eigen-functions, and N H indicates the number of Hermite functions included in thetruncated basis. The coefficients c kl have to be determined.In order to do so, we introduce creation and annihilation operators, satisfying[ a, a † ] = [ b, b † ] = 1 , (3.6)and write x = 1 √ ω (cid:0) a + a † (cid:1) , p x = − i∂ x = i r ω (cid:0) a † − a (cid:1) , (3.7) y = 1 √ ω (cid:0) b + b † (cid:1) , p y = − i∂ y = i r ω (cid:0) b † − b (cid:1) . (3.8)9n terms of these, − L T reads − L T = N R p x + N I p y − i ( Ax − By ) p x − i ( Ay + Bx ) p y − A − λ (cid:0) x − y (cid:1) + λ ω [ X ( x, y ) − X ( y, x )] , (3.9)with the quartic terms X ( x, y ) = − iω (cid:0) x − xy (cid:1) p x , X ( y, x ) = − iω (cid:0) y − x y (cid:1) p y . (3.10)Note that X is independent of ω . Finally, in terms of the creation/annihilationoperators, the FP operator reads − ω L T = − N R (cid:0) a † + a − a † a − (cid:1) − N I (cid:0) b † + b − b † b − (cid:1) + ¯ A (cid:0) a † − a + b † − b + 2 (cid:1) + 2 ¯ B (cid:0) b † a − a † b (cid:1) − A − ¯ λ (cid:2)(cid:0) a † + a + 2 a † a (cid:1) − (cid:0) b † + b + 2 b † b (cid:1)(cid:3) + ¯ λ
12 [ X ( a, b ) − X ( b, a )] , (3.11)where X ( a, b ) = (cid:0) a + a † (cid:1) (cid:0) a † − a (cid:1) − (cid:0) a † + a (cid:1) (cid:0) a † − a (cid:1) (cid:0) b † + b (cid:1) , (3.12)and we introduced the rescaled parameters,¯ A = Aω , ¯ B = Bω , ¯ λ = 6 λω . (3.13)In Ref. [21], where A = B = 0, ω was chosen to be proportional to √ λ , and noadjustable parameters were left on the RHS of Eq. (3.11). As we see below, thereis a great advantage in keeping ω arbitrary.We can now compute the matrix elements with respect to the Hermite func-tions, using the notation | mn i = 1 √ m ! n ! a † m b † n | i , a | i = b | i = 0 , (3.14)where H m ( √ ωx ) = h x | m i , H n ( √ ωy ) = h y | n i . (3.15)The matrix elements are − ω h kl | L T | mn i = (cid:2)(cid:0) N R − ¯ λ (cid:1) (2 m + 1) + (cid:0) N I + ¯ λ (cid:1) (2 n + 1) − A (cid:3) δ k,m δ l,n − (cid:2)(cid:0) N R + ¯ λ − ¯ A (cid:1) f km δ k,m +2 + (cid:0) N R + ¯ λ + ¯ A (cid:1) f mk δ k,m − (cid:3) δ l,n − (cid:2)(cid:0) N I − ¯ λ − ¯ A (cid:1) f ln δ l,n +2 + (cid:0) N I − ¯ λ + ¯ A (cid:1) f nl δ l,n − (cid:3) δ k,m
102 ¯ B (cid:16) √ mlδ k,m − δ l,n +1 − √ knδ k,m +1 δ l,n − (cid:17) + ¯ λ
12 [ X kl,mn − X lk,nm ] , (3.16)with X kl,mn = h f km δ k,m +4 + (2 m + 3 − n ) f km δ k,m +2 + 6( m − n ) δ k,m − (2 m − − n ) f mk δ k,m − − f mk δ k,m − i δ l,n − f km δ k,m +2 − f mk δ k,m − + δ k,m ] [ f ln δ l,n +2 + f nl δ l,n − ] , (3.17)and f km = r k ! m ! . (3.18)Following Ref. [21], the double indices k, l and m, n (all taking values from 0to N H −
1) are converted into single ones, via i = kN H + l + 1 , j = mN H + n + 1 , (3.19)and the inverse k i = ( i − − mod( i − , N H )) /N H , l i = mod( i − , N H ) , (3.20) m j = ( j − − mod( j − , N H )) /N H , n j = mod( j − , N H ) , (3.21)with i, j = 1 , . . . , N H . The matrix elements are denoted as L Tij = h kl | L T | mn i ,and the eigenvalue problem is written as − L Tij v ( κ ) j = κv ( κ ) i . (3.22)We have solved this matrix problem with a FORTRAN90 code using subroutinesprovided by the LAPACK library [27]. Since the matrix size is N H × N H , thereis an upper limit of what is practically feasible. For the maximal number ofHermite functions we have considered, N H = 150, the numerical computationtakes around 36 hours on a standard work station. Convergence can be testedby increasing N H and varying ω (see the detailed discussion below). Consideringthe eigenvalue at (or closest to) 0, the distribution P ( x, y ) can be reconstructedfrom the corresponding eigenvector, as P ( x, y ) = N H X i =1 v (0) i H k i (cid:0) √ wx (cid:1) H l i (cid:0) √ wy (cid:1) . (3.23)Below we drop the subscript ‘0’. 11 I ω N I and ω used, with σ = 1 + i , λ = 1, and 30 ≤ N H ≤ We start with the case of complex noise. The parameters in the action aretaken as σ = 1 + i and λ = 1, and we consider a basis with 30 ≤ N H ≤ ω we used are listed in Table 1. In the limitof large N H the results are expected to be independent of the value of ω . Inpractice however, we find that for finite N H the parameter ω plays the role of atuning parameter: in particular, when ω is too small, there are eigenvalues witha negative real part. This becomes more prominent as N I is reduced, see below.Obviously, in this application this would mean that the FP evolution would notthermalise and display runaway behaviour. Since the CL evolution thermalises(and is obviously independent of the choice of ω ), we expect the real parts ofall eigenvalues to be nonnegative. When the value of ω is increased, we observethat the eigenvalues with a real negative part move into the positive half-planeand the spectrum around the origin converges. Convergence can also be seenby studying the reconstructed probability distribution P ( x, y ), using Eq. (3.23).Interestingly, we always find an eigenvalue consistent with 0. When ω is increasedeven more, convergence properties worsen again. We find therefore that there isan ω interval for which:1. there is an eigenvalue consistent with 0;2. the other eigenvalues are in the right half-plane;3. the reconstructed ground state distribution is stable under variation of N H and ω .The ω interval depends on the parameters and is pushed to larger values as N I is reduced. We have not found a special role for the ω value used in Ref. [21],namely ω = √ λ (in our conventions).We first consider N I = 1, as in Ref. [21]. The smallest 15 eigenvalues areshown in Fig. 3 (left), for several values of ω . For the ω values shown here,all eigenvalues are in the right half-plane and the spectrum around the originis to a good extent independent of ω . The reconstructed distribution P ( x, y ),obtained using the eigenvector corresponding to the eigenvalue at (or closest to)the origin, is shown in Fig. 4 (top). We find a smooth distribution with a doublepeak structure, similar as in Ref. [21]. 12 Re κ -10-50510 I m κ ω = 1.5ω = 2ω = 5ω = 10σ = 1 +i, λ = 1, N I = 1, N H = 150 Re κ -10-50510 I m κ ω = 4ω = 8ω = 12ω = 16σ = 1 +i, λ = 1, N I = 0.01, N H = 150 Figure 3: Eigenvalues of the FP operator − L T for complex noise, with N I = 1(left) and 0.01 (right), magnified around the smallest eigenvalues, for variousvalues of ω , at σ = 1 + i , λ = 1, and N H = 150.Next we reduce the amount of complex noise and consider N I = 0 .
01. Thespectrum is shown in Fig. 3 (right) and the reconstructed distribution in Fig. 4(below). The findings are similar as with N I = 1, but ω has to be increased morein order to find convergence and even then the larger eigenvalues are hard toestablish. The distribution has again two peaks, which are now more pronouncedand rotated in the xy -plane. We note the symmetry P ( − x, − y ) = P ( x, y ). Im-portantly, the distribution is more squeezed in the y direction and the mainfeatures are contained in the interval − . < y < . P x ( x ) and P y ( y ), see Eq. (3.1), on a logarithmicscale, for the case of N I = 1. Besides presenting results for various ω values,we also show the histogram obtained during a CL simulation. We observe anacceptable agreement between the CL results and the solution of the FPE for ω ∼ . ,
2, down to a relative size of 10 − , after which the FP solution can nolonger cope. We interpret this as a manifestation of the truncation. When ω istaken too large, the disagreement occurs for smaller values of x and y .The distributions do not go to zero rapidly but decay as a power, which isclearly visible on a log-log plot. In Fig. 6 we show the distributions multipliedby x k and y k respectively, for k = 4 . ,
5, and 5 .
2, using the CL data. At large | x | and | y | , we observe a power decay with power 5, i.e. P x ( x ) ∼ | x | , P y ( y ) ∼ | y | . (3.24)This suggests that the distribution decays as P ( x, y ) ∼ x + y ) , (3.25)13igure 4: Distribution P ( x, y ) in the xy -plane for complex noise, with N I = 1(top, with ω = 1 .
5) and 0.01 (bottom, ω = 8). Other parameters as in Fig. 3.which we have verified by studying the decay of P r ( r ) = Z π dφ rP ( r cos φ, r sin φ ) , (3.26)which indeed decays as 1 /r . We note that this power decay is in agreementwith the conclusions from the moments above: h z i and h z i are well defined andcan be numerically determined without any problems, while the higher momentsdiverge, which in the CL simulation is reflected in large fluctuations. We now turn to the case where CL appears to work well, i.e. with real noise( N I = 0). The eigenvalues are shown in Fig. 7 for a number of ω values. For14 x P x ( x ) CL ω = 1.5ω = 2ω = 5ω = 10σ = 1 +i, λ = 1, N I = 1, N H =150 y P y ( y ) CL ω = 1.5ω = 2ω = 5ω = 10σ = 1 +i, λ = 1, N I = 1, N H = 150 Figure 5: Partially integrated distributions P x ( x ) (left) and P y ( y ) (right) fordifferent values of ω with complex noise, N I = 1. Other parameters as in Fig. 3.In both cases the noisy (black) data was obtained by a CL simulation. x x k P x ( x ) k = 5.2 k = 5 k = 4.8 σ = 1 +i, λ = 1, N I = 1 y y k P y ( y ) k = 5.2 k = 5 k = 4.8 σ = 1 +i, λ = 1, N I = 1 Figure 6: As above, for x k P x ( x ) and y k P y ( y ) with k = 4 . , , .
2, using the CLdata. The dotted horizontal line is meant to guide the eye. ω < P x and P y under variation of N H and ω , and also compare thosewith the histograms obtained with CL. The results are shown in Fig. 8 for P y ( y )(top) and P x ( x ) (bottom). In the case of P y , convergence as N H is increased isclearly visible (top, left). We note that for the largest N H values the distributionagrees with the result obtained by direct Langevin simulation, indicated withthe black line. The distribution is very well localised and appears to drop to 0around y = 0 .
28. We come back to this below. Convergence as ω is increasedis demonstrated in Fig. 8 (top, right) and we observe that a large value of ω isrequired, ω ∼
50. It is of course expected that the chosen value of ω eventually15 Re κ -10-50510 I m κ ω = 4 ω = 10 ω = 40ω = 50 ω = 60σ = 1 +i, λ = 1, N I = 0, N H = 150 Figure 7: As in Fig. 3, for real noise ( N I = 0).becomes irrelevant, but for finite N H keeping ω as a tuning parameter is essential.The distribution P x ( x ) is shown in Fig. 8 (below) as a function of x (left)and x (right), on a logarithmic scale. In contrast to the case of complex noise,we now find an exponential rather than a power decay. Results from FPE agreewith the CL histogram, independently of the value of ω in this case, but onlydown to a relative size of 10 − ; varying ω does not help in this case (increasing N H probably will). From the CL result, we see that the distribution falls off as P x ( x ) ∼ e − ax , a ∼ . . (3.27)Naively this behaviour can be expected, since for large | x | the original weightbehaves as ∼ exp ( − λx / λ/ .
25. Interestingly this seems to be understandable from aperturbative analysis, see Appendix A.The reconstructed distribution is shown in Fig. 9. This distribution has sim-ilar characteristics as at N I = 0 .
01, except that the two peaks are now verypronounced and the saddle around the origin is much deeper. The peaks liemostly in the y direction and they are therefore clearly visible in P y ( y ). The dis-tribution is squeezed even more than before and its main support is in the region − . < y < .
3. The ripples visible for larger y values are an artefact of thetruncation. In fact, in the next section we will demonstrate that the distributionis strictly 0 when | y | > . σ = 1 + i and λ = 1) thedecay in the case of real noise is manifestly different compared to complex noise.In the latter we found a power decay, resulting in ill-defined moments h z n i when n >
4, while here we find exponential decay in the x direction and, as we willsee below, in the y direction support only inside a strip. As a result there is no16 y P y ( y ) CL N H = N H = N H = N H = N H = N H = N H = σ = 1+ i , λ = 1, N I = 0, ω = 50 y P y ( y ) CL ω = 4 ω = 10 ω = 40ω = 50 ω = 60σ = 1+ i, λ = 1, N I = 0, N H = 150 x P x ( x ) CL ω = 4 ω = 10 ω = 40 ω = 50 ω = 60 σ = 1+ i, λ = 1, N I = 0, N H = 150 x P x ( x ) CL ω = 4 ω = 10 ω = 40 ω = 50 ω = 60 σ = 1+ i, λ = 1, N I = 0, N H = 150 Figure 8: Above: Partially integrated distribution P y ( y ) for several values of N H and ω = 50 (left) and several values of ω and N H = 150 (right). Below: Partiallyintegrated distribution P x ( x ) on a logarithmic scale as a function of x (left) and x (right) for several values of ω and N H = 150. The dotted line on the RHSrepresents P x ( x ) ∼ exp( − ax ) with a = 0 . From the solution of the FPE and the CL process, we conclude tentatively thatfor real noise the distribution is localised in the y direction and has support in astrip around the origin only, with − . . y . .
3. This conclusion can be mademore precise by studying the classical flow diagram and properties of the FPE.This analysis can also be used to find parameter values for which CL breaks downfor real noise (see Sec. 4.3). 17igure 9: Distribution P ( x, y ) in the xy -plane for real noise ( N I = 0) at σ = 1+ i and λ = 1, using N H = 150 and ω = 50. The classical flow diagram is shown in Fig. 10, for σ = 1 + i and λ = 1. Weshow the direction of the classical force by an arrow pointing in the direction( K x ( x, y ) , K y ( x, y )). The arrows are normalised to have the same length. Theclassical force is of course independent of N I . There are three fixed points, where K x = K y = 0: an attractive point at the origin and two repulsive fixed points,determined by σ + λz = 0, or x − y = − A , xy = − B λ , (4.1)yielding ( x, y ) = ( ± . , ∓ .
10) in this case. The flow is directed towards theorigin, provided that | y | is not too large. This can be made more precise bystudying where K y ( x, y ) changes sign. We find that K y ( x, y ) = 0 at y p ( x ) = 2 (cid:18) B λ + x (cid:19) cos (cid:18) α + pπ (cid:19) , p = 1 , , , (4.2)where α = − arctan λAx "(cid:18) B λ + x (cid:19) − (cid:18) Ax λ (cid:19) + π Θ( x ) , (4.3)with Θ( x ) the step function. These lines are indicated in the classical flow dia-gram with full lines. For the parameter values we consider here, the upper andlower curves have extrema at x = ± . y = ∓ . x = ± . y = ∓ . x -2-1012 y Figure 10: Classical flow in the xy -plane, for σ = 1 + i and λ = 1. The at-tractive/repulsive fixed points are indicated with the open/filled circles. The fulllines indicate where K y ( x, y ) = 0. The horizontal dashed lines indicate the stripin which the CL process takes place in the case of real noise.We now realise that along the horizontal dashed lines, which are determinedby the extrema of the centre curve where K y = 0 ( y = ± . − . < y < . . (4.4)This is consistent with the conclusions drawn above from the histograms and theFPE solution of the distribution P ( x, y ). In the presence of complex noise, thisconclusion no longer holds and the entire xy -plane can be explored. It is possible to make the argument based on classical flow presented above rig-orous and show directly from the FPE that the equilibrium distribution P ( x, y )19s strictly zero in strips in the xy -plane, assuming sufficient decay, i.e. K x,y ( x, y ) P ( x, y ) → x and/or y → ±∞ . To achieve this, we note that the FPE takes the form ofa conservation law, i.e.,˙ P ( x, y ; t ) = ∂ x J x ( x, y ; t ) + ∂ y J y ( x, y ; t ) , (4.6)with J x = ( N R ∂ x − K x ) P, J y = ( N I ∂ y − K y ) P, (4.7)which allows us to consider the charge, Q ( y, t ) = Z ∞−∞ dx J y ( x, y ; t ) . (4.8)Specialising now to the equilibrium distribution (and hence dropping the t de-pendence), we find that Q ( y ) is independent of y , provided that the product ofthe drift K x ( x, y ) and the distribution P ( x, y ) drops to zero at large | x | , since ∂ y Q ( y ) = Z ∞−∞ dx ∂ y J y ( x, y ) = − Z ∞−∞ dx ∂ x J x ( x, y ) = − J x ( x, y ) (cid:12)(cid:12)(cid:12) ∞ x = −∞ = 0 . (4.9)We note that the required condition is always satisfied in our case, even in thecase of the power decay. Since Q ( y ) vanishes as y → ±∞ (because J y ( x, y ) does,again relying on the sufficient decay), we find that Q ( y ) = Z ∞−∞ dx ( N I ∂ y − K y ( x, y )) P ( x, y ) = 0 . (4.10)For real noise, this yields therefore the condition Q ( y ) = Z ∞−∞ dx K y ( x, y ) P ( x, y ) = 0 , (4.11)for all y . Since P ( x, y ) is nonnegative, this condition allows us to derive thefollowing useful property: if K y ( x, y ) has a definite sign as a function of x forgiven y , P ( x, y ) has to vanish for this y value. As a function of x , K y ( x, y ) is aparabola with an extremum at x = − B λy (4.12)and a curvature of 6 λy . The value at the extremum is given by F ( y ) ≡ K y ( x , y ) = − λy "(cid:18) y − A λ (cid:19) − A − B λ . (4.13)20 − xy y + −y + −y − Figure 11: The distribution P ( x, y ) is strictly zero in the strips bounded by ± y − and ± y + , provided that 3 A > B and N I = 0.Consider now the case that y is positive (negative). In that case, when F ( y ) > F ( y ) < K y ( x, y ) is strictly positive (negative) and hence P ( x, y ) has tovanish. The zeroes of F ( y ) are given by y ± = A λ ± r − B A ! , (4.14)provided that 3 A − B >
0. Inspection shows that F ( y ) > y − < y < y + and F ( y ) < − y + < y < − y − : hence for these y values, P ( x, y ) = 0.When 3 A − B < F ( y ) has no zeroes and F ( y ) and y have opposite signs. Inthat case, K y ( x, y ) has no definite sign and the reasoning cannot be followed.To summarise, we find the following:1. when 3 A > B , P ( x, y ) = 0 when y − < y < y , as illustrated in Fig. 11;2. when B > A , there are no restrictions on P ( x, y ).In the first case the distribution can in principle be nonzero in the outerregion, y > y . However, once the process is in the inner strip determined by y < y − , it will not be able to leave this strip, due to the nature of the driftterms. Hence there is no objection to putting the distribution to zero also when y > y . We conclude therefore that the equilibrium distribution has support inthe strip determined by y < y − only, in agreement with the reasoning above.Note that P ( x, y ) is therefore a nonanalytic function of y . Of course the value of y − agrees with the boundary determined in the example in the previous section,i.e. with the position of the dashed lines in Fig. 10, as it should be.For vanishing B , the action is real and the distribution is (for real noise)strictly localised on the real axis, y = 0. For small B , the width of the allowed21 y P y ( y ) N I = 0 N I = 0.01 N I = 1 σ = 1+ i, λ = 1 y P y ( y ) B =
B =
B =
B =
B =
B = 2.2B =
B =
B =
B = B = 2.7 B = 2.8 B = 2.9 B = 3.0 σ = 1+ iB, λ = 1, N I = 0 Figure 12: Distribution P y ( y ) for different values of N I (left, with B = 1) and B (right, with N I = 0) at σ = 1 + iB and λ = 1, obtained with CL. On the leftthe vertical line at x = 0 . B = 1 .
7. For larger B values, there is no longer a boundary.region around y = 0 is nonzero and set by y − ∼ B λA . (4.15)Hence increasing the amount of complexity by increasing B results in a broad-ening of the distribution with a width ∼ B . The importance of this controlledincrease has been emphasised earlier in Ref. [28]. The argument presented above breaks down in the presence of complex noise. Inthat case, the process is pushed out in the y direction and the repulsive fixedpoints come into play. Once the repulsive fixed point is crossed, large excursionsin the y direction take place and the distribution is no longer localised. Whenthe amount of complex noise is small, it takes time to notice this, but eventuallyit will happen. There are therefore no strips for complex noise, which also followsfrom the formal derivation above. This is demonstrated in Fig. 12 (left), where P y ( y ) is shown for the values of N I considered above. As shown above, this leadsto power decay, P y ( y ) ∼ / | y | .Interestingly, the derivation above demonstrates that strips are only presentwhen 3 A > B . For larger B values, one may therefore expect a breakdown ofCL with real noise, similar as with complex noise. This is indeed what happens.The distribution P y ( y ) as B is increased is shown in Fig. 12 (right), for realnoise. Note the similarity with the figure on the left. The delocalisation has adetrimental effect on the results of the CL process. This is demonstrated in Fig.22 .75 2 2.25 2.5 2.75 3 B -0.002-0.00100.0010.002 < z n > / n - e x ac t Re
04 and y = ∓ . P ( x, y ) = c ( x + y ) α (4.16)at large x and y , where we found numerically that the power α is consistent with3. Substituting this Ansatz in the FPE (2.10), we find, after some algebra andthe removal of common factors, that α x − y + 2 α ( N R x + N I y )( x + y ) + A (1 − α ) + λ (3 − α )( x − y ) = 0 . (4.17)At large x and/or y the final term dominates: requiring that this term vanishesyields indeed α = 3. This construction assumes that the behaviour at large23 .1 1 10 x, y P x ( x ) P y ( y ) σ = 1+3 i, λ = 1, N I = 0 Figure 14: Above: Distribution P ( x, y ) obtained from the FPE, with ω = 6 and N H = 150. Below: Partially integrated distribution P x ( x ) and P y ( y ) on a log-logscale, obtained from CL. The dotted line shows a power law 1 /x . The verticallines indicate the x and y coordinate of the repulsive fixed point. In both plots, σ = 1 + 3 i , λ = 1, and N I = 0 (real noise).distance is approximately rotationally invariant in the xy -plane and that thereare no preferred directions, which would invalidate the Ansatz and the powercounting above. Based on our numerical evidence, this seems to be the case.We note that the final term in Eq. (4.17) is independent of σ = A + iB and N I ; hence the decay at large distance is independent of the parameters in theaction and of the amount of complex noise. We also note that B has disappearedfrom Eq. (4.17): the reason is that B breaks the invariance under x → − x andindependently y → − y , while the Ansatz is invariant under those.The conclusion is therefore that the decay at large x and y is universal. Ofcourse the presence of complex noise and/or a large value of B > A is essentialin catalysing large excursions, which lead to the power decay. Notably, the powerdecay appears to be unavoidable unless its appearance is strictly forbidden, as in24he case of the strips for real noise and B < A . In order to justify the results obtained with complex Langevin dynamics, it isnecessary that the probability distribution is sufficiently localised in the complex-ified configuration space. Here we have studied properties of this distribution viaa number of methods, in the case of a simple model. Using the insights gatheredfrom classical flow, histograms obtained during the CL process, the criteria forcorrectness and the explicit solution of the FPE, a complete characterisation ofthe distribution can be given.In the case of real noise and provided that B < A , where σ = A + iB , wefound that the distribution is strictly localised, i.e. it has support in a strip in theconfiguration space only, with exponential decay in the real direction. In this caseall moments are well-defined and, relying on the analytical proof of the method,correct results are expected. We also found that the criteria for correctness aresatisfied. In contrast, when the noise is complex or when B > A , the entireconfiguration space is explored. Large excursions are possible due to the presenceof repulsive fixed points and the decay of the distribution changes dramatically.We found strong indications that for large | x | and | y | , the distribution decays asa power, according to P ( x, y ) ∼ x + y ) . (5.1)A consequence of this slow decay is that higher moments are no longer well-defined. As a result, these and the criteria for correctness suffer from largefluctuations during the CL process, an important signal of failure. Here it is im-portant to emphasise that the inclusion of higher moments is essential to observethe breakdown.In this model the FPE can be solved explicitly, via an expansion in a truncatedset of basis functions. However, it is still a nontrivial problem and perhaps thebest way to find the distribution is by brute force, i.e. during the CL simulation.This also has the benefit of being applicable to higher dimensional models. Inthe case of the localised distribution in the strip, the used basis set may not bethe one that is best adapted to the problem and, in hindsight, once it has beendemonstrated that the distribution has support in a strip only, a more suitablebasis can be used. This would however limit the generality of the approach.As an outlook, we note that in the more realistic cases of multi-dimensionalmodels and field theories, the luxury of solving the FPE is typically not avail-able. However, we have demonstrated that the essential insight can already beobtained from a combination of histograms of partially integrated distributionsand the criteria for correctness, which gives a consistent picture of the dynam-ics. These tools are readily available in field theory. Finally, our conclusions are25lso immediately applicable to nonabelian SU( N ) gauge theories, for which gaugecooling provides a means to control the distribution in SL( N, C ), a possibility notpresent in simpler models. Acknowledgments
We thank Denes Sexty and Ion-Olimpiu Stamatescu for discussion. This work issupported by STFC and the Royal Society.
A Perturbative solution of the FP equation
In order to understand the numerical solution for the distribution P ( x, y ) foundabove further, we discuss in this Appendix the perturbative solution of the FPequation (2.10) in the stationary limit. Although it is only of limited use, itprovides some insight, especially along the x axis. A.1 Lowest-order solution
We write the FP operator (2.11) as L T = L T + λL T , (A.1)with L T = N R ∂ x + ( Ax − By ) ∂ x + N I ∂ y + ( Ay + Bx ) ∂ y + 2 A, (A.2)and L T = (cid:0) x − xy (cid:1) ∂ x + (cid:0) x y − y (cid:1) ∂ y + 6 (cid:0) x − y (cid:1) . (A.3)The (normalisable) solution of the lowest-order equation, L T P (0) = 0 , (A.4)is given by P (0) ( x, y ) = N exp (cid:2) − αx − βy − γxy (cid:3) , (A.5)with α = AD (cid:2) ( N R + N I )( A + B ) − A (cid:3) , (A.6) β = AD (cid:2) ( N R + N I )( A + B ) + A (cid:3) , (A.7) γ = A BD , (A.8)where D = ( N R + N I ) ( A + B ) − A , (A.9)26nd N is the normalisation constant,1 N = Z dxdy e − αx − βy − γxy = π p αβ − γ . (A.10)This solution is similar to the one found in the relativistic Bose gas at nonzerochemical potential [20]. It is easy to see that it is the correct solution at leadingorder, by computing (recall that N R − N I = 1 and σ = A + iB ) (cid:10) ( x + iy ) (cid:11) P = Z dxdy P (0) ( x, y )( x + iy ) = 1 σ . (A.11)More generally, one may equate the two expectation values h O ( x ) i ρ = Z dx ρ ( x ) O ( x ) , (A.12) h O ( x + iy ) i P = Z dxdy P ( x, y ) O ( x + iy ) , (A.13)which, assuming that it is possible to shift x → x − iy , yields the relation [29, 30] ρ ( x ) = Z dy P ( x − iy, y ) , (A.14)where the LHS should be independent of N R , I . Evaluating the y integral yieldsin this case ρ (0) ( x ) = N ′ e − S ( x ) , S ( x ) = 12 σx , (A.15)with N ′ = r σ π , (A.16)which is indeed the expected answer. A.2 First-order correction
To compute higher-order corrections, we expand P ( x, y ) = ∞ X k =0 λ k P ( k ) ( x, y ) . (A.17)Higher-order corrections are determined by the inhomogeneous partial differentialequation, L T P ( k ) + L T P ( k − = 0 . (A.18)The homogeneous equation is solved by P (0) . To find the particular solution, wefactor out the leading order solution, P ( k ) = P (0) p ( k ) , (A.19)27with p (0) = 1), and write L T P ( k ) = P (0) L ′ T p ( k ) , L T P ( k ) = P (0) L ′ T p ( k ) , (A.20)with L ′ T = N R [ ∂ x − αx + γy )] ∂ x + ( Ax − By ) ∂ x + N I [ ∂ y − βy + γx )] ∂ y + ( Ay + Bx ) ∂ y , (A.21) L ′ T = (cid:0) x − xy (cid:1) ( − αx − γy + ∂ x )+ (cid:0) yx − y (cid:1) ( − βy − γx + ∂ y ) + 6 (cid:0) x − y (cid:1) . (A.22)Higher-order corrections are then determined by L ′ T p ( k ) = − L ′ T p ( k − . (A.23)For the first-order correction, this yields L ′ T p (1) = 2 αx + 8 γx y − α − β ) x y − γxy − βy − (cid:0) x − y (cid:1) . (A.24)The RHS of Eq. (A.24) is a fourth order polynomial with only even powers. As aparticular solution we may therefore attempt a polynomial of fourth degree, withonly even terms appearing and containing 8 unknown coefficients, p (1) ( x, y ) = c x + c x y + c x y + c xy + c y + c x + c xy + c y . (A.25)Inserting this Ansatz in Eq. (A.24) yields a set of linear equations for the coef-ficients which can be solved. Since the expressions become rather unwieldy, wegive here the results for real noise only, since this is the case of interest.For real noise ( N R = 1 , N I = 0), the parameters in the lowest-order solution(A.5) simplify, and α = A, β = A (cid:18) A B (cid:19) , γ = A B . (A.26)The coefficients of the first-order correction (A.25) are given by c = 0 , c = 12 A (2 A − B ) B (4 A + B ) , (A.27) c = − A + B ) B (4 A + B ) , c = − A (4 A − B ) B (4 A + B ) , (A.28) c = − A A + B ) , c = − A (36 A − B )2 B , (A.29) c = − A (5 A − B ) B (4 A + B ) , c = − A (36 A − A B − B ) B (4 A + B ) , (A.30)28igure 15: Distribution P ( x, y ) in the xy -plane, at first nontrivial order in aperturbative expansion, at σ = 1 + i and λ = 1, for real noise.Hence, to first order, the (normalised) distribution is given by P ( x, y ) = N P (0) ( x, y ) (cid:2) λp (1) ( x, y ) (cid:3) , (A.31)with 1 N = 1 − A + 3 A B + 2 B )2( A + B ) (4 A + B ) λ. (A.32)This distribution satisfies the FP equation to order O ( λ ). It can be checked thatit yields the correct moments to this order, e.g. (cid:10) ( x + iy ) (cid:11) P = Z dxdy P ( x, y )( x + iy ) = 1 σ − λσ + O ( λ ) . (A.33)One may also verify that evaluating ρ ( x ) = Z dy P ( x − iy, y ) (A.34)yields in this case ρ ( x ) = N ′ e − σx (cid:18) − λ x (cid:19) + O ( λ ) , (A.35)with N ′ = r σ π (cid:18) − λ σ (cid:19) , (A.36)as it should be. 29 x P ( x , ) FPEleading order perturbativefirst order perturbative σ = 1+ i, λ = 1, N I = 0 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 y P ( , y ) FPEleading order perturbativefirst order perturbative σ = 1+ i, λ = 1, N I = 0 Figure 16: Comparison between the perturbative distribution and the solutionof the FPE, for P ( x,
0) (left) and P (0 , y ) (right), at σ = 1 + i and λ = 1, for realnoise. For the solution of the FPE, ω = 50 and N H = 150.It is clear that the perturbative distribution is not positive definite and strictlyspeaking only applies when the perturbative correction λp (1) ( x, y ) is small withrespect to 1, i.e. around the origin. However, it can be made positive definite bya simple exponentiation, P ( x, y ) = P (0) ( x, y ) exp (cid:2) λp (1) ( x, y ) (cid:3) , (A.37)which has the same leading order λ dependence. This distribution is normalisablesince the coefficients of the quartic terms are all negative. An example is shownin Fig. 15. We observe a double peak structure, as in the main text.At large y values, the exponentiated construction cannot be correct, since itdecays exponentially rather than be 0 outside the strip found above. In the x direction, however, the perturbative solution gives a surprisingly good descriptionof the decay. Taking y = 0, we find P ( x, ∼ exp( − Ax + c λx ) , (A.38)where, for A = B = 1, c = − /
10. This result is compared with the solutionof the FP equation in Fig. 16 (left), and is seen to agree better than expected.We note that the prefactor 0.3 is also close to what was observed for the inte-grated distribution P x ( x ). In Fig. 16 (right), we also show a comparison with theperturbative expression P (0 , y ) ∼ exp[ − ( β − c λ ) y + c λy ] , (A.39)where β = 3, c = 12 / c = − / A = B = 1). Even though c is positive, it is not large enough to change the curvature. Note that theoscillations visible in the solution of the FPE are due to the finite number ofbasis functions ( N H = 150). 30 eferences [1] G. Parisi, Phys. Lett. B (1983) 393.[2] J. R. Klauder, Stochastic quantization, in: H. Mitter, C.B. Lang (Eds.),Recent Developments in High-Energy Physics, Springer-Verlag, Wien, 1983, p.351; J. Phys. A: Math. Gen. (1983) L317; Phys. Rev. A (1984) 2036.[3] G. Aarts and I.-O. Stamatescu, JHEP (2008) 018 [arXiv:0807.1597 [hep-lat]].[4] G. Aarts, Phys. Rev. Lett. (2009) 131601 [arXiv:0810.2089 [hep-lat]].[5] G. Aarts and K. Splittorff, JHEP (2010) 017 [arXiv:1006.0332 [hep-lat]].[6] G. Aarts and F. A. James, JHEP (2012) 118 [arXiv:1112.4655 [hep-lat]].[7] M. Fromm, J. Langelage, S. Lottini, M. Neuman and O. Philipsen,arXiv:1207.3005 [hep-lat].[8] E. Seiler, D. Sexty and I. -O. Stamatescu, Phys. Lett. B (2013) 213[arXiv:1211.3709 [hep-lat]].[9] G. Aarts, PoS LATTICE (2012) 017 [arXiv:1302.3028 [hep-lat]].[10] G. Aarts, L. Bongiovanni, E. Seiler, D. Sexty and I. -O. Stamatescu, Eur.Phys. J. A (to appear) [arXiv:1303.6425 [hep-lat][.[11] J. Ambjorn and S. K. Yang, Phys. Lett. B (1985) 140.[12] J. Ambjorn, M. Flensburg and C. Peterson, Nucl. Phys. B (1986) 375.[13] J. Berges, S. Borsanyi, D. Sexty and I. O. Stamatescu, Phys. Rev. D (2007) 045007 [hep-lat/0609058].[14] J. Berges and D. Sexty, Nucl. Phys. B (2008) 306 [arXiv:0708.0779[hep-lat]].[15] G. Aarts and F. A. James, JHEP (2010) 020 [arXiv:1005.3468 [hep-lat]].[16] J. M. Pawlowski and C. Zielinski, arXiv:1302.1622 [hep-lat]; arXiv:1302.2249[hep-lat].[17] G. Aarts, E. Seiler and I. -O. Stamatescu, Phys. Rev. D (2010) 054508[arXiv:0912.3360 [hep-lat]].[18] G. Aarts, F. A. James, E. Seiler and I. -O. Stamatescu, Eur. Phys. J. C (2011) 1756 [arXiv:1101.3270 [hep-lat]].3119] P. H. Damgaard and H. H¨uffel, Phys. Rept. (1987) 227.[20] G. Aarts, JHEP (2009) 052 [arXiv:0902.4686 [hep-lat]].[21] A. Duncan and M. Niedermaier, Annals Phys. (2013) 93[arXiv:1205.0307 [quant-ph]].[22] J. R. Klauder and W. P. Petersen, Journal of Stat. Phys. 39 (1985) 53,Print-85-0295 (BTL).[23] H. Okamoto, K. Okano, L. Schulke and S. Tanaka, Nucl. Phys. B (1989)684.[24] G. Aarts, F. A. James, E. Seiler and I. -O. Stamatescu, Phys. Lett. B (2010) 154 [arXiv:0912.0617 [hep-lat]].[25] C. Pehlevan and G. Guralnik, Nucl. Phys. B (2009) 519 [arXiv:0710.3756[hep-th]].[26] G. Guralnik and C. Pehlevan, Nucl. Phys. B (2009) 349 [arXiv:0902.1503[hep-lat]].[27] E. Anderson et al, LAPACK Users’ Guide (Third Ed.), Society for Industrialand Applied Mathematics, 1999.[28] G. Aarts, F. A. James, J. M. Pawlowski, E. Seiler, D. Sexty and I. -O. Sta-matescu, JHEP (2013) 073 [arXiv:1212.5231 [hep-lat]].[29] H. Nakazato, Prog. Theor. Phys. (1987) 20.[30] M. Namiki, I. Ohba, K. Okano, Y. Yamanaka, A. K. Kapoor, H. Nakazatoand S. Tanaka, Lect. Notes Phys. M9