A study of symmetry breaking in a relativistic Bose gas using the contraction algorithm
Andrei Alexandru, Gokce Basar, Paulo Bedaque, Gregory W. Ridgway, Neill C. Warrington
AA study of symmetry breaking in a relativistic Bose gas using the contractionalgorithm
Andrei Alexandru
Department of Physics, George Washington University, Washington, DC 20052 andDepartment of Physics, University of Maryland, College Park, MD 20742
G¨ok¸ce Ba¸sar, Paulo Bedaque, Gregory W. Ridgway, and Neill C. Warrington
Department of Physics, University of Maryland, College Park, MD 20742 (Dated: March 12, 2018)A relativistic Bose gas at finite density suffers from a sign problem that makes direct numericalsimulations not feasible. One possible solution to the sign problem is to re-express the path integralin terms of Lefschetz thimbles. Using this approach we study the relativistic Bose gas both inthe symmetric phase (low-density) and the spontaneously broken phase (high-density). In thehigh-density phase we break explicitly the symmetry and determine the dependence of the orderparameter on the breaking. We study the relative contributions of the dominant and sub-dominantthimbles in this phase. We find that the sub-dominant thimble only contributes substantially whenthe explicit symmetry breaking is small, a regime that is dominated by finite volume effects. In theregime relevant for the thermodynamic limit, this contribution is negligible.
I. INTRODUCTION
The Monte Carlo method is a powerful tool to studyfield theories even when perturbation theory and otheranalytical methods fail. The method is usually appliedin the Euclidean time formalism where the path inte-gral defining the field theory becomes equivalent to astatistical mechanics problem: each field configuration issampled according to a Boltzmann factor exp( − S ), with S being the Euclidean action for the configuration, andthe observables are expressed as correlations of compositefield operators. The Monte Carlo method, however, hasits limitations. When the action is not real, as is oftenthe case in the presence of a chemical potential, this ap-proach fails. A possible solution is to use the real partof the action S R for importance sampling and combinethe imaginary part with the observable, that is, to re-place O with O exp( − iS I ). For theories that exhibit the sign problem , the phase factor fluctuates substantiallyand the signal-to-noise ratio for (cid:104) exp( − iS I ) (cid:105) decreasesexponentially fast as the volume increases. This rendersthe stochastic approach unfeasible.Recently it was suggested that the sign problem can besolved by re-expressing the path integral using Lefschetzthimbles [1]. The idea is to take the original path-integralexpressed in terms of n real degrees of freedom and viewit as an integral over the R n submanifold embedded in C n . Smooth deformations of the integration domain willnot change the value of the integral as long as we do notcross any singularities of the integrand and we keep theintegral finite. While the integral remains the same, thefluctuations of the integrand on some of these manifoldscould be reduced, making this formulation better suitedfor numerical sampling. A particular choice, and in acertain sense the optimal one, is the decomposition of theoriginal integration domain in terms of Lefschetz thimbles,manifolds that have constant S I [1]. If the decomposition involves only one thimble, or if one thimble dominates thepath integral, then the sign problem is solved. This doesnot seem to be the case for theories in 0 + 1 dimensions,even in the continuum limit [2–4], but there is still a hopethat this might happen in the thermodynamic limit or inthe continuum limit for quantum field theories [1]. Wenote that even if this is not the case, there is a possibilitythat alternative methods based on other manifolds mightbe numerically useful [5].In this paper we study the relativistic Bose gas atfinite density using the contraction algorithm. This sys-tem has a complex action leading to a sign problem. Inthe last few years this system was used as a test casefor different approaches to the sign problem: complexLangevin [6], dual variables [7, 8], mean field [9, 10], den-sity of states [11], and Lefschetz thimbles [12, 13]. Ourgoal is to understand the interplay between spontaneoussymmetry breaking and the contribution of non-dominantthimbles, to showcase the contraction algorithm [2] andsome of the optimizations we developed [14], and to testsome ideas about using alternative manifolds as integra-tion domains [5]. In a recent study it was conjectured thatthe tangent space at the critical point of the dominantthimble is equivalent with the integration over the originaldomain [12]. We confirm this conjecture numerically byshowing that the results of our algorithm match very wellthe results from alternative approaches. Using this resultas a benchmark we check the relative contribution of thedominant thimble.The plan of the paper is the following. In Section IIwe review the relevant details about the thimbles anddescribe the algorithm we used to perform the integralover the thimble(s). In Section III we introduce the actionfor the relativistic Bose gas, its discretization and discussthe details relevant for our implementation. In Section IVwe present the results of our method and compare it withresults from other approaches. In Section V we drawconclusions and discuss future directions. a r X i v : . [ h e p - l a t ] J un II. LEFSCHETZ THIMBLES ANDTHE CONTRACTION ALGORITHM
Thermal expectation values of observables in a bosonicfield theory have the path integral expression (cid:104)O(cid:105) = 1 Z (cid:90) Dφ e − S ( φ ) O ( φ ) with Z = (cid:90) Dφ e − S ( φ ) , (2.1)where the integration is taken over real fields φ , and S is the Euclidean action. In practice, to compute thepath integral above, one first regulates the theory byapproximating spacetime by a lattice of points. The effectof this regularization on the path integral is to change theintegration domain from the space of field configurationsto R n , where n is the number of degrees of freedom. When S is real, the observables can be evaluated by a MonteCarlo integration, that is, (cid:104)O(cid:105) is estimated as (cid:104)O(cid:105) ≈ N N (cid:88) a =1 O ( φ a ) , (2.2)where the field configurations φ a are distributed accordingto the probability distribution Pr( φ ) = exp[ − S ( φ )] /Z .The Monte Carlo method does not work when the action S = S R + iS I is complex. A possible solution to thisproblem is the so-called “reweighting” procedure whichamounts re-expressing the original observable as a ratioof observables, which is amenable to a standard MonteCarlo evaluation: (cid:104)O(cid:105) = 1 Z (cid:90) Dφ e − ( S R + iS I ) O = 1 Z/Z R Z R (cid:90) Dφ e − S R (cid:0) e − iS I O (cid:1) = (cid:10) e − iS I O (cid:11) S R (cid:104) e − iS I (cid:105) S R (2.3)where (cid:104)·(cid:105) S R denotes an average with respect to theprobability distribution exp( − S R ) /Z R with Z R = (cid:82) Dφ exp[ − S R ( φ )]. Unfortunately the denominator inthe equation above scales as exp( − βV ∆ f ), where β isthe inverse temperature, V is spatial volume, and ∆ f isthe free energy difference between the system describedby actions S R and S . Consequently, the denominatorgoes to zero exponentially fast as the lattice volumes in-creases, and any reweighted Monte Carlo is dominated bystatistical errors.An elegant, geometric solution to the sign problem wasrecently proposed in [1]. The strategy is to first complexifythe field variables φ , then deform the integration domainfrom R n to another submanifold of C n . A particularlyjudicious choice of submanifold is a set of “Lefschetzthimbles”, which we will denote by {T σ | σ = 1 , ... } . Wefirst define thimbles, then explain their utility in solvingthe sign problem. For every critical point φ c of the action,defined by ∂S∂φ (cid:12)(cid:12)(cid:12)(cid:12) φ c = 0 , (2.4) there is an associated thimble. The thimble T is definedas the set of all points that, when evolved using the downward flow equations dφ i dt = − ∂S∂φ i , (2.5)converge to φ c as t → ∞ . An analogous upward flow isdefined by flipping the sign on the RHS of Eq. 2.5 and anunstable thimble K can be defined with respect to this flow.Denoting the complexified field variables φ i = φ R,i + iφ I,i ,decomposing the effect of the downward flow on its realand imaginary parts and utilizing the Cauchy-Riemannconditions, it can be seen that any trajectory along thedownward flow follows a negative gradient flow of withrespect to S R and a Hamiltonian flow with respect to S I : dφ R,i dt = − ∂S R ∂φ R,i = ∂S I ∂φ I,i ,dφ
I,i dt = − ∂S R ∂φ I,i = − ∂S I ∂φ R,i . (2.6)The decomposition above shows that thimbles arethe multi-dimensional generalization of steepest de-scent/stationary phase paths from asymptotic analysis.The imaginary part of the action is then constant over athimble and it is thus advantageous to deform the regionof integration from R n to thimbles if there is a rapidlyoscillating phase in the integral on R n . This proceduresolves the sign problem. Typically, there are many criticalpoints, and hence many thimbles. Moreover it is a non-trivial task to determine what combination of thimblesis equivalent to the original domain of integration R n .However, a fundamental result of Picard-Lefschetz theoryshows that any integral over R n can be decomposed intointegrals over thimbles [15, 16]: (cid:90) R n dφ e − S ( φ ) O ( φ ) = (cid:88) σ n σ e − iS I ( φ σ ) (cid:90) T σ dφ e − S R ( φ ) O ( φ ) , (2.7)where the summation runs over the critical points φ σ and T σ are the associated thimbles. The integers n σ countthe intersection points between the original integrationcontour R n and the unstable thimble K σ . The n σ can benegative: if the flow takes a volume element around oneof these intersection points in R n into a thimble volumeelement with orientation opposite to that of the thimble,then the intersection point counts negatively toward n σ .Therefore, n σ is the number of points that flow from R n tothe critical point φ σ via the upward flow while preservingorientation, minus the number of intersection points thatreverse orientation.It is in general very difficult to find all critical pointsand the values of the coefficients n σ in the thimble decom-position above. However, assuming none of the n σ arezero, we see from Eq. 2.7 that we can estimate the rela-tive importance of each thimble based on the value of thereal part of the action at the critical point, S R ( φ σ ). Thecritical point with the lowest action, φ ¯ σ , is expected togive the largest contribution to the path integral while thesubdominant ones are suppressed by a factor exp( − ∆ S )with ∆ S = S R ( φ σ ) − S R ( φ ¯ σ ). This estimate, of course,neglects entropy effects and it is valid only to the extentthat the semiclassical expansion is qualitatively correct.When the main thimble dominates the integral averagescan be approximated by (cid:104)O(cid:105) ≈ (cid:90) T ¯ σ dφ e − S R ( φ ) O ( φ ) (cid:46) (cid:90) T ¯ σ dφ e − S R ( φ ) . (2.8)Notice that the phase fluctuations are almost absent inthis case since the imaginary part of the action factors outfor observable averages and the sampling is done usingthe real part of the action. The only remaining phase isthe residual phase , that is the phase associated with thevolume element dφ/ | dφ | , which in general varies smoothlyover the thimbles and can be easily reweighted [2, 13, 17].One of the goals of this paper is to apply the “con-traction algorithm” [2, 5] for a relativistic bose gas atfinite density. We briefly review the algorithm here. Abasic ingredient is the map generated by the upward flow φ n → φ f ( φ n ) dφdt = + ∂S∂φ with φ (0) = φ n and φ f = φ ( T flow ) . (2.9)This is a map from C n to C n . Under this map the valueof S R at every point increases while S I remains fixed.For any integration manifold M n with a finite integral (cid:82) M n dφ exp[ − S ( φ )], it can be shown that the value of theintegral is unchanged if we replace the manifold M n withits image through this map, M f = φ f ( M n ) [5]. Take forthe moment M n to be the original integration contour R n . We have (cid:90) R n dφ e − S ( φ ) = (cid:90) M f dφ f e − S ( φ f ) = (cid:90) R n dφ n det Je − S ( φ f ( φ n )) with J ij = ∂ ( φ f ) i ∂ ( φ n ) j . (2.10)The last step above is derived using a change of variablesfrom φ f to φ n with J being the Jacobian of the map.The contraction algorithm samples the integral usingthe integrand on the RHS of the expression above. Theintegrand is not real, so the sampling is done based onthe Boltzmann factor exp[ − S eff ( φ n )] with S eff ( φ n ) ≡ Re[ S ( φ f ( φ n )) − log det J ]= S R ( φ f ( φ n )) − log | det J | . (2.11)The phase ϕ ( φ n ) ≡ S I ( φ f ( φ n )) − Im log det J isreweighted as in Eq. 2.3. When T flow = 0 the integrationmanifold is unchanged and the phase of the integrandfluctuates rapidly since the original formulation has asign problem. As T flow increases, S R increases and thesampled points concentrate on ever smaller regions of theflowed manifold where S R is small. As these regions are small and S I is preserved by the flow, the phase fluc-tuations on these small sampled regions are also small,alleviating the sign problem. In the T flow → ∞ limit,these sampled regions reduce to isolated points, one foreach contributing thimble. On the flowed manifold, forlarge T flow , the regions with an important statisticalweight form isolated pockets with particularly small S R ,each one corresponding to a particular contributing thim-ble. In this formulation, the full value of the integral isrecovered only when all the pockets are sampled, corre-sponding to the inclusion of all thimbles appearing in thedecomposition. Algorithms involving incremental changesin the field configurations will only sample the pocketclosest to the starting configuration. This is actually howwe employ our algorithm when we are interested in sam-pling only the contribution of one thimble. If we wantto include the contribution of other thimbles we have toalso make proposals that can take us from one pocket toanother. We used this type of large updates previously ina fermionic model [5], and in Section IV we will presenta similar procedure for the model we study in this paper.The contraction algorithm is simply a Metropolis algo-rithm using the variables φ n based on the effective action S eff [ φ f ( φ n ))]. The proposals have to be chosen carefullyto get a reasonable acceptance rate. A number of opti-mizations are possible when the starting manifold M n is not the original integration domain but the tangentmanifold to a thimble of interest. This manifold is alegitimate choice in two cases: either we are interestedin sampling only one thimble in the limit T flow → ∞ [2],or the tangent space to the thimble is equivalent to theoriginal integration domain, which is the case for a classof systems. It was conjectured that this is also true for themodel studied in this paper [12]. We will show evidencesupporting this conjecture.There are two advantages to using the tangent space ofa critical point as the starting manifold instead of R n : wecan make efficient Metropolis proposals [2] and we can usean estimator to take into account the contribution of theJacobian to the effective action [14]. The Jacobian can becomputed exactly by integrating the following differentialequation along the upward flow path, dJdt = HJ , (2.12)where J is the Jacobian matrix and H is the Hessian ofthe action S . The initial condition for J is a matrix thathas as columns a set of vectors that span the tangentspace at the critical point. We use J (0) = (ˆ ρ , . . . , ˆ ρ n ),where ˆ ρ i are the positive “eigenvectors” of the Hessianat the critical point, that is H ˆ ρ i = λ i ˆ ρ i with λ i >
0. Itcan be shown that these vectors span the tangent spaceto the thimble T and that i ˆ ρ i are negative “eigenvectors”that span the tangent space to the unstable thimble K .Integrating the equation above is numerically expensive,since every step of the integrator involves multiplicationof H and J matrices. Fortunately, there is a very goodestimator for log | det J | built out of the “eigenvectors”ˆ ρ i [14]: W = (cid:90) T flow dt n (cid:88) i =1 ˆ ρ † i H ( t )ˆ ρ i . (2.13)Since H is sparse, the cost of computing W is O ( n ),which is much cheaper than the O ( n ) cost of computing J ; the savings are substantial even for small lattices.The algorithm samples configurations with an effectiveaction S (cid:48) eff ( φ n ) = S R ( φ f ( φ n )) − Re W and the difference∆ S = S eff ( φ n ) − S (cid:48) eff ( φ n ) is included as a reweightingfactor. To compute the reweighting factor exactly weneed to integrate the equations for J , but this needs tobe done only for the small subset of configurations usedfor measuring observables, which are separated by a largenumber of Metropolis steps. III. RELATIVISTIC BOSE GAS AT FINITEDENSITY
The Euclidean action of a gas of relativistic spinlessbosons in three (spatial) dimensions at finite density isgiven by S = (cid:90) d x (cid:2) ∂ φ ∗ ∂ φ + ∇ φ ∗ · ∇ φ + ( m − µ ) | φ | + µ ( φ ∗ ∂ φ − φ∂ φ ∗ ) (cid:124) (cid:123)(cid:122) (cid:125) j ( x ) + λ | φ | (cid:3) (3.1)where φ ≡ ( φ + iφ ) / √ j ( x ) is imaginary andis the source of the sign problem in this model. Thisaction is symmetric under global U (1) transformations φ → e iα φ . For values of µ below a critical value of order m the equilibrium state is expected to be U (1) symmetricand (cid:104) φ (cid:105) = 0. For values of µ larger than the critical valueand for low enough temperatures the U (1) symmetry isexpected to be spontaneously broken and (cid:104) φ (cid:105) (cid:54) = 0.In order to study spontaneous symmetry breaking itis necessary to introduce an explicit symmetry breakingterm S → S − h (cid:82) d x ( φ + φ ). This choice breaks theoriginal U (1) symmetry down to the Z sub-group givenby swapping φ x, ↔ φ x, .We use the following discretization of the action: S = a (cid:88) x (cid:34) e µa φ ∗ x +ˆ0 − φ ∗ x a e − µa φ x +ˆ0 − φ x a + (cid:88) ν =1 (cid:12)(cid:12)(cid:12) φ x +ˆ ν − φ x a (cid:12)(cid:12)(cid:12) + m | φ x | + λ | φ x | − h ( φ x, + φ x, ) (cid:35) . (3.2)Setting the lattice spacing to a = 1 and writing the field in terms of real components, φ and φ , we obtain S = (cid:88) x (cid:34)(cid:16) m (cid:17) φ x,a φ x,a − (cid:88) ν =1 φ x,a φ x +ˆ ν,a − cosh µ φ x,a φ x +ˆ0 ,a + i sinh µ (cid:15) ab φ x,a φ x +ˆ0 ,b + λ (cid:0) φ x,a φ x,a ) − h ( φ x, + φ x, ) (cid:35) (3.3)with (cid:15) ab the antisymmetric tensor with (cid:15) = 1. Thederivatives of the action are ∂S∂φ x,a =(8 + m ) φ x,a − (cid:88) ν =1 ( φ x +ˆ ν,a + φ x − ˆ ν,a ) − cosh µ (cid:16) φ x +ˆ0 ,a + φ x − ˆ0 ,a (cid:17) + i sinh µ(cid:15) ac (cid:16) φ x +ˆ0 ,c − φ x − ˆ0 ,c (cid:17) + λ ( φ x,c φ x,c ) φ x,a − h ,∂ S∂φ x,a ∂φ y,b =(8 + m ) δ xy δ ab − (cid:88) ν =1 ( δ x +ˆ ν,y + δ x − ˆ ν,y ) δ ab − cosh µ (cid:16) δ x +ˆ0 ,y + δ x − ˆ0 ,y (cid:17) δ ab + i sinh µ (cid:16) δ x +ˆ0 ,y − δ x − ˆ0 ,y (cid:17) (cid:15) ab + λδ xy [( φ x,c φ x,c ) δ ab + 2 φ x,a φ x,b ] . (3.4)It is very difficult to find all critical points of the (lat-ticized) action. However, critical points with small S R play a special role as their thimbles are more likely todominate the path integral. Critical points with constantfields φ xa = φ are among the ones with smaller valuesof the action and, at the same time, are easy to find. Infact, the global minimum of the real part of the actionrestricted to real values of φ x, and φ x, is a constantfield (assuming h real). This can be seen by looking atEq. 3.1 or its discretized version Eq. 3.2 and noticing thati) the kinetic and gradient term are positive definite andfavor constant fields and ii) the term linear in µ is purelyimaginary and does not contribute to the real part ofthe action. This motivates us to find the constant fieldcritical points. They are given by φ x, = φ x, = φ where φ is one of the roots of the cubic equation:(2 + m ) φ − µφ + λ | φ | φ = h, (3.5)whose three solutions are given, for small h , by: φ = − h ˆ µ − m + O ( h ) ,φ ± = ± (cid:114) ˆ µ − m λ + h µ − m + O ( h ) , (3.6)with ˆ µ = 2 cosh µ −
2. The corresponding values of the Re Φ Im ΦΜ (cid:60) Μ c h (cid:206) (cid:82) Φ (cid:45) Φ (cid:43) Φ Re Φ Im ΦΜ (cid:62) Μ c h (cid:206) (cid:82) Φ (cid:45) Φ Φ (cid:43) Re Φ Im ΦΜ (cid:62) Μ c h (cid:207) (cid:82) Φ (cid:45) Φ Φ (cid:43) FIG. 1. Projections in the constant φ = φ subspace for T , T + and T − , the thimbles attached to φ , φ + and φ − . The left andmiddle graph correspond to h real and the one on the right to a value of h that has a small imaginary component. Thimbles T i are plotted in blue, unstable thimbles K i are plotted with dashed, red lines, and critical points are indicated by black dots. action at these points are: S ( φ ) V = 0 + O ( h ) ,S ( φ ± ) V = − λ (cid:18) ˆ µ − m (cid:19) ∓ h (cid:114) ˆ µ − m λ + O ( h ) , (3.7)where V = ( V /a ) / ( aT ) is the number of lattice sites.We now turn to discuss what we know about the thimbledecomposition of the original integral over real fields. For µ below µ c = cosh − (1 + m /
2) the critical points φ ± areimaginary. In this case the thimbles T ± associated with φ ± do not contribute to the thimble decomposition. Thisis because the real part of the action in R n is boundedfrom below by Re S ( φ ) which is larger than Re S ( φ ± ).Thus no point in the original integration domain canflow into φ ± under steepest-ascent/upward-flow definedin Eq. 2.9. This is illustrated in the the left panel ofFig. 1 which shows a projection of the thimbles T , ± inthe φ = φ plane. We can see that, in the µ < µ c case,the unstable thimbles K ± do not intersect the originalintegration domain, which is the real axis in this figure.For values of µ greater than µ c the situation changes. Inthe case where µ > µ c and h is real, there is a trajectoryof the flow connecting φ and φ + ( φ − ). Indeed, in theconstant field subspace (that is, in the subspace where φ x, = φ x, = φ ) the gradient ∂S∂φ x,a (cid:12)(cid:12)(cid:12)(cid:12) φ x, = φ x, = φ = ( m − ˆ µ ) φ + 2 λφ − h (3.8)points along the constant field subspace (all componentsare equal). Moreover, the gradient is real if φ is real. Thisimplies in the existence of flow trajectories lying entirelyon the real constant field subspace. Since the downwardflow decreases the value of the real part of the action, we conclude that there is a trajectory connecting the real con-stant fields φ and φ + (and another connecting φ to φ − ).Thus, the unstable thimble of φ ( K ) overlaps the thim-ble of φ + ( T + ). The existence of these lines is known as aStokes phenomenon and invalidates the decomposition ofthe integral into integer linear combinations of thimbles.We bypass this problem by making h slightly complex.As can be seen from Eq. 3.7 a complex value of h impliesdifferent values for the imaginary parts of S ( φ + ) , S ( φ − )and S ( φ ). As the flows preserves the imaginary part ofthe action, a complex value of h guarantees that there isno Stokes lines connecting these critical points. The waya complex value of h makes the thimble decompositionwell defined is shown visually in the center and right pan-els of Fig. 1. The figure suggests that the integral overthe real line is equivalent to the integral over T + and T − (with the proper orientations), which would imply that n + = n − = 1 , n = 0. The Fig. 1 only shows a projectionof the whole space but some definite conclusions can bedrawn. For instance, the fact that there is a flow lineconnecting the real plane to φ − shows that K − does in-tersect R n , although the possibility remains that thereare other (non-constant field) points where K − intersects R n . This is a strong but not definitive case that n − = 1.The fact that n + = 1 can be argued even more rigorously.The difference in the real part of the action between φ + and the global minimum on the real plane is proportionalto Im( h ). In the Im( h ) → n + = 1.To the extent that the constant field projection can betrusted we also have n = 0. Unfortunately, the pictureaway from the constant field subspace is much harderto analyze and there is the possibility that T , or eventhat thimbles corresponding to other, non-constant fieldcritical points, may also contribute. Even if the unstablethimble K were to intersect the real integration domainat a point not included in the projection in Fig. 1, andtherefore contribute to the thimble decomposition of theoriginal integral, its contribution is expected to be muchsmaller than the one from T − . In fact, the leading orderestimate (in the semiclassical expansion) for the relativecontribution of two thimbles, say T and T + , is givenby the ratio of Boltzmann factors exp( − [ S ( φ ) − S ( φ + )]).Since the difference in the action at φ and φ + is ap-proximately (ˆ µ − m ) / λ , the contribution of T , if notidentically zero, is exponentially small as V (cid:29) µ (cid:29) m or λ (cid:28) dWdt = Tr (cid:48) H with Tr (cid:48) H ≡ (cid:88) i ˆ ρ † i H ( φ ( t ))ˆ ρ i (3.9)together with the upward flow in Eq. 2.9. These differen-tial equations are integrated using a Cash-Karp integrator,an adaptive step-size integrator of order O (∆ t ) [19]. Forevery integrator step we need to evaluate Tr (cid:48) H ( φ ( t )),which requires the positive “eigenvectors” ˆ ρ i for the Hes-sian at the critical point φ cr H ( φ cr )ˆ ρ i = λ i ˆ ρ i , with λ i > , (3.10)where H ( φ ) x,a ; y,b = ∂ Sφ x,a φ y,b , and ˆ ρ † i ˆ ρ j = δ ij . (3.11)We stress that the estimator involves only the “eigen-vectors” of the Hessian at the critical point and not theones of H ( φ ). This estimator is exact when the action isquadratic and it tracks well the true value of the Jacobianeven when the thimble is curved [14]. We compute the“eigenvectors” before starting the Monte-Carlo simulation.This step is computationally costly but needs to be per-formed only once. To evaluate the estimator, we separatethe φ -dependent part from the Hessian as: H ( φ ) x,a ; y,b = H ( φ cr ) x,a ; y,b + λδ xy (cid:2) ( φ x,c φ x,c ) δ ab + 2 φ x,a φ x,b − ( φ → φ cr ) (cid:3) . (3.12)We have thenTr (cid:48) H ( φ ) = Tr (cid:48) H ( φ cr ) + ∆( φ ) − ∆( φ cr )= (cid:32)(cid:88) i λ i − ∆( φ cr ) (cid:33) + ∆( φ ) , (3.13)where we used the fact that Tr (cid:48) H ( φ cr ) = (cid:80) i λ i and thedefinitions∆( φ ) ≡ (cid:88) x (cid:88) abc R x,ab [ φ x,c φ x,c δ ab + 2 φ x,a φ x,b ] , R x,ab ≡ (cid:88) i (ˆ ρ i ) x,a (ˆ ρ i ) x,b . (3.14)Note that Tr (cid:48) H ( φ cr ), ∆( φ cr ), and R are computed onceat the beginning of the simulation and we only need toevaluate ∆( φ ) along the integration path. This step hascomputational cost of order O ( n ). IV. RESULTS
In this section we present the results of our simulations.We discuss first the exact results, based on simulations onthe thimble tangent plane at the φ + critical point, bothin the symmetric and the broken phase. We then focuson the broken phase and analyze the dependence of theorder parameter as a function of the breaking parameter h to determine the regime where finite volume effects arenot dominant. After that we report the results of singlethimble simulations on the dominant thimble, T + , andcompare them with exact results. Finally, we present amethod to carry out simulations on two thimbles andstudy the relative weight of T − to T + as a function of h .As we discussed in Section II, we can use the con-traction algorithm to get exact results by flowing theoriginal domain of integration R n . However, if the start-ing manifold is the tangent space at the critical point ona thimble then shorter flow times are required to reducethe phase fluctuations and we can use better proposals forthe Metropolis algorithm and fast Jacobian estimators.This choice is justified if we can show that this manifoldis equivalent with the original integration domain. Cristo-foretti et al [12] conjectured this to be the case for therelativistic Bose gas: they argue that the original inte-gration domain is equivalent to the tangent space for thethimble corresponding to the global minimum of S R over R n . They use real values of h and the global minimumcorresponding to the constant field critical point φ for µ < µ c and φ + for µ ≥ µ c . Note that the transitionfrom φ to φ + is smooth since they are equal for µ = µ c .The manifolds T and T + are the tangent manifolds tothe thimbles T and T + at the critical points φ and φ + respectively. Note that when using a slightly complex h , as we use, the points φ and φ + have small complexcomponents, proportional to Im h . The tangent spaces T and T + are not parallel to R n , even in the limit h →
0, soan analytical proof of the conjectured equivalence cannotbe easily established.We carried out simulations using T flow = 0 using as thestarting manifold T for µ < µ c and T + for µ ≥ µ c . Tomake sure that we stay on the manifold we determinethe “eigenvectors” ˆ ρ i and parametrize each point on themanifold as φ = φ cr + (cid:80) i c i ˆ ρ i , where c i are real coefficients.We use c i to represent the points in the tangent space, andthe updates are done one coefficient at a time with a stepdrawn from an uniform probability distribution in theinterval [ − θ, θ ] where θ = ∆ / √ λ i . This type of proposalscales the step size in each direction so that the increasein the action due to displacements in each direction arecomparable, at least in the region where the action is wellapproximated by its quadratic approximation around φ cr .We tune ∆ to get a good acceptance rate. In Fig. 2 we Contraction algorithmComplex Langevin Μ R e (cid:88) n (cid:92) Contraction algorithmComplex Langevin Μ R e (cid:88) n (cid:92) FIG. 2. Dependence of the charge Re (cid:104) n (cid:105) on the chemical potential on a 4 lattice with parameters m = λ = 1 . h = 10 − + i − . In the left panel we plot only the data below µ mean c ≈ .
15, the mean field estimate for the transition tothe symmetry broken phase [9]. Complex Langevin data is taken from [6] and the curve in the left panel from a mean fieldcalculation that is known to agree very well with exact results results [9]. show the results for the charge (cid:104) n (cid:105) = 1 V ∂ log Z∂µ = 1 V (cid:88) x ( δ ab sinh µ − i(cid:15) ab cosh µ ) φ x,a φ x +ˆ0 ,b . (4.1)In these simulations we set h = 10 − + i − and ∆ = 3 . µ we evaluated 5 × updates and performedmeasurements on configurations separated by 1000 up-dates. We compare our results with the results fromComplex Langevin simulations [6], which are known toagree with the results obtained from the dual variablesapproach [7]. Our results agree very well which stronglysupports the conjecture that the tangent plane is equiva-lent with the original integration domain. We note thatwe observed neither instabilities nor runaways in our sim-ulations, in contrast with the experiences reported byCristoforetti et al [12], even when we set h to purely realvalues, as they did in their study.The simulations above are feasible even at T flow = 0because the sign fluctuations are significantly smallerwhen integrating over the tangent space rather than theoriginal integration domain. In Fig. 3 we show the averagephase, Re (cid:10) e − iS I (cid:11) , and its standard deviation as a functionof the chemical potential in the symmetric phase, both forsimulations using R n as a starting manifold and T . Weshow in the graph both the standard deviation and thestandard deviation of the mean of Re (cid:10) e − iS I (cid:11) so that thenumber of sampled points in the simulations do not playa role in the comparison. We see that as we approachthe transition point, the phase fluctuates rapidly whenwe sample points on R n and the number of sampledconfigurations required to measure the observables with agiven precision grows very quickly. In contrast, the phasefluctuations on T are very mild, at least for the latticevolume used in these simulations, and the phase can be easily accounted for by reweighting.As discussed earlier, when the symmetry breaking pa-rameter h goes to zero, the order parameter (cid:104) φ (cid:105) = 1 V ∂ log Z∂h , (4.2)vanishes if the limit h → h is small the fluctuations in φ are large enoughto overcome the potential barrier between φ + and φ − ,∆ S = V h (cid:112) (ˆ µ − m ) / (2 λ ) + O ( h ). Since we discusshere simulations at fixed volume, it is important to de-termine the range of h where these finite volume effects (cid:82) n T Μ R e (cid:60) e (cid:45) i S I (cid:62) FIG. 3. Dependence of the average phase on the chemicalpotential µ below the critical transition on the real plane andthe tangent plane T for a 4 lattice with the parameters m = λ = 1 .
0. The orange bars indicate the standard deviationof Re e − iS I . The phase hardly varies for these parameters, soin this case it is sufficient to shift the integration domain to T to tame the sign problem. R e (cid:88) Φ (cid:92) R e (cid:88) n (cid:92) FIG. 4. Observables as a function of the symmetry breaking parameter h on a 4 lattice (blue) and 6 lattice (red) with m = λ = 1 . µ = 1 .
3. Note that the charge varies slowly as we decrease h , whereas for Re h < ∼ .
02 the order parameter (cid:104) φ (cid:105) becomes dominated by finite volume effects and, as expected, approaches zero as h →
0. As expected, the values of h for whichthere is a significant finite volume effect decreases as the volume is increased. are important. In Fig. 4 we plot the charge and the orderparameter (cid:104) φ (cid:105) as a function of h , as determined fromsimulations over T + . We set the value of the chemicalpotential to µ = 1 .
3, to make sure we are in the high-density, symmetry broken phase. In these simulations theratio between the imaginary and real part of h is keptfixed: Re h/ Im h = 10. We see that for values smallerthan Re h < ∼ .
02 (for the 4 lattice) and Re h < ∼ .
01 (forthe 6 lattice) the finite volume effects become importantand restore the (near) symmetry. For the following calcu-lation we use a value of h = 0 . × (1 + i/
10) to make surewe are away from the region where finite volume effectsdominate.Note that for these simulations, at small values of h ,there is a nearly flat direction in field space that is sampledinefficiently using the proposals discussed earlier. Thiscan be easily fixed. The flat direction corresponds roughlyto a circle in the two-dimensional space spanned by ˆ ρ , the“eigenvector” nearly parallel with the constant field with φ = − φ (the “Goldstone” direction) and “eigenvector”ˆ ρ nearly parallel with the constant field direction with φ = φ (the “Higgs” direction). The updates along theother directions are treated as discussed above, but inthe ˆ ρ , plane we use a polar representation and updatethe angular part with steps of size ∆ / √ λ and the radialpart with steps of size ∆ / √ λ . The polar coordinatesare defined in relation to the center c = 0 and c =( φ − − φ + ) · ˆ ρ /
2, with c , the coordinates in the ˆ ρ , basis. To preserve detailed balance we have to modify theaccept/reject step to take into account the asymmetryin this polar proposal, that is we accept the update withprobability exp[ − S R ( φ (cid:48) ) + S R ( φ )] r (cid:48) /r , where r and r (cid:48) represent the radial coordinates for φ and φ (cid:48) in the polarrepresentation.As we discussed in Section III, to obtain exact resultsfor this model using thimble sampling we have to considera set of thimbles, at the very least T + and T − . One of the questions we want to address here is whether the singlethimble calculation is sufficient to recover the exact result.Note that previous calculations for the relativistic bosongas that used Lefschetz thimbles were either carried outin the tangent plane [12], which as we argued provides theexact result, or without including the symmetry breakingterm in which case the critical point is actually a circle andthe integration was done over the entire set of thimblesattached to this circle [13]. As such, this question has notyet been directly addressed.In Section II we showed that to perform a one thimblecalculation with the contraction algorithm, say T σ , weneed to sample the manifold produced by flowing thetangent space of the thimble at the critical point T σ , inthe limit T flow → ∞ . If we start at the critical point φ σ and perform only small step updates, for large flow timesthe algorithm will only sample T σ . We focus here on thesymmetry broken phase and we carry out a simulationfor µ = 1 .
3. The thimble we sample is T + . As the limit T flow → ∞ cannot be taken in practice, it is necessary todevise an operational definition. We carry out simulationsat increasing flow times and monitor the observables ofinterest and the imaginary part of the action on the sam-pled manifold. As the sampled manifold approaches thethimble the observables should converge to their averageover the thimble and the fluctuations of the imaginarypart of the action around S I ( φ + ) should be reduced dras-tically in amplitude (in the limit T flow → ∞ all points inthe sampled pocket around φ + should have the same S I ).We carried out simulations for µ = 1 . h = 0 . × (1 + i/
10) and measured the order parameter and the chargeas a function of T flow . For each simulation we carried out5 × updates and we adjusted the step size to get anacceptance rate close to 50%. Note that the step sizeneeds to be decreased dramatically, since the updatesare carried out in the parametrization manifold, wherethe distribution of configurations become tightly packed T flow R e (cid:88) Φ (cid:92) T flow R e (cid:88) n (cid:92) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) T flow S I T flow (cid:88) e i (cid:106) (cid:92) FIG. 5. T flow extrapolations for the order parameter and the charge (top row) and the average imaginary part of the action S I ( φ f ) and the residual phase exp iϕ = dφ/ | dφ | which is computed as ϕ = Im log det J (bottom row) on a 4 lattice with m = λ = 1 . h = 0 . i/ µ = 1 .
3. For S I we also indicate the standard deviation to gauge its fluctuations. around the critical point. We have to use anisotropicproposals along the “eigenvectors” in the tangent planeto take into account the fact that the flow compressesthese directions differently, which is done by adjustingthe step size in direction ˆ ρ i to ∆ exp( − λ i T flow ) / √ λ i [2].The measurements were carried out on configurationsseparated by 10 ,
000 updates.In Fig. 5 we show our results. We see that the valuesof (cid:104) φ (cid:105) and (cid:104) n (cid:105) show very little dependence on the flowtime, indicating that the single thimble result is equalto the exact result, at least at the level of the errobars.In the bottom row of the figure we also plot the value ofthe imaginary part of the action for the configurationssampled in our Monte Carlo simulations. It is clear thatthe fluctuations die out quite quickly, and as T flow ap-proaches 1 . T + only. In the bottom left panelof Fig. 5 we plot the residual phase, that is the phaserelated to the curvature of the sampled manifold, as afunction of flow time. In our formulation, this is the partof the reweighting phase that is due to the imaginarypart of the Jacobian, that is Im log det J . In the limit T flow → ∞ the fluctuations of S I vanish, but the residualphase continues to fluctuate. These fluctuations are notexpected to lead to a sign problem, since the thimbles areexpected to be relatively smooth at least in the regionwhere the bulk of the contribution to the integral comesfrom. This is indeed what we find: the phase fluctua-tions are very close to one for all intermediate manifolds and asymptote to (cid:104) exp( i Im log det J ) (cid:105) = 0 . et al [13] for µ = 1 . (cid:104) exp( i Im log det J ) (cid:105) = 0 . T + ( h = 0 in that calculation), soit is not clear that a geometric measure like the averageresidual phase should be identical.We conclude that in the regime where the finite volumeeffects are small, the partition function is dominated bythe contribution of one thimble T + . However, it is stillworthwhile investigating the contribution due to otherthimbles for two reasons. First, if we are interested instudying the system when the symmetry breaking termis small and the finite volume effects are important, itis likely that we need to include the contribution fromother thimbles. Secondly, while we showed that the twoobservables we measured seem to be saturated by thecontribution of T + , it is possible that these observablesare special in the sense of being insensitive to removingthe samples due to other thimbles.Semiclassical arguments suggest the contribution ofsubdominant thimbles vanishes. If we consider the ratioof the two contributions to the partition function Z + Z − ≡ (cid:82) T + | dφ | e − S R ( φ ) (cid:82) T − | dφ | e − S R ( φ ) (4.3)0we expect that it goes to infinity very quickly as Re h increases. At leading order in the semiclassical expansionthis ratio is Z + Z − ≈ e − [ S R ( φ + ) − S R ( φ − )] ≈ e hV √ (ˆ µ − m ) / (2 λ ) , (4.4)which justifies our expectation. The next-to-leading orderestimate includes the gaussian fluctuations Z + Z − ≈ e − [ S R ( φ + ) − S R ( φ − )] (cid:118)(cid:117)(cid:117)(cid:116) (cid:102) det H ( φ − ) (cid:102) det H ( φ + ) , (4.5)with H ( φ ) the hessian at φ and (cid:102) det H ( φ ) being the prod-uct of positive “eigenvalues” λ i as defined in Eq. 3.10.The next to leading order ratio increases more slowlythan the leading order estimate and it is possible thatthe semiclassical arguments could break down. In theremainder of this section we will study numerically therelative contribution due to T − and compare it with thesemiclassical approximation.We present here a method to sample two thimbles inthe context of the contraction algorithm. Assuming thatthe original integration domain decomposes into a sumover T + and T − , the average observable is given by (cid:104)O(cid:105) = n + (cid:82) T + dφ e − S ( φ ) O ( φ ) + n − (cid:82) T − dφ e − S ( φ ) O ( φ ) n + (cid:82) T + dφ e − S ( φ ) + n − (cid:82) T − dφ e − S ( φ ) . (4.6)Notice that both thimbles appear in both the numeratorand the denominator, which makes a decomposition of (cid:104)O(cid:105) into a straightforward sum of integrals over thimblesnot possible. The integers n + and n − count the number ofintersection points between the parametrization manifoldand the unstable thimbles K + and K − . Each of them lieswithin a different pocket in the parametrization manifold.In our simulations on the tangent plane T + we have onlyseen evidence for one intersection point with K + at φ + and one intersection point close to φ − ( φ − is not includedin T + but it is nearby), which we assume is associatedwith K − . Here we will discuss a method to sample onlythese two pockets, implicitly assuming that | n ± | = 1; thesign of n ± is automatically taken into account correctlyby the flow.To sample T + and T − simultaneously, we use the con-traction algorithm with large T flow and set T + as theparametrization manifold, which we assume to be equiv-alent with R n . The only difference is that instead ofstarting close to φ + and making only small proposals,we interweave these updates with large proposals thatmove us from the vicinity of φ + to the neighborhood of φ − . Since we parametrize the configurations in T + usingcoefficients in the ˆ ρ i -basis, we implement this proposal byflipping the sign of the coefficient corresponding to theeigenvector ˆ ρ , that is nearly parallel with the constantfield φ = φ configuration. In the limit that T flow → ∞ this process will sample the two thimbles according to the Semiclassical approx LOSemiclassical approx NLOContraction algorithm Z (cid:43) (cid:144) Z (cid:45) FIG. 6. The relative weights of thimble contribution asdetermined from sampling with the contraction algorithm,compared to the predictions from semiclassical approximation. probability density d Pr ± ( φ ) | dφ | = e − S R ( φ ) (cid:46) (cid:90) T + ∪T − | dφ | e − S R ( φ ) . (4.7)As before, the residual phase needs to be folded in theobservable in the reweighting step. Additionally, there isanother fluctuating phase, even in the T flow → ∞ limit,since the imaginary part of the action is different on thetwo thimbles ( S I ( φ + ) (cid:54) = S I ( φ − )). In the regime we studyhere, both these phases lead to mild fluctuation and thereis no problem reweighting them.Note that for computing the reweighting factor we donot need to identify the thimble associated with any ofthe sampled points. The flow automatically produces theright result. In fact, the identification might not evenbe possible for small T flow since the separation into T ± contributions is not sharply defined. However, as weincrease T flow it is easy to identify the associated thimble h R Z /Z LO NLO5 × − . × − . × − . × − . × − . × − . .
02 10 / . × . × . / . × . × TABLE I. Distribution of points among the two lowest actionthimbles T + and T − over the course of a 10 × step MonteCarlo on a 4 lattice with parameters m = λ = 1 . µ = 1 . T flow = 1 .
0. For large values of h not a single transition occursbetween T + and T − over the entire course of a Monte Carlosimulation. The results are compared with the leading order(LO) and next to leading order (NLO) semiclassical predictions. T + since theyconcentrate very sharply around φ ± . We run a set of10 × updates for µ = 1 . m = 1, λ = 1 and a seriesof values for h , while we keep as before Re h/ Im h = 10.We set T flow = 1, so that the associated thimble canbe identified easily and the ratio of sampled points in T + and T − can be computed. The results are plottedin Fig. 6. We see that for small values of h , where T − has a non-negligible contribution, this ratio is actuallyvery close to the predictions of the next-to-leading ordersemiclassical approximation. We conclude that the non-gaussian fluctuations do not affect the sampling ratiosignificantly and that the semiclassical arguments can betrusted. In Table I we record the measured values for thisratio. We also include the results for two simulations withlarge values of h , similar to the ones we discussed earlier inthis section, and we see that the subdominant thimble T − is never visited, in agreement with semiclassical prediction.We conclude that for values of h where the finite volumeeffects are small, the subdominant thimbles are indeedunlikely to contribute significantly. V. CONCLUSIONS
We analyzed the relativistic Bose gas using contractionalgorithm, presenting the first application of the contrac-tion algorithm to a quantum field theory. In studying thehigh density broken phase of the theory, we noticed theexistence of a Stokes line and fixed this problem by usinga complex value of the symmetry breaking parameter h . We verified that the results obtained for the chargedensity agreed with previous calculations when available.We then focused on the order parameter (cid:104) φ (cid:105) , which issensitive to spontaneous symmetry breaking. We firstdetermined the values of h for which finite volume ef-fects, which tend to restore the symmetry, are small. Theresults from complex Langevin calculations agree with our tangent plane calculations, lending support to theconjecture that the tangent plane is indeed equivalent tothe real plane. In contrast to [12] we found no runawaytrajectories, as expected if the tangent plane is in factequivalent to the real plane.We then showed how to use the contraction algorithmto isolate the contribution of a single thimble. The resultsobtained from the single thimble agree with the resultsobtained from the tangent plane within a few percent,indicating that the contributions from other thimbles arenegligible for the parameters explored. We generalizedthe contraction algorithm to perform calculations overtwo thimbles. We found that the contribution from thesub-dominant thimble follows closely the semiclassicalestimates (negligible at large h but comparable to theleading one at smaller h ).Our calculations point towards some obvious extensionsand generalizations. The most pressing one is perhapsthe extension to the thermodynamic limit, with the goalof determining whether the contribution of any otherthimble survives this limit. Also, our group has recentlydeveloped similar technology to that described in thispaper to study the real time dynamics of a simple quantummechanical model [20]. We look forward to combiningthe experience acquired in the present paper to study thereal time dynamics of the Bose gas with an eye towardsthe computation of transport coefficients. ACKNOWLEDGMENTS
A.A. is supported in part by the National Science Foun-dation CAREER grant PHY-1151648. A.A. gratefullyacknowledges the hospitality of the Physics Departmentat the University of Maryland where part of this workwas carried out. G.B., P.B., G.R and N.C.W. are sup-ported by U.S. Department of Energy under Contract No.DE-FG02-93ER-40762. [1] M. Cristoforetti, F. Di Renzo, and L. Scorzato (Aurora-Science), Phys. Rev.
D86 , 074506 (2012), arXiv:1205.3996[hep-lat].[2] A. Alexandru, G. Basar, and P. Bedaque, Phys. Rev.
D93 , 014504 (2016), arXiv:1510.03258 [hep-lat].[3] H. Fujii, S. Kamata, and Y. Kikukawa, JHEP , 078(2015), [Erratum: JHEP02,036(2016)], arXiv:1509.08176[hep-lat].[4] H. Fujii, S. Kamata, and Y. Kikukawa, JHEP , 125(2015), arXiv:1509.09141 [hep-lat].[5] A. Alexandru, G. Basar, P. F. Bedaque, G. W. Ridgway,and N. C. Warrington, (2015), arXiv:1512.08764 [hep-lat].[6] G. Aarts, Phys. Rev. Lett. , 131601 (2009),arXiv:0810.2089 [hep-lat].[7] C. Gattringer and T. Kloiber, Nucl.Phys. B869 , 56 (2013),arXiv:1206.2954 [hep-lat].[8] C. Gattringer and T. Kloiber, Phys. Lett.
B720 , 210(2013), arXiv:1212.3770 [hep-lat]. [9] G. Aarts, JHEP , 052 (2009), arXiv:0902.4686 [hep-lat].[10] O. Akerlund, P. de Forcrand, A. Georges, and P. Werner,Phys. Rev.
D90 , 065008 (2014), arXiv:1405.6613 [hep-lat].[11] L. Bongiovanni, K. Langfeld, B. Lucini, R. Pellegrini, andA. Rago, (2016), arXiv:1601.02929 [hep-lat].[12] M. Cristoforetti, F. Di Renzo, A. Mukherjee,and L. Scorzato, Phys. Rev.
D88 , 051501 (2013),arXiv:1303.7204 [hep-lat].[13] H. Fujii, D. Honda, M. Kato, Y. Kikukawa, S. Komatsu,and T. Sano, JHEP , 147 (2013), arXiv:1309.4371 [hep-lat].[14] A. Alexandru, G. Basar, P. F. Bedaque, G. W. Ridgway,and N. C. Warrington, (2016), arXiv:1604.00956 [hep-lat].[15] V. Vassiliev, Applied Picard-Lefschetz Theory (AmericanMathematical Society, 2002).[16] J. Milnor,
Morse Theory (Princeton University Press,Princeton, NJ, 1963). [17] A. Mukherjee, M. Cristoforetti, and L. Scorzato, Phys.Rev. D88 , 051502 (2013), arXiv:1308.0233 [physics.comp-ph].[18] In some singular cases this region can shrink to a lineor other manifold with dimension larger than zero butsmaller than n . This possibility can be excluded by notic-ing that, for small Im( h ), the global minimum of S R on R n is very close to φ + , close enough that the quadraticapproximation to the action applies. Within the gaussianapproximation the flow equations can be solved analyti- cally and the requirement that the region with S R smallerthan S R ( φ + ) shrinks to a point can be translated on acondition over the “eigenvectors” ˆ ρ i . We have verified thiscondition numerically.[19] J. R. Cash and A. H. Karp, ACM Trans. Math. Softw.16