Complex Langevin simulations and the QCD phase diagram: Recent developments
CCP3-Origins-2020-08 DNRF90NT@UW-20-27
EPJ manuscript No. (will be inserted by the editor)
Complex Langevin and the QCD phase diagram: Recentdevelopments
Felipe Attanasio , Benjamin J¨ager , and Felix P.G. Ziegler Department of Physics, University of Washington, Box 351560, Seattle, Washington 98195-1560, USA CP3-Origins & Danish IAS, Department of Mathematics and Computer Science, University of Southern Denmark, Campusvej55, 5230 Odense M, DenmarkReceived: date / Revised version: date
Abstract
In this review we present the current state-of-the-art on complex Langevin simulations and theirimplications for the QCD phase diagram. After a short summary of the complex Langevin method, wepresent and discuss recent developments. Here we focus on the explicit computation of boundary terms,which provide an observable that can be used to check one of the criteria of correctness explicitly. We alsopresent the method of Dynamic Stabilization and elaborate on recent results for fully dynamical QCD.
PACS.
Strongly coupled quantum matter encompasses some ofthe most interesting problems in modern physics. MonteCarlo methods are commonly used to perform numericalsimulations of theories not amenable to perturbative ex-pansions. These methods typically rely on the path in-tegral formulation of the theory in Euclidean space-timeto have Boltzmann-like weights, which can be interpretedas probability distributions. This allows the generation offield configurations distributed according to these weights,whence observables can be sampled.Real-time theories, QCD at finite baryon density, andnon-relativistic bosons, amongst others, however, have com-plex actions and, therefore, complex weights in their pathintegrals. This forbids a probabilistic interpretation of thepath integral measure. Moreover, this poses a big numer-ical challenge as oscillatory contributions from the gener-ated configurations must cancel precisely in order to giveaccurate answers. This is known as the sign problem .In this work, we focus on the complex Langevin (CL)method. It is an extension of stochastic quantisation [1],where real dynamical variables are allowed to take com-plex values. After some early works [2, 3], it was realisedthat the method was plagued by runaway solutions andconvergence to wrong limits [4–7]. More recently, complexLangevin experienced a revival [8–14], with studies focus-ing on understanding its properties and why it sometimesfailed [15–17]. This progress led to the use of adaptive stepsize for the numerical integration [18], which has improvedthe numerical stability and reduced the problem of run-away solutions. In addition, criteria for correctness [19,20]that allow a posteriori checks were formulated. Anotherbyproduct of the resurgence of complex Langevin was the invention of the gauge cooling method [21, 22], inspiredby gauge fixing [10], to limit excursions on the complexmanifold in a gauge invariant way. Many more studies fol-lowed, applying the method to SU(3) spin models [23], theThirring model [24, 25], random matrix theories [26–28],and QCD with staggered quarks [29], with a hopping ex-pansion [30], and in the limit of heavy-dense quarks [31].Additional uses of complex Langevin outside of QCDrelated models include superstring-inspired models [32]and quantum many-body studies: rotating bosons [33,34],spin-orbit coupled bosons [35], fermions with repulsive in-teractions [36], non-zero polarisation [37, 38], mass imbal-ance [39,40], and to determine their virial coefficients [41].
The goal is to generate field configurations with a complexmeasure e − S D φ ≡ ρ D φ (1)where φ generically represents all fields in the theory, and S is a complex action defined on a real manifold M . Thismeasure is replaced by a positive one, P D φ R D φ I , definedon the complexification M c of M . This is the equilibriummeasure of the Langevin process on M c . When the criteriafor convergence, outlined in [15, 19], are met observablescalculated with either measure have the same expectationvalue.The Langevin process is given by dφ R = K R dt + dw , (2) dφ I = K I dt , (3) a r X i v : . [ h e p - l a t ] M a y Felipe Attanasio et al.: Complex Langevin and the QCD phase diagram: Recent developments where t is known as the Langevin time, and dw is a Wienerprocess normalised as (cid:104) dw (cid:105) = 2 dt . The drifts are givenby K R = − Re [ ∇ φ S [ φ R + iφ I ]] , (4) K I = − Im [ ∇ φ S [ φ R + iφ I ]] . (5)The choice to add the Wiener process only for the realfield is arbitrary and can be generalised. However, studieshave shown that it is more beneficial to add noise onlyto the real part of the forces, see e.g. [42]. The processis said to produce correct results if the expectation valuefor a generic observable O for asymptotically long times, (cid:104)O(cid:105) ∞ , agrees with the ‘correct’ expectation value, calcu-lated with the original complex measure (cid:104)O(cid:105) ∞ = (cid:104)O(cid:105) c ≡ (cid:90) D φ O ( φ ) e − S . (6)The complexification of the original manifold implies mak-ing all fields complex-valued. In the case of gauge theoriesthis means relaxing the unitary constraint of the gaugelinks, thus allowing the full space of non-singular matricesto be explored. For QCD this implies that the standardcolour group SU(3) is extended to SL(3 , C ), which is nota compact group. Group elements can be arbitrarily faraway from SU(3), thus contradicting the criteria for cor-rectness. The method of gauge cooling [21, 22] was con-structed as a means to bring the evolution closer to theunitary manifold in a gauge invariant way.Gauge cooling consists of a series of gauge transform-ations constructed to reduce the distance to SU(3) in asteepest descent fashion. A recent discussion on how gaugecooling stabilises complex Langevin dynamics can be foundin [43]. In that work, the authors carry out analytical andnumerical studies of the effects of gauge cooling on SU(2)and SU(3) Polyakov chains. They point out four main ef-fects of gauge cooling: (i) the removal of a large numberof redundant degrees of freedom, (ii) some components ofthe drift no longer have any effect on the dynamics, (iii)emergence of additional drift terms towards the unitarymanifold supporting stability of the simulation, and (iv)the introduction of singularities in the drift. (i) - (iii) sta-bilise the Langevin dynamics and support requirementsfor the criteria of correctness, while implications of (iv)remain still unknown. Correct convergence of the complex Langevin techniqueholds if the expectation value of an observable O over theprobability distribution P ( φ R , φ I ; t ) agrees with that cor-responding to the time-evolved complex measure ρ ( φ ; t ) (cid:104)O(cid:105) P ( t ) = (cid:104)O(cid:105) ρ ( t ) . (7)The important step in the proof of correctness is the in-troduction of an interpolation function connecting the leftand right-hand side of (7) F O ( t, τ ) = (cid:90) D φ R D φ I P ( φ R , φ I ; t − τ ) O ( φ R , φ I ; τ ) . (8) Here the interpolating parameter τ is such that 0 ≤ τ ≤ t .For τ = 0 (8) reproduces the LHS of (7) and for τ = t the RHS is recovered assuming P ( φ R , φ I ; 0) = ρ ( φ ; 0), seeequation (27) in [19] for a derivation. Correct convergenceof the complex Langevin process requires the interpolat-ing function to be τ -independent. This can be proven tohold in the absence of boundary terms arising from theintegral (8) (most prominently in φ I -direction). This willbe further discussed in the next section.The formal argument of correctness relies on the holo-morphicity of the action and hence of the drift. However,the justification of correctness can be extended also to thepresence of meromorphicity. For instance, in QCD zerosof the fermion determinant give rise to poles in the drift.Correctness can be ensured if the distribution P vanishessufficiently fast close to poles [44,45]. In practice, this canbe verified a posteriori . It should be noted, however, thata pole lying inside the distribution P can lead to a separ-ation of configuration space and to ergodicity problems.Complex Langevin studies of chiral random matrixtheory [26, 46] and of effective Polyakov line models [47]have identified the branch cut of the logarithm in the ef-fective action as a source of failure of correct convergence.The ambiguity of the complex logarithm induces a wind-ing of the CL trajectory around a pole in the drift. How-ever, simulations of full QCD indicate that the winding isnot relevant for the correctness of the method there [45].The criteria for correctness have to be checked forevery observable. Care has to be taken in the presenceof poles. Results are affected by the interplay of the poleorder and the behaviour of the distribution P and theobservables around the pole.In [20, 48] the criteria of correctness are formulatedfrom a slightly different but equivalent view point. There,the authors point out that the exponential (power-law)fall-off behaviour of the probability distribution associ-ated with the drift indicates directly if CL works (fails).It is argued that this constitutes a necessary and sufficientcriterion. The proof of the key identity (6) is facilitated interms of two ingredients. (i) formulating the correctnesscriteria (7) in terms of a discretized Langevin time ex-pansion ( ε -expansion) of the time-evolved observable andthe probability distribution P ( φ R , φ I ; t ), and (ii) relaxingthe assumption that the radius of convergence τ in thecorresponding series expansion of the time evolution (8)is infinite to being finite. Correctness of the CL processfollows if the ε -expansion is valid. The latter stands andfalls with the exponential decay of the drift histogram.Moreover, it is interpreted that – which is relevant for thenext section – the appearance of boundary terms arisingfrom an integration by parts in the τ -derivative of the in-terpolation function is related to the breakdown of the ε -expansion. The criterion is probed in numerical simu-lations using simple models [20], gauge theories [48] andfor full QCD [49]. For a comparison of the drift criterionwithin a recent study on correctness in terms of boundaryterms see [50]. elipe Attanasio et al.: Complex Langevin and the QCD phase diagram: Recent developments 3 The condition of (8) being τ -independent can be writtenas ∂∂τ F O ( t, τ ) = lim Y →∞ B O ( Y ; t, τ ) = 0 , (9)where Y is a cutoff in the non-compact directions and themiddle term is a boundary term left from the integrationby parts B O ( Y ; t, τ ) = (cid:90) D φ R [ K I ( φ R , Y ) P ( φ R , Y ; t − τ ) O ( φ R + iY ; τ ) − K I ( φ R , − Y ) P ( φ R , − Y ; t − τ ) O ( φ R − iY ; τ )] . (10)This term has been recently investigated in [50] for asimple one plaquette model with U(1) symmetry. In thatmodel, (10) is non-zero and spoils the correctness of ex-pectation values. Moreover it has been shown that theboundary term can be determined stochastically from theLangevin process. Comparisons with numerical solutionsof the Fokker-Planck equation (FPE), which describes theevolution of P , have been performed. This is however verydifficult in higher dimensions. It is shown that the additionof a regulator term to the action is capable of reducing theboundary terms.In [51], studies of boundary terms have been deepenedin the U(1) one plaquette model and for the Polyakovchain. This has been also extended to field theories suchas the three-dimensional XY model as well as HDQCD.In the case of gauge theories, the boundary terms appearto be related to the distance to the unitary manifold, typ-ically measured via the unitarity norm [21]. Moreover, ithas been shown that boundary terms provide an estimateof the error of the Langevin process compared to the cor-rect results. The quantification of the latter relies on theassumption that the boundary term for the observablesof interest is maximal at τ = 0. From the analysis of theFPE this was shown to hold in the U (1) one-plaquettemodel for some classes of observables. The error estima-tion demands the calculation of ‘higher order’ boundaryterms which are numerically more expensive. The useful-ness of this measure is currently seen to depend on themodel being simulated and is under further investigation.There are indications that the aforementioned assumptionfor the τ = 0 behaviour is not valid for HDQCD.Table 1 summarises results from complex Langevinsimulations of HDQCD from [51]. It is worth noticingthat the magnitude of B is typically an order of mag-nitude larger than the error across a wide range of inversecouplings. The errors were computed by comparing theCL results with reweighting simulations. Table 2 showsresults for the three-dimensional XY model, obtained us-ing the dual worldline formulation [52] and with complexLangevin. The CL results were corrected using the bound-ary term analysis in [51]. The explicit computation of theboundary terms may hence serve to correct simulationswith non-zero contribution from them. Table 1.
Boundary terms and the error of the complexLangevin simulations for the spatial plaquette in a HDQCDsimulations with volume 6 , µ = 0 . N F = 1, and κ = 0 . β B CL error5 . − . . . − . . . − . − . . − . − . Table 2.
Comparison of the worldline and the corrected com-plex Langevin simulations for the three-dimensional XY model,as shown in [51]. β µ worldline corrected CL0 . − − . − . . . − . − . . . − . − . . − − . − . . . − . − . . . − . − . In state-of-the-art full QCD simulations the boundaryterms have to be monitored. In order to keep them smallthe CL trajectory is cut off when a prescribed threshold ofthe unitarity norm is reached [53]. Moreover, decreasing β can cause an increase in the boundary terms. In ref. [31] it was noticed that for some simulations theLangevin process would initially converge to the correctvalue (comparing with reweighting, when applicable) andthen slowly drift to an incorrect one, despite the use ofgauge cooling. The departure from the correct result al-ways coincided with the increase in the distance of the fieldconfigurations from the unitary manifold. The method ofdynamic stabilisation (DS) was then proposed. The ideais to add a term to the Langevin drift itself, which aimsto (a) be small in comparison to the drift originating fromthe action; (b) vanish in the na¨ıve continuum limit; (c)affect only the non-compact directions of fields; (d) beSU(3) gauge invariant.The additional force proposed in [54] is K x,µ → K x,µ + iα DS M x , (11)where M x is proportional to a power of the unitarity norm.This choice of force, however, is non-holomorphic and thusviolates the criteria of correctness. The real parameter α DS controls the strength of the DS term. Given the complex-ity of gauge theories, it is difficult to predict its effectson the CL simulations a priori , except in two limitingcases: when α DS is very small, the DS term will have anegligible effect and the dynamics should remain essen-tially unchanged; conversely, for large values of the con-trol parameter the DS force heavily suppresses excursions Felipe Attanasio et al.: Complex Langevin and the QCD phase diagram: Recent developments
Table 3.
Data comparing HMC vs. complex Langevin (CL)simulations at vanishing chemical potential for two observ-ables: plaquette and chiral condensate ψψ . These simulationsused four flavours of na¨ıve staggered fermions with β = 5 . a m q = 0 .
025 as shown in [54].plaquette ψψ Volume HMC CL HMC CL6 . . . . . . . . . . . . . . . . into the non-unitary directions of SL(3 , C ), thus effectivelyre-unitarising the theory. The optimal values for α DS arefound in a region where expectation values of the observ-ables are least sensitive to it. Calculating the boundaryterms explicitly could provide a non-heuristic way of de-termining optimal values of α DS .Agreement between complex Langevin and HMC atzero chemical potential can be found when Dynamic Sta-bilisation is used. Table 3 shows the accuracy that canbe achieved. At vanishing chemical potential there is nosign problem, so that CL and HMC simulations shouldagree. However, the bi-linear noise scheme (see, e.g., [29])provides real drifts only on average. Thus a small sourceof complexity is always present, even though the theory it-self is real. Without Dynamic Stabilisation this eventuallyleads to the failure of CL. When the quark chemical potential is sufficient for theformation of bound states, the fermion determinant ex-hibits zeros. When the Langevin process explores regionsaround those zeros the drift becomes near-singular, lead-ing to unstable simulations. In ref. [55] the fermion matrixis changed by the addition of a term iαψ ( x )( γ ⊗ γ ) ψ ( x )to the Lagrangian density, where the γ ’s act on spinorand flavour indices, respectively. This deformation, for α large enough, moves the eigenvalue distribution of the fer-mion matrix away from zero.Studies have been performed using a lattice of volume4 ×
8, coupling β = 5 .
7, quark mass am = 0 .
05, andchemical potential 0 . ≤ aµ ≤ .
7. Histograms of theLangevin drift show a power-law behaviour for α < . α ∼ .
6. These observationsindicate that an extrapolation from the deformed to theoriginal manifolds should only use points simulated with0 . < α < .
6. Extrapolated values for the number densityand chiral condensate simulated with CL were comparedto RHMC simulations in the phase-quenched ensemble asa function of the chemical potential. Both observables alsoshow a steeper dependence on µ in the CL simulations,which is qualitatively consistent with the expectation that,in the thermodynamic limit at zero temperature, physical *8, Symanzik, β =3.72-stout. N F =4, m=0.02 T ≈
260 MeVm π ≈
525 MeV ∆ ( p / T ) µ /TTaylor exp. 6th orderTaylor exp. 4th orderTaylor exp. 2nd orderCLE6th order fit Figure 1.
The pressure difference as a function of the chemicalpotential as shown in [56]. The simulation parameters are listedin the figure. observables are independent of µ for µ < m N /
3, for fullQCD, or µ < m π /
2, in the phase-quenched case.
Simulations of the QCD phase diagram with fully dynam-ical quarks are underway. Results have been reported forhigh [56], low [49], and zero [57] temperature regimes. Theformer two studies used four flavours of staggered quarks,while the latter used two flavours. Additionally, a study ofthe deconfinement transition in QCD with heavy pions canbe found in [58], where two flavours of Wilson quarks havebeen considered. All works have employed gauge coolingto reduce large explorations of non-unitary directions.
In ref. [56] the complex Langevin simulations were en-hanced by stout smearing [59] to smooth the gauge con-figurations. To achieve this, the smearing procedure hadto be generalized to SL(3 , C ) link variables. In addition,tree-level Symanzik improved gauge action was used fora volume of 16 ×
8. This allows a comparison betweenthe standard Wilson plaquette action and the improvedsetup. The quarks are heavier than in nature to keep thesimulations numerically feasible, resulting in a pion massin the range of 500 to 700 MeV. A comparison betweenstandard Taylor expansion and complex Langevin simula-tion allows an additional consistency check. Figure 1 (leftpanel of figure 7 in [56]) shows the pressure difference usinga Taylor expansion up to 6th order as well as direct resultsfrom complex Langevin simulations. The agreement is re-markably good, so that complex Langevin can be used todetermine thermodynamic quantities, especially at hightemperatures. elipe Attanasio et al.: Complex Langevin and the QCD phase diagram: Recent developments 5
In refs. [49, 60] the authors study a lower temperaturesetup with β = 5 . am = 0 .
01 for two lattice volumes:8 ×
16 and 16 ×
32. In order to ensure reliability of theirresults, they used the criterion based on the distributionof the Langevin drift [20], described at the end of section2.2. They found that the average quark number exhib-its a plateau at (cid:104) N (cid:105) = N f N c N s = 24 as a function ofthe chemical potential. The observed plateau was under-stood qualitatively in terms of a picture of free fermionsat zero temperature: due to the discreteness of the lat-tice momenta, N f N c N s is the maximum number of zeromomentum quarks that can exist until the chemical poten-tial is large enough to excite the first non-zero momentumstates. This interpretation is possible due to the smallnessof the gauge coupling considered, making a picture of freefermions valid. For a larger volume, they found that theplateau shifts to smaller values of µ . This is expected, asthe discrete momentum states become closer. In a complementary study of QCD at zero temperatureand finite chemical potential the authors of ref. [57] useda setup of β = 5 . and β = 5 . , both cases with quark mass am = 0 .
025 and two fla-vours of staggered quarks for their CL simulations, aug-mented with gauge cooling. They found that, in general,simulations at smaller coupling ( β = 5 .
7) produced res-ults closer to the correct ones, when those were available,also reporting that the average unitarity norm decreasesas the continuum limit is approached, in accordance withthe findings of [31]. However, the expected transition fromhadronic to nuclear matter at µ ≈ m N /
3, where m N isthe nucleon mass, is not seen. No signs of new, exoticphases of matter, such as colour-superconductors werefound for µ ≥ m N /
3, and there is no indication of dif-ferences between the full theory and its phase-quenchedapproximation for aµ ≥ .
5, but not large enough wheresaturation is dominant.It is argued in ref. [57] that the discrepancy betweensimulation results and physical expectations could be dueto a number of issues complex Langevin simulations face.One possibility is that CL is known to converge to phase-quenched results in some random matrix theories [28], al-though this has not been observed in HDQCD. Another isthat CL produces correct results for non-singular observ-ables, but those used in ref. [57] do have poles.Lastly, for small temperatures and finite µ , the CLfaces severe problems because the eigenvalues of the Diracmatrix approach zero. Therefore the matrix becomes ill-conditioned and standard algorithms such as the conjug-ate gradient method break down. The latter representsthe backbone of calculating the fermionic drift force. Inref. [61] the method of selected inversion has been sugges-ted. It is a purely algebraic technique based on the LUdecomposition where a subset of the elements of the in- Table 4.
The fitted curvature and T c (0) for two flavour Wilsonfermions with N s = 12 and 16 using two different methodstaken from [58]. The curvature has the form T c ( µ ) = T c (0) − κ µT c (0) . More details on the two methods can be found insection IIIA and IIIB of [58].Method N s κ T c (0)/MeVfit B3 12 0 . . . . × − . verse matrix is computed, thus making it cheaper than afull inversion. In ref. [58] the phase diagram of QCD in the T - µ planewas investigated with two flavours of Wilson quarks, forchemical potentials up to µ ∼ T using CL. The study wascarried out with a relatively high pion mass of m π ≈ . κ of the transition line. It has also beennoticed that κ has a non-monotonic behaviour as a func-tion of the quark mass. Despite the transition being asmooth crossover for the parameters considered the crit-ical temperature can be determined relatively well usingthe Binder cumulant. A summary of the results is shownin Table 4. Recent developments on the complex Langevin methodhave shown promising progress. An important step hasbeen to go beyond the analysis of the correctness criteriain simple models to a practical approach, which is ap-plicable to field theories. With the determination of theboundary terms the error on the complex Langevin pro-cess can be estimated and, in some cases, compensatedfor. Challenges however remain in the development of thisnovel approach for full QCD.From a different angle, dynamic stabilization providesa viable technique for the CL to be applied to the regionsof lower temperatures and medium densities. Addition-ally, analysis of boundary terms, both numerically andanalytically, may be able to provide a consistent way ofoptimising the control parameter of DS.First results for CL simulations of the phase diagram offull QCD have recently appeared. Despite the lack of evid-ence for the expected transition from hadronic to nuclear
Felipe Attanasio et al.: Complex Langevin and the QCD phase diagram: Recent developments matter at zero temperature, results at low, but finite, tem-perate are encouraging and can be understood physically.The agreement with Taylor expansion at high temperat-ures is a major step towards the ultimate goal of simulat-ing the QCD phase diagram using the complex Langevinmethod.
We are grateful for discussion and collaboration with GertAarts and Ion-Olimpiu Stamatescu. We thank D´enes Sextyfor providing one of the figures to this manuscript. Thework of F.A. was supported by US DOE Grant No. DE-FG02-97ER41014 MOD27.
References
1. G. Parisi, Y.s. Wu, Sci. Sin. , 483 (1981)2. F. Karsch, H.W. Wyld, Phys. Rev. Lett. , 2242 (1985)3. P.H. Damgaard, H. H¨uffel, Phys. Rept. , 227 (1987)4. J. Ambjorn, S.K. Yang, Phys. Lett. B165 , 140 (1985)5. J.R. Klauder, W.P. Petersen, J. Stat. Phys. , 53 (1985)6. H.Q. Lin, J.E. Hirsch, Phys. Rev. B , 1964 (1986)7. J. Ambjorn, M. Flensburg, C. Peterson, Nucl. Phys. B275 ,375 (1986)8. J. Berges, I.O. Stamatescu, Phys. Rev. Lett. , 202003(2005), hep-lat/0508030
9. J. Berges, S. Bors´anyi, D. Sexty, I.O. Stamatescu, Phys.Rev.
D75 , 045007 (2007), hep-lat/0609058
10. J. Berges, D. Sexty, Nucl. Phys.
B799 , 306 (2008),
11. C. Pehlevan, G. Guralnik, Nucl. Phys.
B811 , 519 (2009),
12. G. Aarts, I.O. Stamatescu, JHEP , 018 (2008),
13. G. Aarts, Phys. Rev. Lett. , 131601 (2009),
14. G. Guralnik, C. Pehlevan, Nucl. Phys.
B822 , 349 (2009),
15. G. Aarts, E. Seiler, I.O. Stamatescu, Phys. Rev.
D81 ,054508 (2010),
16. G. Aarts, F.A. James, JHEP , 020 (2010),
17. G. Aarts, K. Splittorff, JHEP , 017 (2010),
18. G. Aarts, F.A. James, E. Seiler, I.O. Stamatescu, Phys.Lett.
B687 , 154 (2010),
19. G. Aarts, F.A. James, E. Seiler, I.O. Stamatescu, Eur.Phys. J.
C71 , 1756 (2011),
20. K. Nagata, J. Nishimura, S. Shimasaki, Phys. Rev.
D94 ,114515 (2016),
21. E. Seiler, D. Sexty, I.O. Stamatescu, Phys. Lett.
B723 ,213 (2013),
22. G. Aarts, L. Bongiovanni, E. Seiler, D. Sexty, I.O. Stam-atescu, Eur. Phys. J.
A49 , 89 (2013),
23. G. Aarts, F.A. James, JHEP , 118 (2012),
24. J.M. Pawlowski, C. Zielinski, Phys. Rev.
D87 , 094503(2013),
25. J.M. Pawlowski, C. Zielinski, Phys. Rev.
D87 , 094509(2013),
26. A. Mollgaard, K. Splittorff, Phys. Rev.
D88 , 116007(2013),
27. A. Mollgaard, K. Splittorff, Phys. Rev.
D91 , 036007(2015),
28. J. Bloch, J. Glesaaen, J.J.M. Verbaarschot, S. Zafeiro-poulos, JHEP , 015 (2018),
29. D. Sexty, Phys. Lett.
B729 , 108 (2014),
30. G. Aarts, E. Seiler, D. Sexty, I.O. Stamatescu, Phys. Rev.
D90 , 114505 (2014),
31. G. Aarts, F. Attanasio, B. J¨ager, D. Sexty, JHEP , 087(2016),
32. J. Nishimura, A. Tsuchiya, JHEP , 077 (2019),
33. T. Hayata, A. Yamamoto, Phys. Rev.
A92 , 043628 (2015),
34. C. Berger, J. Drut, PoS
LATTICE2018 , 244 (2018)35. F. Attanasio, J.E. Drut, Phys. Rev. A , 033617 (2020),
36. A.C. Loheac, J.E. Drut, Phys. Rev.
D95 , 094502 (2017),
37. A.C. Loheac, J. Braun, J.E. Drut, Phys. Rev.
D98 , 054507(2018),
38. L. Rammelm¨uller, A.C. Loheac, J.E. Drut, J. Braun, Phys.Rev. Lett. , 173001 (2018),
39. L. Rammelm¨uller, W.J. Porter, J.E. Drut, J. Braun, Phys.Rev.
D96 , 094506 (2017),
40. L. Rammelm¨uller, J.E. Drut, J. Braun (2020),
41. C. Shill, J. Drut, Phys. Rev. A , 053615 (2018),
42. G. Aarts, P. Giudice, E. Seiler, Annals Phys. , 238(2013),
43. Z. Cai, Y. Di, X. Dong, Commun. Comput. Phys. , 1344(2020),
44. J. Nishimura, S. Shimasaki, Phys. Rev.
D92 , 011501(2015),
45. G. Aarts, E. Seiler, D. Sexty, I.O. Stamatescu, JHEP ,044 (2017),
46. K. Splittorff, Phys. Rev.
D91 , 034507 (2015),
47. J. Greensite, Phys. Rev.
D90 , 114507 (2014),
48. K. Nagata, J. Nishimura, S. Shimasaki, JHEP , 004(2018),
49. S. Tsutsui, Y. Ito, H. Matsufuru, J. Nishimura, S. Shima-saki, A. Tsuchiya (2019),
50. M. Scherzer, E. Seiler, D. Sexty, I.O. Stamatescu, Phys.Rev.
D99 , 014512 (2019),
51. M. Scherzer, E. Seiler, D. Sexty, I.O. Stamatescu, Phys.Rev.
D101 , 014501 (2020),
52. D. Banerjee, S. Chandrasekharan, Phys. Rev.
D81 , 125007(2010),
53. M. Scherzer,
PhD thesis (Heidelberg, 2019)54. F. Attanasio, B. J¨ager, Eur. Phys. J.
C79 , 16 (2019),
55. K. Nagata, J. Nishimura, S. Shimasaki, Phys. Rev. D ,114513 (2018),
56. D. Sexty, Phys. Rev.
D100 , 074503 (2019),
57. J.B. Kogut, D.K. Sinclair, Phys. Rev.
D100 , 054512(2019),
58. M. Scherzer, D. Sexty, I.O. Stamatescu (2020),
59. C. Morningstar, M.J. Peardon, Phys. Rev.
D69 , 054501(2004), hep-lat/0311018
60. Y. Ito, H. Matsufuru, J. Nishimura, S. Shimasaki,A. Tsuchiya, S. Tsutsui, PoS
LATTICE2018 , 146 (2018),
61. J. Bloch, O. Schenk, EPJ Web Conf. , 07003 (2018),, 07003 (2018),