Complex Langevin Simulationsof nonzero density QCD

Dénes Sexty
Wuppertal University
Hadrons and Hadron interactions – Long term workshop
Kyoto, 4th of March, 2015.
Collaborators: Gert Aarts, Erhard Seiler, Ion-Olimpiu Stamatescu
Felipe Attanasio, Lorenzo Bongiovanni, Benjamin Jäger,Zoltán Fodor, Sándor Katz
1. Sign problem in lattice QCD
2. Complex Langevin equation and gauge cooling
3. Phase diagram of HDQCD
4. kappa and kappa_s expansion
5. Full QCD

Motivations
Phase diagram of QCD matter

QCD sign problem
Euclidean SU(3) gauge theory with fermions:ā D Ψ exp(−S E [ Aaμ ]− Ψ
̄ D E ( Aμa ) Ψ )Z =∫ DAμ D Ψ

Integrate out fermionic variables, perform lattice discretisation

Aμa ( x , τ) → U μ ( x , τ )∈SU (3) link variables
D E ( A) → M (U ) fermion matrix
Z =∫ DU exp(−S E [U ])det (M(U))
det ( M (U ))>0 → Importance sampling is possible

Non-zero chemical potential
For nonzero chemical potential, the fermion determinant is complex

det ( M (U ,−μ ∗ ))=(det ( M (U ) ,μ)) ∗Sign problem

Naïve Monte-Carlo breaks down

QCD sign problem
det ( M (U ,μ))∈ℂ for μ>0
Z =∫ DU exp(−S E [U ])det (M(U))
Path integral with complex weight

Only the zero density axis is directly accessible
to lattice calculations using importance samplingReweighting
det M (μ)
F
DU e det M (μ) F ∫ DU e R∫
R〈 F 〉μ ==−S∫ DU e det M (μ) ∫ DU e−S R det M (μ )
R
−S E

−S E

E

E

〈 F det M (μ)/ R 〉 R=〈 det M (μ)/ R 〉 R〈〉

R=det M (μ=0), ∣det M (μ)∣, etc.

det M (μ)
Z (μ)
V==exp − Δ f (μ , T )
R
ZR
T
R

(

)

Δ f (μ , T ) =free energy difference
Exponentially small as the volume increases〈 F 〉μ → 0/0

Reweighting works for large temperatures and small volumes
Sign problem gets hard at

μ /T ≈1

Evading the QCD sign problem
Most methods going around the problem work only for  = B /3T
(Multi parameter) reweighting

Barbour et. al. '97; Fodor, Katz '02

Analytic continuation of results obtained at imaginary μ
2Taylor expansion in (μ / T )Lombardo '00;
de Forcrand, Philipsen '02;
D'Elia Sanfilippo '09; Cea et. al. '08

de Forcrand et al. (QCD-TARO) '99; Hart, Laine, Philipsen '00;
Gavai and Gupta '03;Allton et al. '05 ; de Forcrand, Philipsen '08,...

Canonical Ensemble, denstity of states, curvature of critical surface,subsets, fugacity expansion, SU(2) QCD, G2 QCD, dual variables, worldlines, ….

A Direct Method: Complex Langevin
Use analiticity, expand integrals to the complex plane

Stochastic quantisation
Recent revival:
Aarts and Stamatescu '08
Bose Gas, Spin model, etc. Aarts '08, Aarts, James '10 Aarts, James '11
Proof of convergence:
Aarts, Seiler, Stamatescu '11
QCD with heavy quarks:
Seiler, Sexty, Stamatescu '12
Kappa Expansion:
Aarts, Seiler, Sexty, Stamatescu 1408.3770
Full QCD with light quarks: Sexty '14

Stochastic Quantization
Given an action

Parisi, Wu (1981)

S ( x)dx∂S=− 
d∂xStochastic process for x:

Gaussian noise ⟨ η( τ)⟩=0⟨ η( τ) η(τ ' )⟩=δ( τ−τ ')Averages are calculated along the trajectories:
T1⟨O ⟩=limT →∞ ∫ O( x ( τ))d τ=
T 0−S ( x ) O( x )∫edx−S( x)
edx∫

Fokker-Planck equation for the probability distribution of P(x):∂P ∂ ∂P∂S=P
=−H FP P∂ ∂ x ∂ x∂xReal action 

for real action theLangevin method is convergent

positive eigenvalues

Langevin method with complex action
Klauder '83, Parisi '83, Hueffel, Rumpf '83,
Okano, Schuelke, Zeng '91, ...
applied to nonequilibrium: Berges, Stamatescu '05, ...

The field is complexifieddx∂S=− 
d∂xreal scalarcomplex scalarlink variables: SU(N)SL(N,C)

compact

non-compact

det (U )=1,
Analytically continued observables

1
1∫ P comp ( x )O ( x )dx= Z ∫ P real ( x , y )O ( x +iy ) dx dyZ〈 x 2 〉real → 〈 x 2− y 2 〉 complexified“troubled past”: Lack of theoretical understanding
Convergence to wrong results
Runaway trajectories

U + ≠ U −1

Proof of convergence
If there is fast decay

P ( x , y )→ 0 as y →∞and a holomorphic action

S ( x)

then CLE converges to the correct result[Aarts, Seiler, Stamatescu (2009)
Aarts, James, Seiler, Stamatescu (2011)]

Non-holomorphic action for nonzero density
S =S W [U μ ]+ln Det M (μ)

measure has zeros (Det M=0)
complex logarithm has a branch cut
meromorphic drift
Is it a problem for QCD?

[see also: Mollgaard, Splittorff (2013), Greensite(2014)][QCD and poles: Aarts, Seiler, Sexty, Stamatescu 2014]

Non-real action problems and CLE
1. Real-time physics“Hardest” sign probleme[Berges, Stamatescu (2005)][Berges, Borsanyi, Sexty, Stamatescu (2007)][Berges, Sexty (2008)][Anzaki, Fukushima, Hidaka, Oka (2014)][Fukushima, Hayata (2014)]

iS M

Studies on Oscillator, pure gauge theory

2. Theta-Term

μν

μ νθρ

S =F μ ν F +i Θ ϵ

F μ ν F θρ[Bongiovanni, Aarts, Seiler, Sexty, Stamatescu (2013)+in prep.]Θ real → complex action, ⟨Q⟩ imaginaryΘ imaginary → real action, ⟨Q⟩ realOn the lattice
μ νθρ

Q=ϵ

F μ ν F θ ρ → ∑ x q( x)Not topological
Cooling is needed
Θ L bare parameter needs
renormalisation

Θ imaginary → use real Langevin or HMC
Θ real → use complex Langevin3. Bose gas with non-zero chemical potential[Aarts (2009)]

4. XY-model, SU(3) Spin model[Aarts and James (2010-2012)]

5. Bose gas in rotating frame[Hayata and Yamamoto (2014)]

6. Random matrix theory[Mollgaard and Splittorff (2014)]

This Talk:
QCD with static quarks
Hopping expansion of QCD expansion to all orders
Full QCD with light quarks[Seiler, Sexty, Stamatescu (2013)][Aarts, Bongiovanni, Seiler, Sexty, Stamatescu (2013)][Sexty (2014)][Aarts, Seiler, Sexty, Stamatescu (2014)]

Measure
on real axis

Gaussian Example
S [ x ]=σ x 2 +i λ xCLEd( x +i y )=−2 σ ( x +iy )−i λ + ηdτP ( x , y )=e22−a( x−x 0 ) −b( y− y0 ) −c (x− x 0 )( y− y 0 )]Gaussian distribution
around critical point∂ S ( z)=0∂z z

Thimble ż=−∂ z S ( z)Straight lines
starting from z0

Measureon thimble

0

σ=1+i

λ=20

Gauge theories and CLElink variables: SU(N)SL(N,C)

compact

non-compact

det (U )=1,

U + ≠ U −1

Gauge degrees of freedom also complexify

Infinite volume of irrelevant, unphysical configurations
Process leaves the SU(N) manifold exponentially fast
already at μ≪1

Unitarity norm:

Distance from SU(N)∑i Tr (U i U i+ )2+
∑ij ∣(U U −1)ij∣
Tr (U U + )+Tr (U −1 (U −1 ) + )≥2 N

Minimize unitarity norm∑i Tr (U i U i+ −1)

Distance from SU(N)
Gauge transformation atchanges 2d link variablesx

U μ ( x )→exp(−α ϵ λ a G a ( x ))U μ ( x )U μ ( x−aμ )→U μ ( x−aμ )exp(α ϵ λ a G a ( x ))Steepest descent

Dynamical steps are interspersedwith several gauge cooling steps

The strength of the cooling is determined bycooling steps
gauge cooling parameter α

Empirical observation:
Cooling is effective forβ>β mina<a max

but remember, β→∞in cont. limit

a max ≈0.1−0.2 fm

The effect of gaugecooling

Smaller coolingexcursions into complexified
manifold“Skirt” develops
small skirt gives correct result

Heavy Quark QCD at nonzero chemical potential (HDQCD)
Hopping parameter expansion of the fermion determinant
Spatial hoppings are dropped
2
Det M (μ)=∏ x det (1+C P x )2 det (1+C ' P−1)x

P x = ∏τ U 0 ( x + τ a 0 )C =[2 κ exp(μ)]NτC '=[2 κ exp(−μ)]N

S =S W [U μ ]+ln Det M (μ)
Studied with reweighting

R=e

De Pietri, Feo, Seiler, Stamatescu '07∑ x C Tr P x +C ' Tr P −1

CLE study using gaugecooling[Seiler, Sexty, Stamatescu (2012)]

τGauge cooling stabilizes the distribution
SU(3) manifold instable even at μ=0

Fermion density:

1 ∂ ln Zn=N τ ∂μaverage phase:⟨exp(2 i ϕ)⟩=⟨Det M (μ)
Det M (−μ)⟩

det (1+C P )=1+C 3 +C Tr P +C 2 Tr P −1

Reweigthing is impossible atSign problem is absent at
small or large μ

6≤μ /T ≤12 , CLE works all the way to saturation

Comparison to reweighting

6 4 lattice , β=5.9
6 4 lattice , μ=0.85
Discrepancy of plaquettes at β≤5.6
a skirted distribution developsa(β=5.6)=0.2 fm

Large lattice:phase transition clearly visiblefor β>β min

Mapping the phase diagram[Aarts, Attanasio, Jäger, Seiler, Sexty, Stamatescu, in prep.]

fixed β=5.8 → a≈0.15 fm
κ=0.12
onset transition at μ=−ln (2 κ)=1.43
N t∗83 lattice
N t =2.. 28

Temperature scanning

Exploring the phase diagram of HDQCD

Onset in fermionic density
Silver blaze phenomenon

β=5.8 κ=0.12 N f =2 N t =2. ..24

Polyakov loop
Transition to deconfined statePolyakov loop susceptibility

Hint of first order deconfinement and first order onset transitionκ Expansion using the loop expansion
M =1−κ Q=1−R−κ s S

Wilson fermionsS=∑i 2 Γ−i U i ( x) δ y , x +i +2 Γ+i U −1i ( y ) δ y , x−i

Spatial hoppings

R=2 κ e μ Γ−4 U 4 ( x )δ y , x+4 +2 κ e−μ Γ+4 U −1
4 ( y )δ y , x−4

Temporal hoppings

(

n

)

(

ls

n
s
Det M =exp(Tr ln M )=exp −Tr ∑ κ Q =exp −Tr ∑C κ Lcnssum for distinct paths= ∏C det (1−κl Lc )Static limit

κ →0 , μ →∞ , ζ=2 κ e μ=const
Only Polyakov loops contributeWilson fermionsΓ+ν Γ−ν =0)

no backtracking

Caculation of the first few orders
Is possible using loop expansion

with full gauge action[Bender et al. (1992)][Aarts et al. (2002)][De Pietri, Feo, Seiler, Stamatescu (2007)]κ 2 corrections

with strong coupling expansion[Fromm, Langelage, Lottini, Philipsen (2012)][Greensite, Myers, Splittorff (2013)][Langelage, Neuman, Philipsen (2014)]

κ 4 corrections

expansions with complex Langevin[Aarts, Seiler, Sexty, Stamatescu 1408.3770]

κ expansion
M =1−κ Q=1−R−κ s Sn(n
Det M =exp(Tr ln M )=exp −Tr ∑ κ Q
n)Contribution to Drift term:

K μ , x , a=Tr ( ∑ κ Qnoise vector ηK μ , x , a=η ∗ D μ , x , a Q s

κ s expansionnn−1

Dμ , x , a Q )with s=−∑ κ n Q n−1 η

(

n

κs S
κs S nDet M =Det (1− R) Det 1−= Det (1−R )exp −Tr ∑1−R
n (1−R )n(Contribution to Drift term:)using noise vector
analitically (same as LO HDQCD))

κ expansionn(n
Det M =exp(Tr ln M )=exp −Tr ∑ κ QnNumerical cost:

No poles!

N multiplications with Q

Q=R + κ S with R + ∝ eμ

κ s expansion

)

(

n

κs S
κs S nDet M =Det (1− R) Det 1−= Det (1−R )exp −Tr ∑1−R
n (1−R )n
Numerical cost:(μbad convergence at high))N multiplications with S and (1−R)−1

Temporal part analytically

better convergence properties

Calculation of high orders of corrections is easy
Explicit check of the convergence to full QCD
Convergence to full QCD with no poles

non-holomorphicity of the QCD action
is not a problem

Converges at low temperatures(large lattices)

Convergence radiusκ≈0.14 ?

Extension to full QCD with light quarks [Sexty (2014)]
Z =∫ DU e

QCD with staggered fermionsM ( x , y )=m δ( x , y )+ ∑ν

−S G

det M

ην δ μ−δ μ−1(e U ν ( x )δ( x + a ν , y )−eU ν ( x−a ν , y )δ( x−a ν , y ))
2 aν
ν4

ν4

Still doubling present N_F=4

Z =∫ DU e−S G(det M )N F/4

NF
S eff =S G −ln det M
4

Langevin equation
U '=exp ( i λ a (−ϵ Da S[U]+ √ ϵ ηa ) ) U

Drift term:−D a S [U ]= K G + K F

K Gax ν =−D ax ν S G [U ]
K

F
ax ν

NF
NF
=
D ax ν ln det M =Tr ( M −1 M ' ν a ( x , y , z ))44
M ' ν a ( x , y , z )=D a z ν M ( x , y)Extension to full QCD with light quarks[Sexty (2014)]
QCD with fermionsZ =∫ DU e

−S G

det M

Additional drift term from determinant

K

F
ax ν

NF
NF
=
D ax ν ln det M =Tr ( M −1 M ' ν a ( x , y , z ))44

Noisy estimator with one noise vector
Main cost of the simulation: CG inversion
Inversion cost highly dependent on chemical potential
Eigenvalues not bounded from below by the mass
(similarly to isospin chemical potential theory)

Unimproved staggered and Wilson fermions
Heavy quarks: compare to HDQCDLight quarks: compare to reweightingZero chemical potential
Drift is built from random numbers

real only on average

Cooling is essential already for small (or zero) muCLE and full QCD with light quarks [Sexty (2014)]
Physically reasonable results

Non-holomorphic action
poles in the fermionic drift
Is it a problem for full QCD?

So far, it isnt:
Comparison with reweighting
Study of the spectrum
Hopping parameter expansion

Comparison with reweighting
for full QCD
Reweighting from ensemble at

R=Det M (μ=0)[Fodor, Katz, Sexty (in prep.)]

Sign problem

Sign problem gets hard around⟨exp(2 i ϕ)⟩=⟨μ /T ≈1−1.5

det M (μ)
det M (−μ)⟩

Spectrum of the Dirac Operator

N F =4 staggered

Massless staggered operator at μ=0 is antihermitianSpectrum of the Dirac Operator

N F =4 staggered

Conclusions
Direct simulations of QCD at nonzero density using complexified fields
Complex Langevin Equations
Recent progress for CLE simulations
Better theoretical understanding (poles?)Gauge cooling
Kappa expansion
Two novel implementations with CLE: kappa and kappa_s
Calculations at very high orders are feasible
Convergence checked explicitly
Shows that poles give no problem in QCD
Phase diagram of HDQCD mapped out
First results for full QCD with light quarks
No sign or overlap problem
CLE works all the way into saturation region
Comparison with reweighting for small chem. pot.
Low temperatures are more demanding

Backup slides

Conclusion
QCD = HQCD for quark mass > 4/a(For large mass) HQCD is qualitatively similar to QCD

Phasequenched vs full−S
Z=∫ dU e |det M|gin phasequenched P=P−1
in full theory, inv. Polyakov loop rises first
Reweighting form PQ theory better than Reweighting from μ=0 ?

Nonzero value when:colorless bound states
formed with P or P'
1 quark:meson with P'
2 quark:
Baryon with P

P' has a peak before P
Large chemical potential: all quark states are filled
No colorless state can be formed
P and P' decays again

Spectrum of the Dirac Operator
Large chemical potential, towards saturation
Fermions become “heavy”