Real-time methods for spectral functions

In this paper we develop and compare different real-time methods to calculate spectral functions. These are classical-statistical simulations, the Gaussian state approximation (GSA), and the functional renormalization group (FRG) formulated on the Keldysh closed-time path. Our test-bed system is the quartic anharmonic oscillator, a single self-interacting bosonic degree of freedom, coupled to an external heat bath providing dissipation analogous to the Caldeira-Leggett model. As our benchmark we use the spectral function from exact diagonalization with constant Ohmic damping. To extend the GSA for the open system, we solve the corresponding Heisenberg-Langevin equations in the Gaussian approximation. For the real-time FRG, we introduce a novel general prescription to construct causal regulators based on introducing scale-dependent fictitious heat baths. Our results explicitly demonstrate how the discrete transition lines of the quantum system gradually build up the broad continuous structures in the classical spectral function as temperature increases. At sufficiently high temperatures, classical, GSA and exact-diagonalization results all coincide. The real-time FRG is able to reproduce the effective thermal mass, but overestimates broadening and only qualitatively describes higher excitations, at the present order of our combined vertex and loop expansion. As temperature is lowered, the GSA follows the ensemble average of the exact solution better than the classical spectral function. In the low-temperature strong-coupling regime, the qualitative features of the exact result are best captured by our real-time FRG calculation, with quantitative improvements to be expected at higher truncation orders.

In this paper we develop and compare different real-time methods to calculate spectral functions.These are classical-statistical simulations, the Gaussian state approximation (GSA), and the functional renormalization group (FRG) formulated on the Keldysh closed-time path.Our test-bed system is the quartic anharmonic oscillator, a single self-interacting bosonic degree of freedom, coupled to an external heat bath providing dissipation analogous to the Caldeira-Leggett model.As our benchmark we use the spectral function from exact diagonalization with constant Ohmic damping.To extend the GSA for the open system, we solve the corresponding Heisenberg-Langevin equations in the Gaussian approximation.For the real-time FRG, we introduce a novel general prescription to construct causal regulators based on introducing scale-dependent fictitious heat baths.Our results explicitly demonstrate how the discrete transition lines of the quantum system gradually build up the broad continuous structures in the classical spectral function as temperature increases.At sufficiently high temperatures, classical, GSA and exact-diagonalization results all coincide.The real-time FRG is able to reproduce the effective thermal mass, but overestimates broadening and only qualitatively describes higher excitations, at the present order of our combined vertex and loop expansion.As temperature is lowered, the GSA follows the ensemble average of the exact solution better than the classical spectral function.In the low-temperature strong-coupling regime, the qualitative features of the exact result are best captured by our real-time FRG calculation, with quantitative improvements to be expected at higher truncation orders.

I. INTRODUCTION
Non-perturbative studies in thermal field theory are most commonly done within the imaginary time formalism [1], where time is analytically continued to replace t → −iτ with the Euclidean time variable τ ∈ [0, β], where β = 1/T is the inverse temperature of the equilibrium system.Along the imaginary-time axis one introduces (anti-)periodic boundary conditions for the fields.This is the basis for the powerful path-integral representation of the partition function (see standard texts, e.g.[2,3]).If one is interested in real-time quantities such as response functions or transport coefficients, however, this analytic continuation has to be undone.Because the Euclidean time interval is compact at finite temperature, the information is incomplete without further input, and the underlying Wick rotation cannot simply be inverted by back-substituting τ → it.But even at vanishing temperature, this typically leads to ill-posed inverse numerical problems.
Methods working directly in real-time have therefore become increasingly popular over the last decades, despite their generally higher computational requirements.They avoid the analytic-continuation problem and allow studying off-equilibrium systems.The price is that standard Monte-Carlo simulations for ab-initio studies as in Euclidean spacetime are not possible, since importance sampling fails due to the sign problem.Other computational methods are needed for real-time computations.
Here we are particularly interested in real-time methods for spectral functions, which generally contain the complete spectrum of quasi-particle, multi-particle and collective excitations of a system contributing to a given correlation function, often including transport coefficients in particular low-frequency limits as relevant e.g. for hydrodynamic descriptions.Spectral functions can be formally defined as expectation values of the unequal-time commutators of the corresponding fields.After Fourier transformation, they represent the density of states containing all possible excitations in that channel.Under certain circumstances functional equations for originally Euclidean correlation functions such as Dyson-Schwinger [4][5][6][7] or FRG flow equations [8][9][10][11][12] can be analytically continued back to the real-frequency domain before they are being solved.For phenomenological scattering theory and resolvent based real-time approaches to calculate spectral functions, see e.g.Refs.[13][14][15].The flexibility to include dissipative or diffusive dynamics and the different dynamic critical behavior have so far been beyond these approaches, however.
In this work, we therefore focus on genuine realtime methods as more general and versatile alternatives.In particular, we employ three different nonperturbative real-time methods for calculating spectral functions in a simplified test system consisting of a single self-interacting bosonic degree of freedom, the quartic anharmonic oscillator, which can also be considered as the λφ 4 theory of self-interacting real scalar fields in (0 + 1) dimensions, coupled to an external heat bath which is modeled as an ensemble of harmonic oscillators after Caldeira and Leggett [16][17][18][19].This rather simple system has the advantage that its spectral function can be calculated exactly using a standard discretization of the Schrödinger equation, which allows to qualitatively discuss and quantitatively compare in detail the different approaches and approximation schemes.
This paper is organized as follows.In Section II we start by briefly summarizing the discretization scheme and formulas needed for an exact calculation of the spectral function.The following three Sections III to V describe in detail the three real-time calculation methods for the spectral function, i.e. the classical-statistical simulations in Section III, the Gaussian state approximation in Section IV, and the real-time FRG in Section V. Our results from the different methods are presented, compared and discussed in detail in Section VI, and our conclusions are given with a brief outlook on possible further studies in Section VII.Several appendices are added with further technical details and derivations especially for the Heisenberg-Langevin equations of motion in the GSA, and our regulators and truncation scheme for the real-time FRG flows.

II. QUARTIC ANHARMONIC OSCILLATOR
We consider a single quartic anharmonic oscillator of unit mass, defined by the Hamiltonian in thermal contact with an external heat bath in equilibrium.The heat bath is modeled as an ensemble of harmonic oscillators, which is generally known as the Caldeira-Leggett model in the literature [16][17][18][19], and which we will explain in more detail in Section IV.
The anharmonic oscillator can be interpreted as a selfinteracting single-component real scalar field theory in (0 + 1) dimensions.It serves here as a benchmark system for comparing different methods for calculating spectral functions, since its spectrum can be numerically determined with essentially arbitrary precision using a discretization of the Schrödinger equation on a lattice.Therefore, we will refer to the Schrödinger discretization method as the exact-diagonalization solution.
The spectral function is defined as the real distribution given by the thermal expectation value of the commutator of two Heisenberg (field) operators taken at unequal times as follows [42,43], where the average is taken over the canonical ensemble e −β Ĥ /Z at temperature T = 1/β with the parti-tion function Z = Tr e −β Ĥ .To obtain a real distribution also in the frequency domain from the real and odd ρ(−t) = −ρ(t), a factor of 2πi is commonly absorbed in the definition of its Fourier transform, which is then positive, ρ(ω) ≥ 0, for ω > 0, also odd ρ(−ω) = −ρ(ω), and normalized according to Without dissipation, this spectral function may be expressed as a sum over energy-eigenstates (see for example Chapter 6.2 of Ref. [2]), (2.5) For the non-interacting theory (λ = 0), i.e. the harmonic oscillator with frequency ω 0 , this reduces to When the free oscillator is coupled to an Ohmic heat bath, one describes an open quantum system and the retarded Green function G R 0 acquires an additional damping term with constant γ > 0, corresponding to G R 0,γ −1 = −(ω 2 − ω 2 0 + iγω), and the spectral function becomes This of course describes the collisional broadening of the free spectral function ρ 0 (ω) due to the Ohmic heat bath, together with the frequency shift of the damped harmonic oscillator from the poles in G R 0,γ (ω) at (2.8) The interaction with the environment in the open system introduces the new parameter γ which essentially quantifies how strongly the single quartic oscillator (2.1) couples to the degrees of freedom in the external heat bath.The damping γ then causes the spectral function to be broadened even in the harmonic case (with λ = 0), because the system particle can decay into heat-bath excitations, corresponding to a finite lifetime ∼ 1/γ.Hence, the thermal expectation value in (2.2) is understood over a reduced density operator where interactions with the environment x R l 2 E L t q C K M 1 V C c q R m p p C L 3 z + n s n g i Y O F s 5 U n R W X B c L A Y 4 y 4 p u a 0 F X V q 3 F I B J I T 2 + x g o T Q 9 U k A Q I q Z / h k 6 q g x + g p f q z f 5 L R e 6 I e C L d i f r 9 Q d H 2 b s B / t v J e 9 n K O u j G T i 7 2 G r 9 I q d l c x v m Y o N a O 8 q x 2 Y 0 + N 4 0 x A a J G 5 h Z q y G Z 3 A K L q K S r B j v 1 p L w K 8 i U u J K m / g p h 1 f o n x m e S m u X s o i R k r q p 3 e Y S + E + u k B E u t C h x 4 j a o 0 q Z e W 4 O 5 6 m j s u a r n D h R b z 1 X N B X Y a J 0 l w y Q 0 w J 5 b R o c z w + D T M p t R Q 5 q J w L a L g i m k p q S o 9 u Q R W h 1 E + X n m 4 k 5 O 3 + 2 R q Y y P w B 4 e 1 C 7 F x C R U p b V z W 1 C d J H a e C z M C o g 9 4 h S L I i w m b R + F A h q t g v j P q x s k x X s A Z x 1 2 9 p n s J 8 J w + + 0 w 9 J 0 9 a 6 o e V S Q E y N 6 z A g f A s n I 0 b Q 2 s c L 4 B a i z I d 5 D g t M p q l 8 9 1 P M 9 U T o K z D b V C o W 2 R D i u d z e B P 6 / c 9 r v 5 R 9 6 2 d d + 5 / j o 5 n B 2 0 Q v 0 E r 1 G O f q I j t E X d I K G i K H v 6 G e j E f + l R n O / 2 W t m 6 9 C d x k 3 O c 7 R h z c F v q n g i C g = = < / l a t e x i t > |ni $ |n + 1i < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 j L K c g K X Z 2 Z Q q n Q 9 4 Q 7 S C 4 I r G D j u B J w 3 B q g s B Z y V s 8 + J P 7 s A Y 7 l W 3 9 y i g Z G k E 8 V r z q i L 0 H h b E K W 5 q k A 5 0 l B T 6 i t f v G c y e O L g y t n a k 7 K 2 Y D h Y j H G P V N w 2 g i 6 s W w j A p J Q e X 2 O F i a F q k g A B t T N 8 M n X U G H 2 J r 9 X b / V s u 9 E L A 4 + 1 u n v W P D v P 9 I / y 3 U 2 T 5 0 r r o x k 7 G O 6 1 f p N J s L u N 8 T F B r h 0 X e u J G n x n E m I H T I 3 E J D 2 Y x O Y B h d R S X Y k V + u J e D X E a l w r U 3 8 l M N L 9 M 8 M T 6 W 1 C 1 n G S E n d 1 K 5 z C f w n V 8 o I l 1 p U O H F 3 q M q m X m u D u f p w 5 L l q 5 g 4 U W 8 1 V z w V 2 G i d J c M U N M C c W 0 a H M 8 P g 0 z K b U U O a i c B 2 i 4 J J p K a m q P L k A 1 o R h M V p 6 u F u Q d 7 t k a m M j 8 H s H j Q u x c Q U 1 q W x c 1 t Q n S R 2 n g s z A q L 3 s A C R Z E u F u 0 f h Q I e r Y L w z 7 s b J M V 7 A C c c + v a Z 7 C f L c I v t s P S d P O q q H l U k B M j e s w I H w H J y N G 0 M b H C + A W o s w H R Q F X m E x T + d 6 n m O u J 0 J d g 1 q l U L L I h x H O 5 v Q n 8 f + e 0 n x U f s v x r v 3 t 8 e H M 4 m + g l e o X e o A J 9 R M f o C z p B A 8 T Q d / S z 1 Y r / U q u 9 2 8 7 a + S p 0 o 3 W T 8 w L d s f b R b 6 / u I g w = < / l a t e x i t > |ni $ |n + 3i < l a t e x i t s h a 1 _ b a s e 6 4 = " e N e M n Z N 5 M v M y U G x S G c s p / 2 a h l 2 8 = " > A A A D f 3 i c d V N L b x M x E H Y a H i W 8 U j h y s U h Q k a D R b i D Q S h w q c e F Y J N J W i q P K 6 5 1 N r P i x s p 2 m k e s / y Y 1 n : < l a t e x i t s h a 1 _ b a s e 6 4 = " N o F C 9 C y a P u F H 0 l a J P b m 9 A U 0 V u q M = " > A M u e J M g k Q A j W R p 9 3 8 x 8 n k d W C W 5 s k n x v 7 D X v 3 L 1 3 v / U g e v j o 8 Z P 9 9 s H T c 1 P O N Y M h K 0 W p L z N q Q H A F Q 8 u t g M t K A 5 W Z g I t s 9 q H m L 6 5 B G 1 6 q L 3 Z Z w V j S K 8 U L z q g N 0 K T 9 k a i S q x y U J R X V W b l w 6 R s m v S M W F t Y U j m S F A c 3 B Y I w 7 J O e m E n R p 7 F I A J p l 0 O D n y H e / x p B 0 n 3 d 7 J I H l 9 g v 9 2 0 m 6 y s h h t 7 G x y 0 P h J 8 p L N Z V B m g h o z S p P K j h 3 V l j M B P i o e / A K x c F R 0 = < / l a t e x i t > 0, < l a t e x i t s h a 1 _ b a s e 6 4 = " x o g M N / F 5 B R 5 M u e J M g k Q A j W R p 9 3 8 x 8 n k d W C W 5 s k n x v 7 D X v 3 L 1 3 v / U g e v j o 8 Z P 9 9 s H T c 1 P O N Y M h K 0 W p L z N q Q H A F Q 8 u t g M t K A 5 W Z g I t s 9 q H m L 6 5 B G 1 6 q L 3 Z Z w V j S K 8 U L z q g N 0 K T 9 k a i S q x y U J R X V W b l w 6 R s m v S M W F t Y U j m S F A c 3 B Y I w 7 J O e m E n R p 7 F I A J p l 0 O D 3 y H e / x p B 0 n 3 d 7 J I H l 9 g v 9 2 0 m 6 y s h h t 7 G x y 0 P h J 8 p L N Z V B m g h o z S p P K j h 3 V l j M B P i x F X u 9 s Y 8 X 2 r m y n T W T 5 P X g a r u U V e A h u i A s H 8 C Z F I g F G s j T 6 v p n 5 P D 9 5 L b i x S f K l F V 2 7 f u P m r Z 3 b 7 T t 3 7 9 3 f 7 e w 9 O D b V X D M Y s k p U + j S n B g R X M L T c C j i t N V C Z C z j J Z 2 8 b / u Q c t O G V + m C X N Y w l P V O 8 5 I z a A E 0 6 G V E V V w U o S 2 q q 8 2 r h 0 p d M e k c s L K w p H c l L A 5 q D w Z g U l T X e 4 0 k n T n r Z w S B 5 c Y D / d t J e s r I Y X d n R Z K / 1 I + S y u Q w q T F B j R m l S 2 7 G j 2 n I m w L f J 3 E B N 2 Y y e w S i 4 i k o w Y 7 d q z u M n A S l w W e n w l M U r 9 M 8 M R 6 U x S 5 m H S E n t 1 G x z D f h P L p c B z i t R 4 I b b o A r T a G 1 9 z J a D s e O q n l t Q b P 2 v c i 6 w r X A z W F x w D c y K Z X A o 0 z y 0 h t m U a s p s G H + b K L h g l Z R U F Y 6 c A 6 v 9 K B 2 v P B y n 5 P l T M j V A g X B s 3 R r S g t S O a c g N h z f 0 0 h Q U m 0 6 Z 8 9 0 3 I d U R U F 6 C 3 q a Z Y Y L 0 P 5 / L 7 J v D / n e O s l 7 7 q J e + z + H B w d T g 7 6 B F 6 j J 6 h F L 1 G h + g d O k J D x N B H 9 A l d o s / R Z f Q 1 + h Z 9 X 4 d G r a u c h 2 j D o p + / A K a q E b 0 = < / l a t e x i t > . . .
n : < l a t e x i t s h a 1 _ b a s e 6 4 = " N o F C 9 C y a P u F H 0 l a J P b m 9 A U 0 V u q M = " > A M u e J M g k Q A j W R p 9 3 8 x 8 n k d W C W 5 s k n x v 7 D X v 3 L 1 3 v / U g e v j o 8 Z P 9 9 s H T c 1 P O N Y M h K 0 W p L z N q Q H A F Q 8 u t g M t K A 5 W Z g I t s 9 q H m L 6 5 B G 1 6 q L 3 Z Z w V j S K 8 U L z q g N 0 K T 9 k a i S q x y U J R X V W b l w 6 R s m v S M W F t Y U j m S F A c 3 B Y I w 7 J O e m E n R p 7 F I A J p l 0 O D n y H e / x p B 0 n 3 d 7 J I H l 9 g v 9 2 0 m 6 y s h h t 7 G x y 0 P h J 8 p L N Z V B m g h o z S p P K j h 3 V l j M B P i o e / A K x c F R 0 = < / l a t e x i t > 0, < l a t e x i t s h a 1 _ b a s e 6 4 = " x o g M N / F 5 B R 5 M u e J M g k Q A j W R p 9 3 8 x 8 n k d W C W 5 s k n x v 7 D X v 3 L 1 3 v / U g e v j o 8 Z P 9 9 s H T c 1 P O N Y M h K 0 W p L z N q Q H A F Q 8 u t g M t K A 5 W Z g I t s 9 q H m L 6 5 B G 1 6 q L 3 Z Z w V j S K 8 U L z q g N 0 K T 9 k a i S q x y U J R X V W b l w 6 R s m v S M W F t Y U j m S F A c 3 B Y I w 7 J O e m E n R p 7 F I A J p l 0 O D 3 y H e / x p B 0 n 3 d 7 J I H l 9 g v 9 2 0 m 6 y s h h t 7 G x y 0 P h J 8 p L N Z V B m g h o z S p P K j h 3 V l j M B P i x F X u 9 s Y 8 X 2 r m y n T W T 5 P X g a r u U V e A h u i A s H 8 C Z F I g F G s j T 6 v p n 5 P D 9 5 L b i x S f K l F V 2 7 f u P m r Z 3 b 7 T t 3 7 9 3 f 7 e w 9 O D b V X D M Y s k p U + j S n B g R X M L T c C j i t N V C Z C z j J Z 2 8 b / u Q c t O G V + m C X N Y w l P V O 8 5 I z a A E 0 6 G V E V V w U o S 2 q q 8 2 r h 0 p d M e k c s L K w p H c l L A 5 q D w Z g U l T X e 4 0 k n T n r Z w S B 5 c Y D / d t J e s r I Y X d n R Z K / 1 I + S y u Q w q T F B j R m l S 2 7 G j 2 n I m w L f J 3 E B N 2 Y y e w S i 4 i k o w Y 7 d q z u M n A S l w W e n w l M U r 9 M 8 M R 6 U x S 5 m H S E n t 1 G x z D f h P L p c B z i t R 4 I b b o A r T a G 1 9 z J a D s e O q n l t Q b P 2 v c i 6 w r X A z W F x w D c y K Z X A o 0 z y 0 h t m U a s p s G H + b K L h g l Z R U F Y 6 c A 6 v 9 K B 2 v P B y n 5 P l T M j V A g X B s 3 R r S g t S O a c g N h z f 0 0 h Q U m 0 6 Z 8 9 0 3 I d U R U F 6 C 3 q a Z Y Y L 0 P 5 / L 7 J v D / n e O s l 7 7 q J e + z + H B w d T g 7 6 B F 6 j J 6 h F L 1 G h + g d O k J D x N B H 9 A l d o s / R Z f Q 1 + h Z 9 X 4 d G r a u c h 2 j D o p + / A K a q E b 0 = < / l a t e x i t > . . .
O 2 e 9 b v q m m 3 z u x S e D z e G 0 0 B P 0 F L 1 A K X q L T t A n d I q G i K G v 6 B u 6 Q T / 2 b p q N Z q s Z r U P 3 G p u c Q 7 R l z c N f e c 0 V a Q = = < / l a t e x i t > n : < l a t e x i t s h a 1 _ b a s e 6 4 = " N o F C 9 C y a P u F H 0 l a J P b m 9 A U 0 V u q M = " > A M u e J M g k Q A j W R p 9 3 8 x 8 n k d W C W 5 s k n x v 7 D X v 3 L 1 3 v / U g e v j o 8 Z P 9 9 s H T c 1 P O N Y M h K 0 W p L z N q Q H A F Q 8 u t g M t K A 5 W Z g I t s 9 q H m L 6 5 B G 1 6 q L 3 Z Z w V j S K 8 U L z q g N 0 K T 9 k a i S q x y U J R X V W b l w 6 R s m v S M W F t Y U j m S F A c 3 B Y I w 7 J O e m E n R p 7 F I A J p l 0 O D n y H e / x p B 0 n 3 d 7 J I H l 9 g v 9 2 0 m 6 y s h h t 7 G x y 0 P h J 8 p L N Z V B m g h o z S p P K j h 3 V l j M B P i J z A x V l M 3 o F o + A q K s G M 3 a p h j 1 8 E J M d F q c N T F q / o e / A K x c F R 0 = < / l a t e x i t > 0, < l a t e x i t s h a 1 _ b a s e 6 4 = " x o g M N / F 5 B R 5 M u e J M g k Q A j W R p 9 3 8 x 8 n k d W C W 5 s k n x v 7 D X v 3 L 1 3 v / U g e v j o 8 Z P 9 9 s H T c 1 P O N Y M h K 0 W p L z N q Q H A F Q 8 u t g M t K A 5 W Z g I t s 9 q H m L 6 5 B G 1 6 q L 3 Z Z w V j S K 8 U L z q g N 0 K T 9 k a i S q x y U J R X V W b l w 6 R s m v S M W F t Y U j m S F A c 3 B Y I w 7 J O e m E n R p 7 F I A J p l 0 O D 3 y H e / x p B 0 n 3 d 7 J I H l 9 g v 9 2 0 m 6 y s h h t 7 G x y 0 P h J 8 p L N Z V B m g h o z S p P K j h 3 V l j M B P i J z A x V l M 3 o F o + A q K s G M 3 a p h j 1 8 E J M d F q c N T F q / l a t e x i t > 1, < l a t e x i t s h a 1 _ b a s e 6 4 = " B 1 O Q t p T A T q k d 1 B x Q 8 V 4 P t a w h H j s = " > A x F X u 9 s Y 8 X 2 r m y n T W T 5 P X g a r u U V e A h u i A s H 8 C Z F I g F G s j T 6 v p n 5 P D 9 5 L b i x S f K l F V 2 7 f u P m r Z 3 b 7 T t 3 7 9 3 f 7 e w 9 O D b V X D M Y s k p U + j S n B g R X M L T c C j i t N V C Z C z j J Z 2 8 b / u Q c t O G V + m C X N Y w l P V O 8 5 I z a A E 0 6 G V E V V w U o S 2 q q 8 2 r h 0 p d M e k c s L K w p H c l L A 5 q D w Z g U l T X e 4 0 k n T n r Z w S B 5 c Y D / d t J e s r I Y X d n R Z K / 1 I + S y u Q w q T F B j R m l S 2 7 G j 2 n I m w L f J 3 E B N 2 Y y e w S i 4 i k o w Y 7 d q z u M n A S l w W e n w l M U r 9 M 8 M R 6 U x S 5 m H S E n t 1 G x z D f h P L p c B z i t R 4 I b b o A r T a G 1 9 z J a D s e O q n l t Q b P 2 v c i 6 w r X A z W F x w D c y K Z X A o 0 z y 0 h t m U a s p s G H + b K L h g l Z R U F Y 6 c A 6 v 9 K B 2 v P B y n 5 P l T M j V B C N x + v 7 Y + C B d Q k s K E Y U 1 d s x j L q S A z 0 G q / 1 w d J V o T f L B o a F a I M e n 6 U h c q y 2 e U a x F 1 H C m 5 q Q Z f G L g V g 0 o S 5 O P U u z r z v B r 2 1 o O F S Q E g N 4 9 A g X B s 3 R r S g t S O a c g N h z f 0 0 h Q U m 0 6 Z 8 9 0 3 I d U R U F 6 C 3 q a Z Y Y L 0 P 5 / L 7 J v D / n e O s l 7 7 q J e + z + H B w d T g 7 6 B F 6 j J 6 h F L 1 G h + g d O k J D x N B H 9 A l d o s / R Z f Q 1 + h Z 9 X 4 d G r a u c h 2 j D o p + / A K a q E b 0 = < / l a t e x i t > . . .
have been traced out.In the limit γ → 0 + we recover the usual ε-prescription of the retarded/advanced propagators, and effectively consider a system in infinitesimal contact with an external heat bath, described by the canonical density operator ρ = e −β Ĥ /Z.
Applying the same collisional broadening to the spectral function of the anharmonic oscillator in (2.5), one analogously obtains with ∆E mn = E m − E n .The only assumption here is that the width γ is not affected by the anharmonicity, in particular, that it does not acquire a frequency dependence for λ = 0. 1 Otherwise this is an exact expression which we will use for our benchmark calculations.
To compute the spectral function of the anharmonic oscillator with this constant broadening we discretize the coordinate x on a sufficiently large interval and solve the eigenvalue problem for the corresponding finite Hamiltonian matrix, obtained from (2.1) in coordinate space, by exact diagonalization.We then verify that the interval is large enough to cover the support of all relevant eigenfunctions at a given temperature and that the discretization is fine enough to obtain precise results sufficiently 1 Although this is reasonable for small γ in the Ohmic bath, it will not hold in more realistic situations with ultraviolet cutoff ω D as in the Drude model for the bath, where memory effects will necessarily occur on time scales shorter than ω −1 D , inducing a frequency dependent damping γ(ω) on scales ω ∼ ω D [19].
high up in the spectrum.For the parameters λ, γ and T used in Sec.VI typical interval sizes are x ∈ [−20 , 20] (in units of 1/ √ ω 0 ) with ∼ 3000 grid points.An example is shown in Figure 1 where we plot the spectral function of the quartic anharmonic oscillator at a rather large coupling of λ/ω 3 0 = 4 for a temperature T /ω 0 = 1, and with a comparatively small damping of γ/ω 0 = 0.03, from Eq. (2.9).In particular, the small damping allows resolving the individual transitions: (a) between adjacent energy levels |n and |n + 1 in the main peak which are split up because they are no-longer equidistant when λ = 0, (b) transitions across three levels between |n and |n + 3 in the second sequence of peaks at higher frequencies which are also due to the sizable λ > 0, and (c) analogous transitions across five levels |n and |n + 5 , here for frequencies above 8ω 0 .The finite temperature manifests itself in the contributions from the individual transitions with n ≥ 1 in each sequence which vanish for T → 0 where only the corresponding ground-state transitions |0 ↔ |1 , |3 , |5 , . . ., survive.In the results section below we will implicitly assume all dimensionful quantities to be quoted in the appropriate units of ω 0 .

III. CLASSICAL SPECTRAL FUNCTIONS
In the classical-statistical limit [20][21][22][23][24], the full Heisenberg equations of motion are truncated by approximating all quantum mechanical expectation values of products of operators by the corresponding products of their expectation values, i.e.Ô1 . . .Ôn → Ô1 . . .Ôn , which can be formalized via the real-time path integral formulation for classical-statistical systems [17,44,45].One then arrives at the well-known Langevin equations of motion describing the purely classical dissipative dynamics [17,18], where X(t) and P (t) are the expectation values of x(t) and p(t) in coherent states, and the stochastic fluctuating force ξ is given by a Gaussian white noise with zero mean and variance 2γT , As mentioned in the previous section, the spectral function ρ can be defined as the thermal expectation value of the commutator of two operators, cf.(2.2) for our oscillator here.It is related to the corresponding time-ordered Green's function G T via where F (t, t ) is the statistical two-point function defined as the expectation value of the corresponding anti-commutator [42,43], i.e. here, Since these equilibrium two-point functions all depend only on the time difference t−t , their Fourier transforms F (ω), ρ(ω) depend on a single frequency ω, where for ρ(ω) we use the definition in Eq. (2.3), but for the even function F (−ω) = F (ω) the conventional form so that In thermal equilibrium one furthermore applies the periodicity or Kubo-Martin-Schwinger (KMS) condition, in the decomposition of Eq. (3.5), in order to derive the fluctuation-dissipation relation (FDR), e.g.see [46], where we have used the special convention for the definition of ρ(ω) in (2.3), and n B (ω) = 1/(e βω −1).In the classical limit T ω we approximate coth(ω/2T ) ≈ 2T /ω.The classical FDR from (3.9) then relates the corresponding classical two-point functions, In the time domain, undoing the Fourier transform, this reads Furthermore, because the statistical two-point function is in the classical limit given by the purely thermal correlation function the spectral function (3.11) can be written as [23] where P = Ẋ is the conjugate momentum which has zero mean in the thermal ensemble, P (t) β = 0.Because of time-reversal invariance of the thermal expectation values, the two terms in (3.13) are the same, and the explicit anti-symmetrization in this definition of ρ c (−t) = −ρ c (t) can be introduced without loss.Evaluating Eq. (3.13) provides a straightforward way of calculating the spectral function in the classical-statistical limit [20][21][22][23][24].
IV. GAUSSIAN-STATE APPROXIMATION A. Closed system Before considering the coupling to an environment, we first briefly discuss the Gaussian state approximation (GSA) for a closed system.The GSA is obtained by truncating the full Heisenberg equations of motion for the canonically conjugate Heisenberg operators x(t) and p(t): The equations of motion for the expectation values can be obtained by averaging equations (4.2a) and (4.2b) over some density operator ρ describing the mixed initial state of the ensemble.These equations then contain expectation values of the form x2 (t) and x3 (t) , whose evolution equations in turn include expectation values of even higher-order combinations of x and p.This leads to an infinite hierarchy of equations that cannot be solved analytically or numerically without further approximations.Moreover, to deal with expectation values of products of x and p we follow Ref. [26] and introduce the Wigner transform of the density matrix in position eigenstates, w(x, p) = dy e −ipy x + y/2|ρ|x − y/2 , ( which allows expressing the expectation values of symmetrized products of x and p in the form of classical phase space integrals, such as e.g. To truncate the infinite set of equations given by (4.2a) and (4.2b), the density matrix is itself approximated by a Gaussian, and can therefore be characterized by a Gaussian Wigner function likewise [25], Here, the parameters X ≡ x , P ≡ p describe the center of the Gaussian wave packet in coordinate and momentum space.As such they are not necessarily the expectation values in coherent states yet, here.The symmetrized connected expectation values characterize the dispersions of the wave packet, and N is a normalization factor.Equations (4.2a) and (4.2b) are then averaged over the Gaussian state with the Wigner function (4.5).Applying Wick's theorem as needed, one then obtains d dt X = P, (4.6a) To evolve the dispersions σ xx , σ xp , and σ pp , the Heisenberg equations for the corresponding symmetrized operator products are employed Averaging equations (4.7a) -(4.7c) over the Gaussian state with the Wigner function (4.5), applying Wick's theorem, and subtracting the disconnected contributions, the remaining equations of motion are obtained as where we have introduced the curvature of the potential The equations of motion (4.6a), (4.6b) and (4.8a) -(4.8c) can be integrated numerically using a symplectic leapfrog algorithm, as described in Appendix A 3, to obtain a complete description of the Gaussian state at any time.
We conclude this section with some general remarks on the formal structure of the GSA which will be particularly relevant for the systematic construction of a thermal equilibrium state in Section IV B 2 below.
• In general, we call a (possibly mixed) state ρ Gaussian, if its Wigner transform (4.3) has the form of a Gaussian probability distribution (4.5) for some X, P, σ xx , σ xp , σ pp , with obvious generalization to an arbitrary number of degrees of freedom, where the σ's are replaced by the covariance matrix Σ.
• Note that a Gaussian state ρG defined in this way is not necessarily a pure state.This can be seen most easily by calculating the von Neumann entropy [26].
One observes that the purity of a Gaussian state is related to the determinant det Σ of the covariance matrix Σ.This determinant is a product of pairs of symplectic eigenvalues f k , one per bosonic degree of freedom (dof).With N dof of them, Heisenberg's uncertainty relation then implies that det Σ = On the other hand, the von Neumann entropy vanishes and the Gaussian state is pure, if and only if f k = 1/2 for all dof's.For a single degree of freedom as above, for example, we have and restricting to pure Gaussian states therefore defines a non-linear subset G ⊂ H = L 2 (R) of the full Hilbert space (of square-integrable functions) which can be parametrized by a 4-dimensional manifold with X, P ∈ (−∞, ∞), σ xx , σ pp ∈ (0, ∞) and constrained by σ xx σ pp ≥ 1/4.To ensure that the von Neumann entropy vanishes, the offdiagonal variance is then fixed up to a sign by • To describe more general Gaussian states ρG , we again follow Ref. [26] and define the set of mixed Gaussian states in terms of those density operators that can be written as mixtures of the pure Gaussian states in G, with probabilities p(ψ) that are Gaussian likewise.

B. Caldeira-Leggett model
In order to study the dynamical properties of thermal equilibrium states, we introduce a coupling between the system, our anharmonic oscillator we will also refer to as the particle, and the environment consisting of an ensemble of harmonic oscillators, which models a fixedtemperature heat bath.Such a model in the canonical operator formalism as well as in the functional path integral formulation (after Feynman and Vernon) has been discussed frequently in the literature, see for example Refs.[18,19,[47][48][49][50][51][52][53][54][55].Often Born and Markov approximations are employed, leading to master equations which are easy to solve but not generally accurate.Exact solutions also have been obtained analytically [55].However, this is unfortunately not the case for the anharmonic oscillator.
The total Hamiltonian under consideration then consists of those of the system S, the heat bath B together with their interactions I, and reads [18,19] where φs , πs denote the coordinate and the conjugate momentum of the heat-bath oscillator with index s, ω s is its eigenfrequency, and g s the coupling constant of its linear coupling to the coordinate x.The quadratic term in ĤI serves to exactly compensate the bath-induced (negative) shift of the oscillator frequency squared, that would otherwise arise.This guarantees that ω 2 0 is the physically measured natural frequency of the noninteracting system oscillator with damping.Completing the square, we can then write the interaction-plus-bath part as

Heisenberg-Langevin equations of motion
Introducing a spectral function to describe the ensemble of bath modes by which corresponds to a positive-definite spectral density, one can eliminate the heat bath from the full set of Heisenberg equations to derive the quantum equations of motion of the Caldeira-Leggett model for a general heat bath described by J(ω), see e.g.Refs.[18,19], with an operator-valued fluctuating force Here η(t) is defined such that it only acts in the bath's Hilbert space, and we have furthermore introduced the damping kernel For the simplest case of an Ohmic bath with damping constant γ and a sharp cutoff at ω = Λ, we have and thus Note that the 'transient term' −γ(t)x(0) in Eq. (4.16) corresponds to a sudden initial shift in the thermal distribution when the particle is connected to the bath [18] at t = 0.It can therefore safely be omitted in all calculations as long as we are interested in times Λt 1, or analogously for ω D t 1 with γ(t) = γω D exp{−ω D t} in the Drude model for the bath [19].
Equations (4.15a), (4.15b) together are commonly referred to as the Heisenberg-Langevin equations (HLEs) in the literature [55,56].They generalize the classical Langevin equations (3.1), (3.2) which describe the dissipative dynamics of the expectation values of position and momentum.The corresponding HLEs, on the other hand, encode the full dynamics of the quantummechanical operators in the Heisenberg picture [56].A common approximation is to replace the operators in the HLEs by their expectation values and the quantum noise by a classical colored-noise source, which results in the so-called quasiclassical Langevin equations [18].While these can provide a reasonable description of nearly harmonic systems [57,58], other important general features of the HLEs such as the Heisenberg uncertainty principle for the operators x and p are lost.Finally, the classical Langevin equations are obtained in the Markovian limit in which all memory effects are disregarded and the noise becomes local in time (white noise limit).
As a first step towards a solution of the Heisenberg-Langevin equations in the GSA, we assume that at the initial time t = 0 the bath is in equilibrium, while the system particle is described by some density matrix ρS with Gaussian Wigner function as defined in Eq. (4.5), i.e. we assume that the initial state is given by ρ0 ≡ ρ(t 0 ) = ρS ⊗ ρB . (4.20) We defer the precise specification of the equilibrium density matrix ρB of the heat bath to Section IV B 2 where we discuss some further subtleties associated with its GSA description.
Our goal here is to formulate a Langevin-type equation in the Gaussian state formalism, which preserves more of the features of the HLE than the quasiclassical or classical approach, thus having extended applicability.First of all, we consider the most general Gaussian Wigner function W that describes the entire system of oscillator and heat bath, (4.21) where the vector ζ = (x, p, ..., ϕ s , π s , ...) ∈ Γ denotes a point in the full phase space Γ of the system, Z(t) = (X(t), P (t), ..., Φ s (t), Π s (t), ...) are the expectation values of the corresponding Heisenberg operators, and represents the corresponding covariance matrix.Note that Σ will in general also contain cross-correlations such as σ xϕs between the oscillator particle and the heat bath, which encode quantum entanglement.In order to translate (4.15a), (4.15b) into corresponding equations of motion within the GSA [25], in a systematic and thermodynamically consistent way, we adopt the following procedure: (i) Average the Heisenberg equations of motion (4.15a), (4.15b) together with (4.16) using the general Gaussian state (4.21) and write down the resulting equations for all dynamic quantities contained in Z and Σ. Evaluate correlation functions using Wick's theorem.
(ii) Integrate out the heat-bath degrees of freedom by solving the equations of motion obtained from step (i) for the bath oscillator coordinates Φ s , Π s and the system-bath cross-correlations σ xϕs , σ pϕs , σ xπs , σ pπs , under the assumption that the full solution X(t), P (t), σ xx (t), σ xp (t), σ pp (t) is already known.Use initial conditions of the form , to obtain the 2-point correlators, but leave the initial phase-space variables Φ s (0), Π s (0) of the heat bath arbitrary.Insert these formal solutions into the five remaining equations of motion for the Gaussian particle.
(iii) Notice that the remaining five equations of motion for the particle only depend on the initial conditions of the bath, Φ s (0), Π s (0), Σ bath (0), where Σ bath (0) (the lower right corner of Σ(0) in (4.23)) describes all the 2-point correlators of the bath oscillators' phase-space variables.Therefore, to introduce thermal fluctuations, all we need to do is distribute the initial expectation values of the bath oscillators according to a thermal distribution (to be specified below), e.g.via a fluctuating force term.
Since the quantum average over the Gaussian Wigner function (4.21) in step (i) according to the prescription in (4.4) is entirely different from the thermal average over the initial conditions in step (iii), we denote the former by • • • and the latter by Before we can investigate the solution of the equations of motion from steps (i) and (ii), we first have to take a closer look at how to represent the 'thermal initial state of the bath' which enters in step (iii) in the Gaussian approximation, i.e. how are Φ s (0), Π s (0) and Σ bath (0) distributed in a thermal-equilibrium state at a given temperature T .This is surprisingly subtle, however, as discussed in the next subsection.

Gaussian thermal equilibrium state
Since the bath is described as an ensemble of harmonic oscillators in Gaussian mixed states, we still have to specify what precisely we mean by a 'thermal equilibrium state' for a single bath oscillator in the harmonic case with λ = 0, where the Gaussian approximation becomes exact.Constructing a thermal ensemble of Gaussian states that models a quantum canonical state at temperature T , along the lines described for mixed Gaussian states at the end of Section IV, is not entirely trivial.We have to express the full quantum-mechanical mixed thermal state of a harmonic oscillator, described by the density operator (with spectral representation in the energyeigenstates |n of the harmonic oscillator), that acts on the full Hilbert space H, in terms of a ρG acting on G ⊂ H to describe a Gaussian mixed state, by a density operator of the form (4.10).In general, a mixed thermal state describing a canonical ensemble is not of this form, so this involves an approximation.For the harmonic oscillator, however, it can be easily verified from the definition in (4.3) that the density operator ρHO in Eq. (4.24) does have a Gaussian Wigner function, which is given by where is its thermal distribution function.To express the canonical equilibrium ensemble at temperature T , represented by the mixed-state density operator ρHO in Eq. (4.24), as a mixed Gaussian state of the form (4.10), we now use our classical phase-space variables X = x α and P = p α , here restricted again to the expectation values of x and p in coherent states These states are characterized by the complex variable α whose real and imaginary parts are given by X and P , We therefore denote these coherent states simply by the two real phase-space variables in the following, i.e.
The coherent states of the harmonic oscillator have minimal uncertainty with On the other hand, for the mixed equilibrium state in the canonical ensemble from (4.24) at temperature T = 1/β, with Ô β = Tr Ô ρHO , one has x β = 0 and p β = 0, and the full variances σ (β) in the thermal state are readily computed (or read off from (4.25)) as, where n B (ω 0 ) = 1/(exp(βω 0 ) − 1) is the Bose-Einstein distribution.We therefore see here explicitly that the full thermal widths σ xx and σ pp of the oscillator's phasespace variables can be split into purely thermal or 'classical' parts σ c xx , σ c pp plus the purely 'quantum' parts σ 0 xx , σ 0 pp from minimal uncertainty, as noted in Ref. [26], i.e. with The minimal-uncertainty variances are already included in each coherent pure state.To define a mixed Gaussian state ρG with the variances of the thermal equilibrium ensemble, we therefore only include the classical thermal widths of (4.32) in the incoherent sum [26] One may directly verify that such a thermal state is indeed a stationary solution of the harmonic Gaussian equations of motion (4.6a), (4.6b) and (4.8a) -(4.8c).In general, it maps Gaussian states to Gaussian states in the Hilbert space, as required.This finishes the discussion of the Gaussian thermal state of a single harmonic oscillator, and we now continue to model the entire system consisting of our anharmonic oscillator coupled to an ensemble of harmonic oscillators, in mixed Gaussian states with thermal variances as described here, reintroducing the heat-bath index s.

Heisenberg-Langevin equations in the GSA
We now turn to the form of the Heisenberg-Langevin equations of motion (4.15a), (4.15b) evaluated in the Gaussian approximation.Steps (i) and (ii) from the quantum equations of motion (4.15a), (4.15b) with the Ohmic bath (4.18), on time scales Λt 1 for a sufficiently large cutoff Λ, cf.(4.19), lead to for the center (X, P ).The derivation is the same as the one for the classical Langevin equations of motion, except for the application of Wick's theorem to the 3-point correlator x3 (t) .The fluctuating force term ξ(t) is given by the expectation value of the quantum stochastic force ξ(t) from (4.16), where the transient initial shift vanishes if the expectation value of x(0) does, It thus only depends on the initial conditions of the bath oscillators' phase-space variables Φ s (0), Π s (0).Although no-longer operator valued and hence classical, this GSA noise ξ(t) is colored in general, however, as we will discuss in Subsection IV B 5 below.
For the Gaussian widths, the analogous averaging of step (i) leads to Here we have already assumed the Ohmic heat bath J Λ (ω) in the limit Λ → ∞ where the memory integrals collapse, cf.(4.19).The initial delta-functions can then safely be neglected.
On the other hand, the irreducible correlators of the oscillator particle's position and momentum with the fluctuating force operator, x(t)η(t) and p(t)η(t) , both contain a logarithmically divergent contribution d(Λ) ∼ ln(Λ/C 0 ) where C 0 is the t = 0 initial value of the timedependent curvature of the potential, We will show in Appendix A 1 that this divergence can be absorbed by a formally infinite but time-independent shift of the particle's momentum width σ pp → σ pp − d(Λ) in such a way that it cancels from both equations (4.36b) and (4.36c) which then describe the time dependence of the ultraviolet-finite part of σ pp , together with finite σ xp and σ xx at all times.Note that the divergence of σ pp with Λ → ∞ is an unavoidable effect of the unrealistic assumption of an Ohmic bath without ultraviolet (UV) cutoff.It can be interpreted as corresponding to the bath continuously 'measuring' the position of the particle with arbitrarily high 'resolution' without UV cutoff for Λ → ∞ [18,49].Finally, for the evaluation of the irreducible correlators x(t)η(t) and p(t)η(t) according to rules (ii) and (iii), we need to make an additional adiabatic approximation as explained explicitly also in Appendix A 1. In this adiabatic approximation we assume that we can average the curvature C of the potential in (4.37) over time scales that are large compared to the relaxation time of the heat bath.The heat-bath oscillators are then considered as the fast degrees of freedom that can adjust to slow changes in the curvature C(t).This adiabatic approximation then yields for the Gaussian widths, where for the Ohmic bath J Λ (ω) with Λ → ∞ the fluctuating force x(t)η(t) after ultraviolet subtraction yields CF − ∆ K .The first contribution is obtained from upon inserting C(t) together with an equally slowly varying frequency assuming weak damping.The second contribution to the fluctuating force is the ultraviolet subtracted one, given by which determines the relevant (ultraviolet-finite) part of the fluctuating force p(t)η(t) on the ultraviolet subtracted σ pp as well, and which vanishes when C ≡ C 0 is used at all times in the static limit, see Appendix A 1. Together with Eqs.(4.34a) and (4.34b), these equations for the evolution of the Gaussian widths, from Eqs. (4.38a) -(4.38c), constitute the full set of equations of motion for our particle in the anharmonic potential, and in contact with an external heat bath, i.e. the Heisenberg-Langevin equations within the GSA in our adiabatic approximation.
The static approximation is obtained from Eqs. (4.38a) -(4.38c) by simply using the time-independent C = C 0 , for which we have ∆ K (C 0 ) = 0.In this case, the equations for the widths, cf.(A.26a) -(A.26c), can be solved independently of those for the coordinates X, P , see Appendix A 1. The asymptotic behavior of the solution, cf.Eq. (A.28), uniquely fixes Therefore, with the static solution, there are only two effects remaining of the GSA in comparison with the classical time-evolution.These are (a) a time-dependent shift of the oscillator frequency in Eq. (4.34b), with the stationary value (4.42) of σ xx (t) for sufficiently late times, and (b) a modified colored noise ξ(t), which we will specify in Subsection IV B 5 below.At very high temperatures, the frequency shift is negligible, and the noise becomes white again, such that the correct classical limit is guaranteed to be recovered in this static approximation.
In the other direction, to go beyond the adiabatic approximation, one could in principle include the feedback of the time dependence of the curvature C(t) in (4.37) on the off-diagonal variances between system particle and heat-bath oscillators via post-adiabatic corrections in the spirit of time-dependent perturbation theory, in the future, as briefly outlined in Appendix A 1 as well.For the results presented below, we have either used the static approximation with constant C 0 , cf.Eqs.(A.26a) -(A.26c), or the adiabatic approximation in Eqs.(4.38a) -(4.38c), for comparison.

Initial Conditions
In the adiabatic approximation the thermal equilibrium value C 0 is an inital condition that has to be know beforehand.It can be obtained by maximizing the von Neumann entropy at fixed energy.
For a mixed Gaussian state of our single bosonic degree of freedom the von Neumann entropy S = −Tr(ρ ln ρ) can be written in terms of the symplectic eigenvalue which is related to the pair of eigenvalues λ ± = ±if of ΣΩ, where Ω is the symplectic matrix for the canonically conjugate variables x and p. Thermal equilibrium implies σ xp = 0 and we thus have The von Neumann entropy can then be written as Because we also have X = 0 and P = 0 in thermal equilibrium, the expectation value of energy in the Gaussian state reduces to We can thus express σ pp in terms of E and σ xx and write Because the entropy (4.45) increases monotonically with f , it reaches its maximum when f does, which is the case when ∂f 2 /∂σ xx = 0.This yields The temperature T is now introduced using .
Working out the partial derivatives with respect to f and σ xx , we thus obtain with C 0 = ω 2 0 + λ 2 σ xx in the interacting case for λ = 0. Before we start a simulation at a given the temperature, we can therefore calculate C 0 numerically via (4.50).In the static approximation this is then fixed, and so is σ xx in the HLEs (4.34a) and (4.34b) for X and P .
The underlying initial conditions for the widths σ ϕsϕs and σ πsπs of the heat-bath oscillators in Eq. (4.23) correspond to the thermal harmonic-oscillator variances Eqs.(4.30a) -(4.30c), with additional off-diagonal couplings σ xϕs between system and bath suddenly switched on at t = 0 as explained in more detail in Appendix A 1 b.
One crucial point left to mention here, however, is that beyond the static approximation, the widths σ xx , σ pp , σ xp actually do evolve non-trivially in time, even in the adiabatic approximation, when C(t) is assumed to vary slowly in time.This is because the relaxation time for the widths of the system particle to approach their stationary limits is given by 1/γ, and this relaxation time is in general not negligible compared to the characteristic time scale δt of the variations δC(t).Assuming, in the adiabatic approximation, that the heat-bath dof's are fast compared to this characteristic time δt is totally different from assuming that 1/γ is.In fact, for small damping γ we expect to have 1/γ δt 2π/ω s for the relevant high frequencies that dominate the Ohmic bath.We will further comment on this in Section IV B 6 below, after elaborating on the colored noise needed in either case.

Colored Noise
Step (iii) in our approach to modelling the heat bath in the GSA by the quantum mechanical expectation value ξ(t) ≡ ξ(t) of the stochastic quantum force from Eq. (4.35) requires specifying initial conditions for the thermal correlations of the bath oscillator expectation values Φ s = φs α and Π s = πs α in coherent states.In particular, the discussion leading to (4.33) implies that their initial thermal variances are given by the classical variances (4.32) for each heat-bath oscillator, Eq. (4.35) then yields With the definition of the spectral density of the heat bath in Eq. (4.14) we can thus finally represent the unequal-time correlations of ξ(t) in the form, for an arbitrary spectral distribution J(ω) of oscillators in the bath.For the Ohmic bath (4.18) in the limit Λ → ∞ the integral can be solved analytically, yielding in the time domain, where the terms in brackets approach 2δ(t − t ) for T → ∞, i.e. the classical-statistical limit with a noise term ξ(t) in the equations of motion which is Gaussian and local in time [23].For numerical purposes and completeness, in frequency domain the noise in (4.54) corresponds to where the white-noise limit is recovered from the classical Rayleigh-Jeans law with n B (ω) → T /ω at high temperatures.Therefore, ξ(t) represents a 'colored' noise term with a Gaussian autocorrelation, but with the quantum contribution subtracted, in the sense that ξ(t) vanishes identically at T = 0.This is due to the fact that in our framework, zero-point fluctuations are already naturally taken into account by the Gaussian widths, σ 0 xx , σ 0 xp , σ 0 pp and therefore do not contribute to the fluctuation of the mean coordinate and momentum.

Extracting the Gaussian spectral function
In order to extend the definition of the (antisymmetrized) classical-statistical spectral function in Eq. (3.13) to the correct quantum spectral function that respects the fluctuation-dissipation relation with the colored-noise distribution of the heat bath in the GSA, first note that the FDR (3.9) must be replaced by which allows us to use the classical-statistical extraction scheme for the spectral function, i.e. using Eq.(3.13) in the GSA as well, and just rescale the result to obtain the corresponding quantum spectral function.
Having defined the extraction scheme of the spectral function in the GSA, we can further elaborate on the problems that occur when we keep the full time evolution of the Gaussian widths (4.38a) -(4.38c).Since the GSA is only an approximation to the infinite hierarchy of the time evolution of higher moments of the (in principle exact) Wigner quasi-probability distribution w(x, p) in phase space, to quadratic order [25], one can no-longer guarantee that classical and quantum dynamics are strictly divided into the evolution of the expectation values X, P , and the second order moments σ xx , σ pp , σ xp , respectively.Therefore some classical contributions to the spectral function are also contained in the time evolution of the widths.Extracting the spectral function naïvely as in (4.57) is therefore not sufficient, if the quantum corrections are highly non-Gaussian by themselves, when the full time evolution of the widths is included.In this case, one would actually need some improved procedure to correctly extract these non-Gaussian contributions contained in the second-order unequal-time correlators.This undesirable effect is explicitly demonstrated for sufficiently large anharmonicity λ in Section VI.
In summary, a complete simulation eventually comprises the following steps: (i) Generate a random realization of the stochastic force ξ(t) distributed according to the colored noise correlations in (4.54).
(iii) Calculate the spectral function from the particular time history of the expectation values via (4.57).
(iv) Finally, repeat steps 1 through 3 and average the spectral function over all realizations of the stochastic force to obtain the thermal equilibrium spectral function.
More details of steps (i) and (ii) are described in Appendices A 2 and A 3, respectively.

V. REAL-TIME FRG
As another possibility for real-time calculations we have also explored the functional renormalization group (FRG) [27][28][29] on the closed time path [38].In the FRG one aims to compute the effective action Γ, which is obtained from the generating functional Z of the theory [3].Assuming the so-called effective average action Γ Λ is known at some initial energy-scale Λ in the ultraviolet (UV), the essential idea is to construct the full Γ step-by-step by 'interpolating' Γ k from the microscopic UV-action Γ Λ to the macroscopic action Γ in the infrared (IR).This is achieved by introducing a parameter k corresponding to the energy scale down to which the theory is valid.Using an auxiliary device called the regulator R k , one suppresses both thermal and quantum fluctuations of modes with ω < k.In particular, assuming that at sufficiently high energies (and momenta, where S ) the theory behaves classically, one may start the interpolation in the UV with Γ Λ = S, the classical bare action, see Eq. (5.27) below.(This can be shown rigorously by a saddle-point approximation, where the regulator term acts as a δ-functional in the limit k → ∞ [36].)One then computes the flow of the effective average action through 'theory space' [29] until it reaches the full macroscopic effective action Γ k→0 = Γ of the theory.Ideally, one would solve the flow from k → ∞ to k = 0, but in almost all practical applications it is sufficient to start with some large but finite UV cutoff k = Λ, and to stop at some sufficiently small finite value k = k IR in the IR.

A. Flow Equation
The FRG flow is determined by the equation that was derived by Wetterich [27,28] in the imaginary-time formalism, but here formulated on the closed time path (CTP) [38], for a real scalar field φ T (x) = (φ c (x), φ q (x)) in D = d + 1 spacetime dimensions and Keldysh space, where we adopt the convention of Ref. [37] and define the Keldysh rotation to be measure-preserving, i.e.
for the classical (or average) field φ c and the quantum (or response) field φ q , and vice-versa for the inverse transformation.By φ ± we denote the field components that live on the forward (+) and backward (−) parts of the closed time path.This convention has the convenient property that for a symmetric potential V (φ) the potential term −V (φ + ) + V (φ − ) on the closed time path is invariant under the interchange φ c ↔ φ q , i.e. does not distinguish between classical and quantum field components.For a detailed introduction to the formalism and a derivation of the flow equation see for example Refs.[38] and [36].
The flow in Eq. (5.1) on the right hand side includes the regulator R k and Γ (2) where we have already use a notation defined below, in Eq. (5.5), for brevity.R k and Γ k both have the form of self-energies on the closed time path, i.e. 2 × 2 matrices in Keldysh (c, q) space.Therefore, • denotes 2×2 matrix multiplication, and the trace also implies integration over adjacent coordinates.We denote functional derivatives of the effective average action as where Greek indices from the beginning of the alphabet denote CTP indices, α 1 , ..., α n ∈ {c, q}.Correspondingly, Γ k [φ] denotes the tensor of rank n containing all functional derivatives w.r.t. the classical and quantum fields.The full scale-dependent propagator G k in front of the background field expectation value φ is given by in compact matrix notation, or explicitly at the scale-dependent minimum φ 0,k , which satisfies the quantum equations of motion δΓ[φ 0,k ] = 0.The superscripts A, R, K and K denote the advanced, retarded, Keldysh and anomalous components, respectively.For the expressions in an arbitrary background field configuration φ(x) see for example Ref. [38].
Being the imaginary part of the retarded propagator, the spectral function ρ k (ω) can be computed from the retarded 2-point function in the usual way, (5.8)

B. Causal regulators
The need of respecting causality in the process of constructing regulators for the real-time FRG was already mentioned in Ref. [40], where the advantages of using such a causal regulator were manifest in the results for dynamical critical exponents.In this section, we construct such a causal regulator step-by-step in the 0 + 1 dimensional case.The construction is based on considering the regulator term as an additional self-energy with the causal matrix structure of the Keldysh action and can be generalized to field theories in higher dimensions.For our single real degree of freedom φ T (t) = (φ c (t), φ q (t)), we first add a term to the Keldysh action of the form, with being a 2 × 2-matrix on the CTP.We call a regulator causal if it maintains the causal structure of the corresponding Keldysh action, as e.g.defined in Chapter 2.7 of Ref. [17].Most importantly, a causal regulator has to ensure that the retarded (advanced) propagator has to stay retarded (advanced) during the flow.
Here, we propose a new quite general construction of causal regulators which proceeds as follows: We assume that the bilinear term ∆S k [φ] really is the result of a coupling of the field φ to an ensemble of Gaussian degrees of freedom which have been integrated out, e.g.see Chapter 3.2 of Ref. [17].In fact, the other way round, any term quadratic in the fields can be linearized via Hubbard-Stratonovich transformation to replace it by a linear coupling to the Gaussian Hubbard fields representing the ensemble.Therefore, the assumption that ∆S k [φ] is the result of integrating an ensemble of Gaussian degrees of freedom should still be fairly general.As such, after Fourier transform, the corresponding retarded/advanced components can readily be written as spectral integrals so that the spectral density J k (ω) ≥ 0 for ω > 0 of the fictitious Gaussian ensemble is given by the imaginary part of the retarded bath propagator D R bath (ω) which is here represented, after the Gaussian integration of the bath, by as can explicitly be checked from (5.11).Note that for a single oscillator with frequency ω k and coupling g k in the bath, it is here normalized to ωJ k (ω) = πg 2 k δ(ω − ω k ) + δ(ω + ω k ) , and we generally have J(−ω) = −J(ω).This fixes the retarded/advanced components of the regulator.
Since it is furthermore desirable to keep a system in thermal equilibrium during the flow, if it was in equilibrium initially, one may also require the symmetry which is a sufficient condition for thermal equilibrium on the closed time path [59], where the transformation T β is defined as for our real degree of freedom φ(t).We can now insert the general ansatz for the regulator term and check that this condition is satisfied, if This implies that our fictitious Gaussian ensemble should then represent a heat bath at the same temperature T as that of the equilibrium system that is being regulated.
Our construction of causal regulators for thermal equilibrium systems therefore starts at specifying suitable FRG scale k dependent spectral densities J k (ω) to represent some fictitious heat bath.Here, we specifically use an analytic spectral density of the form (5.15) which has the desired regulating property by giving rise to an Ohmic bath with damping constant γ k = k/2 in the IR, while it rapidly goes to zero towards the UV, (5.16) The suppression of low-frequency modes is thus realized by the coupling to the fictitious heat bath with an FRGscale dependent damping constant which starts out large, of the order of Λ in the UV, and vanishes with k → 0 towards the IR.Inserting (5.15) into the spectral representation (5.11) explicitly yields for the corresponding causal heat-bath regulator in our example, R R/A HB,k (ω) = (5.17) with the Dawson function (5.18) Real and imaginary part of the retarded regulator are shown in Figure 2. Note, however, that its feature of representing a causal Green function unavoidably entails, In particular, whenever we construct a causal regulator in this way, the real part must have a zero crossing.It must start out negative for ω k in the UV, and with a single zero crossing as here, it thus turns positive for ω k in the IR.Note that a positive real part at ω = 0 corresponds to a negative mass-squared shift which can affect the propagator poles for sufficiently large k in an uncontrollable and unphysical way, a feature that is at least unpleasant for an FRG regulator.The easiest way out here seems to add a frequency independent counter-term that has the form of a Callan-Symanzik regulator.Such a frequency independent term is certainly causal, and it can offset R R/A HB,k (ω) by a k-dependent constant such that the real part of the resulting regulator stays strictly negative during the flow.We therefore introduce an additional positive parameter α to define (5.20) Since the absolute value of the real part of R R/A HB,k (ω) is monotonically decreasing for ω > 0, we restrict α from below by requiring α > R R/A HB,k (0)/k 2 .The effect of this can be understood in two equivalent ways: First, as every bath ensemble, cf.Eq. (4.12), the regulator heat-bath necessarily induces a negative shift of the system particle's mass or oscillator frequency squared, by −∆ω 2 HB (k), which here is k-dependent as our regulator bath is and can explicitly be written in the form With our analytic spectral density (5.15), for example, we have ∆ω 2 HB (k) = k 2 / √ 4π, and α must be chosen large enough to compensate this mass shift, αk 2 ≥ ∆ω 2 HB (k).This is necessary for the theory to remain causal during the flow in the first place; see the discussion that follows below.Secondly, in order to regularize all infrared modes we must have a negative real part of the regulator at ω = 0, requiring the strict inequality αk 2 > R R/A HB,k (0).
To further illustrate how the analytic structure of the propagators is changed by the heat-bath regulator, and how a suitable Callan-Symanzik counter-term can solves the issue, it is constructive to consider a spectral density for the regulator bath based on the Drude model, cf.Sec.IV B and Ref. [19], as a simpler alternative which yields with 2γ = ω D = k, (5.22) Using Eq. (5.11) again, we then obtain explicitly, where we already included the Callan-Symanzik counterterm for which we expect to require α > 1/2 = ∆ω 2 HB /k 2 .To see how the regulator affects the poles during the flow, we consider an exemplary retarded (advanced) propagator of the form The two are related by the symmetry The simplicity of the regulator (5.23) allows to derive analytic expressions for the poles ω p = ω p (k) by solving the cubic equation To keep the analytic structure intact, when the regulator is switched on, the retarded (advanced) propagator must only have poles in the lower (upper) half plane.For simplicity, we assume in our illustration here that ω 2 0 and γ stay constant during the flow.The resulting poles ω p (k) of the retarded propagator are shown in Figure 3 for different choices of α.For k = 0 the physical poles of the retarded propagator are located at in the complex plane, as represented by the black dots in Figure 3, corresponding to the expected quasi-particle excitations.We see that these two poles ω p (k) are located symmetrically around the real-part-zero axis.They move upwards with k towards the IR, coming from lower values of their imaginary parts, and thus always stay in the lower half plane and never cause problems with causality, even for vanishing Callan-Symanzik term with α = 0.For non-vanishing k > 0 there is a third pole with vanishing real part moving along the imaginary axis, however, which is entirely due to the regulator.In the time domain, it represents the purely relaxational contribution [18] here arising from our regulator heat bath.The crosses in the origins in Fig. 3 mark the point where it disappears for k → 0 + in the IR.The corresponding FRG Poles !p (k) of the retarded propagator (5.23) with the Drude regulator (5.22) in the complex plane as a function of k, with !0 = 1, = 0.5 and three di↵erent values of ↵.The black dots mark the quasi-particle poles (5.25) of the propagator at k = 0, where the regulator vanishes.They move around with k in the lower half-plane, but never cross the real axis.The crosses at the origin mark the points where the regulator-induced third poles disappear with k !0 in the IR.For ↵ < 1/2 this relaxational regulator pole moves into the upper half-plane at a finite value of the FRG scale k and the regulator thus violates causality at large k towards the UV.
Following the flow backwards towards the UV we see that for ↵  1/2 the imaginary part of the relaxational regulator pole first moves to smaller values.Eventually, however, it turns around to increase again towards the UV, for ↵ < 1/2 without bound.In this case it thus always crosses the real axis and moves into the upper half-plane (where a retarded self-energy should be analytic) so that causality is violated by the regulator at finite FRG scale k.For ↵ = 1/2 it turns around as well, but approaches 0 for k ! 1 in the UV and never moves into the upper half-plane.This is the liming case where ↵ is chosen precisely such that the regulator has a root at ! = 0.For larger values ↵ > 1/2 the imaginary part of the relaxational regulator pole decreases monotonically and the regulator never violates causality.Its real part has no zero-crossings anymore, and always leads to a pos-FIG.3. Flow of the imaginary parts Im !p (k) of the regulatorinduced relaxational poles in the retarded propagators of Fig. 2 over the FRG scale k.Here, ↵ = 1/2 is the limiting case, i.e. causality is always violated at large k for ↵ < 1/2.
itive mass/frequency shift, because the Callan-Symanzik counter-term is large enough to compensate the negative shift in the squared mass/frequency by ! 2 HB from the heat bath regulator.

C. Truncation for the E↵ective Average Action
In analogy to Ref. [28], we use the vertex expansion around the scale-dependent minimum 0,k (x) up to order Q, since it was proven that such a truncation gives rise to qualitative structures such as the collisional broadening and further resonance frequencies in the spectral function, corresponding to 1 $ 3 processes.
We will now briefly summarize the truncation and the di↵erences to the one presented in Ref. [28].For the quartic oscillator, the minimum 0,k (x) ⌘ 0 is kindependent because of the inversion symmetry of the e↵ective action and the assumption that no spontaneous symmetry breaking occurs.We consider a vertex expansion up to sixth order in the field , which may be explicitly written as FIG. 3. Trajectories of the poles ωp(k) of the retarded propagator (5.24) with the Drude regulator (5.23) in the complex plane, moving with the FRG scale k from the UV towards the IR, with ω0 = 1, γ = 0.5 and three different values of α.The black dots mark the quasi-particle poles (5.26) of the propagator at k = 0, where the regulator vanishes.They move with k in the lower half-plane as indicated by the arrows, but never cross the real axis.The crosses at the origin mark the points where the regulator-induced third poles disappear with k → 0 in the IR.Their flows are indicated by arrows as well.For α < 1/2 this relaxational regulator pole violates causality: It starts in the upper half-plane and only crosses at a finite value of the FRG scale k during the flow into the lower where it has a turning point (at the end of the line) before it moves up again and disappears in the origin.scale k dependence of its imaginary parts Im ω p (k) are shown for the same three values of α in Figure 4.
Following these from the IR towards the UV we see that for α ≤ 1/2 the imaginary part of the relaxational regulator pole first moves to smaller values.Eventually, however, it turns around to increase again towards the UV, for α < 1/2 without bound.In this case it thus always crosses the real axis and moves into the upper half-plane (where a retarded self-energy should be analytic) so that causality is violated by the regulator at finite FRG scale k.For α = 1/2 it turns around as well, but approaches 0 − for k → ∞ in the UV and never moves into the upper half-plane.This is the liming case where α is chosen precisely such that the total regulator in (5.23) has a root at ω = 0.For larger values α > 1/2 the imaginary part of the relaxational regulator pole remains strictly negative, and the regulator never violates causality (here, for α = 1 it decreases monotonically, in fact).Its real part has no zero-crossings anymore, and always leads to a positive mass/frequency shift, because the Callan-Symanzik counter-term is large enough to over-compensate the negative shift in the squared mass/frequency caused by ∆ω 2  HB from the heat-bath regulator.

C. Truncation of the Effective Average Action
Before we discuss the details of our truncation, we first give the explicit expression for the effective average action at the starting point k = Λ, where it equals the bare Keldysh action, Γ Λ = S.For the anharmonic oscillator (2.1) coupled to an external Ohmic heat bath with damping constant γ the Keldysh action is given by [17] S where we already used the Fourier transform of the quadratic part in the action for convenience.With this general structure of the bare action in mind, one possibility to truncate the Wetterich equation for the full effective average action is a functional Taylor expansion in terms of the 1-PI n-point vertex functions [28].Using the origin in field space, φ 0 = φ c 0 , φ q 0 = 0, as the expansion point, and the abbreviation such that the Q-point vertex is the highest one that is taken into account.
In this work, we use this vertex expansion (5.28) up to the order Q = 6 in combination with a loop expansion of the corresponding right hand sides of the flow equations. 3For the latter we adopt an ordering scheme tied to the vertex expansion in such a way that the highest n-point function (with n = Q) is assumed to be given by a frequency (and momentum) independent (but scale k dependent) vertex, while successively higher loops are included for the lower n-point functions.Specifically, we include (Q − n)/2-loop structures for the n-point functions with n = 2, . . .Q.With Q = 6 here, this amounts to taking into account the 2-loop structure of the 2-point functions, the 1-loop structure of the 4-point functions, and the scale-dependent constant 6-point vertex without substructure (corresponding to the order zero in the loop expansion).
Note that our combined vertex and loop-structure expansion at the order Q = 4 would essentially only yield a mass resp.frequency shift of the main peak in the spectral function, corresponding to the 0 ↔ 1, 1 ↔ 2 and higher one-step transitions.In order to describe effects such as collisional broadening or higher resonance excitation frequencies in the spectral function, as e.g.corresponding to 0 ↔ 3 or 1 ↔ 4 transitions, one needs non-local (here meaning frequency dependent) corrections of oneloop form in the 4-point function [36].In the combined scheme we adopt here, with Q = 6 this implies self energies of two-loop structure, and it then automatically also includes the local but k-dependent 6-point vertex which leads to a further quantitative improvement.
To explain the truncation in more detail, we start with the formal expression for the effective average action, where we denote spacetime integrations over x in short 3 In the symmetric phase of φ 4 -theories, without spontaneous symmetry breaking or tunneling in quantum mechanics for ω 2 0 > 0, all odd n-point functions vanish identically, and the minimum of the effective average action is fixed at φ 0 = (φ c 0 , φ q 0 ) = 0, independent of k.At this point it is convenient to follow Ref. [36] and to introduce the shorthand notations for convolutions of propagators and regulator insertions with fixed outgoing legs (c, q), (q, c) and (c, c), respectively.These are counterparts of the retarded, advanced and Keldysh propagators G R k , G A k and G K k with all possible ways of inserting one regulator term ∂ k R k in between.

6-Point Function and Effective Potential
Working out general flow equations for n-point couplings by the diagrammatic method can rather cumbersome, especially in the case of the 6-point coupling that we are interested in.In fact, it is much more convenient to first consider the flow equation for the scale dependent force from the effective potential V k (ϕ), here defined as a function of the rescaled classical field ϕ ≡ φ c / √ 2 by and generally valid with any ansatz for the effective average action Γ k [φ], where the prime denotes ordinary differentiation w.r.t. the constant classical field variable ϕ.In a general non-equilibrium situation, the potential might be spacetime x = (x 0 , x) dependent.This is not the case, however, for a spatially homogeneous system in thermal equilibrium.The definition in (5.32) is motivated by the form of the potential term S V [φ c , φ q ] in the bare Keldysh action on the closed time path, given by for a general potential V (ϕ) in the Lagrangian of the theory, where −V (ϕ) is the force in the classical field equations obtained from the φ q → 0 limit.For the spacetime independent vertices we are interested in the flow equations for the Taylor coefficients V (n) k (0) of the scale dependent effective potential V k (ϕ) used in (5.32), when expanded around a possibly likewise scale dependent minimum ϕ 0,k .In our case, ϕ 0,k ≡ 0 for all k, and we are interested, in particular, in the sixth order Taylor coefficient which is precisely our 6-point coupling constant, µ k ≡ V (6) k (0).Using Eq. (5.32) to define the scale dependent effective potential (up to a constant), we can then furthermore relate the desired Taylor coefficients to the spacetime integrals of the corresponding n-point functions, via (5.34)where we have also used the exchange symmetries of the n-point functions, Γ ...αβ... k (. . ., x, y, . . . ) = Γ ...βα... k (. . ., y, x, . . .), (5.35) which are valid specifically for a real scalar field theory [60], in order to combine equivalent terms in (5.34).
To obtain the flow equation for derivative of the effective potential, we thus have to project the Wetterich equation (5.1) on constant classical field configurations accordingly.To achieve this, we first take the functional derivative with respect to the quantum field φ q (x) on both sides of the Wetterich equation, and then set φ q = 0, φ c = const.which diagrammatically corresponds to the equation, see e.g.[37], (5.36) Due to the functional derivative, the flow of the zeropoint energy V k (0) is lost, of course.This in generally true on the closed-time path, however, where the Keldysh action contains no information on the zero point energy either, because the contributions to a constant offset in V (ϕ) from the forward and backward branches exactly cancel, cf.Eq. (5.33).For the flow of the higher Taylor coefficients of the effective potential, we can set φ q = 0 in (5.36), but we need to maintain the dependence on the constant classical field ϕ.This implies that one would need a partially field-dependent full 3-point vertex function in the loop diagram on the r.h.s. of (5.36) which obeys its own flow equation involving successively higher n-point functions as usual.At this point we employ a localvertex approximation in the sense that we neglect possible spacetime dependent substructures but maintain the required field dependence in the local part.Consistency with Eq. (5.34) then requires us to use, where the dependence on the constant classical field ϕ will be needed for the higher-order derivatives later on.Using this local-vertex approximation and the definition in (5.31c), the flow equation for the effective potential from (5.32) and (5.36) becomes, (5.38) The notation B K ϕ,k on the right indicates that analogously field dependent propagators G ϕ,k ≡ G k [ √ 2ϕ, 0] are being used inside the loop.To further illustrate the technique needed for further derivatives w.r.t. the classical field expectation value ϕ, first consider for example the fully field-dependent retarded propagator G R k [φ c , φ q ].For the successive Taylor coefficients of the effective potential we can again set the classical field φ c (x) = √ 2ϕ to its constant expectation value and the quantum field to zero, φ q (x) = 0.For the partially field dependent , 0] this implies that we can relate the ordinary partial derivative w.r.t. the constant field expectation value ϕ to its functional derivative w.r.t. the classical field φ c (x), To evaluate the functional derivative inside the integral, we make use of (5.7b) which, after applying the func-tional chain and product rules, directly tells us that (5.40) The dots in the arguments of the intermediate 3-point function indicate that the middle argument is fixed by the external differentiation point x, whereas the first and the third argument are convoluted with those of the propagators at the outgoing and incoming legs.After this functional derivative we can now insert our local-vertex approximation from Eq. (5.37) which then from (5.39) yields, This is the final derivative relation for the retarded propagator that we intended to illustrate here.Together with the FDR, we can now construct arbitrarily high ϕderivatives of the retarded, advanced and Keldysh propagators in the local-vertex approximation, ϕ,k (ω), and (5.42c) (5.42d) These relations form the basis for a set of recurrence relations, because we can now obtain flow equations for the higher derivatives of the scale dependent effective potential from (5.38) by iterating these relations (5.42a) -(5.42d).Setting the classical field variable ϕ to the expansion point (here at ϕ = 0) afterwards, then finally results in corresponding flow equations for the Taylor coefficients V (n) k (0) which are related to the n-point coupling constants via (5.34).The resulting flow equation for the sixth order Taylor coefficient

4-Point Function
In a real-time φ 4 theory there are three different types of 4-point vertices at 1-loop level, namely (a) the classical φ c φ c φ c φ q vertex, (b) the quantum φ c φ q φ q φ q vertex and (c) the 'anomalous' φ c φ q φ c φ q vertex [36].The former vertices (a), (b) already exist at tree level in the bare Keldysh action (5.27), and acquire (non-local) cor-rections during the FRG flow.In contrast, the anomalous vertex (c) does not exist at tree level and is first generated at 1-loop order.Since in our truncation scheme we want the flow of the 4-point function to be 1-loop exact, we have to consider all three vertices (a), (b) and (c).
We start with a few general remarks on how to truncate the flow equation for a 4-point function Γ The flow equation for each of these 4-point functions is obtained from a corresponding forth-order functional derivative of the Wetterich equation (5.1).This is straightforward but tedious, so it is not explicitly repeated here.The result, evaluated at the origin in field space (see e.g.[61]), can be compactly summarized [38] as follows, (5.43) At this level, when using scale-dependent local vertices on the right-hand side of the flow, there is a natural separation into s, t and u-channel contributions, corresponding to the first three diagrams in (5.43).We furthermore omitted the explicit labelling of the legs in the right-most diagram containing the 6-point function which is local anyway, at the order Q = 6 of our truncation scheme, and the order of the external legs hence irrelevant.With scale-dependent local vertices inside the flow, we can now apply our loop expansion by splitting the 4point function Γ αββ α k into a sum of s, t and u-channel contributions each of which are local in two of the three relative coordinates x − x , x − y and x − y , . Note that such an s, t and u-channel splitting is not generally possible beyond one-loop: Although (5.43) formally has a one-loop structure, with full 4 and 6-point vertices on the right-hand side of (5.43), it would contain structures of arbitrarily high loop order.In general, i.e. with nonlocal vertex functions inside the loops, Equation (5.44) would therefore represent an additional approximation.
Here we stick with scale-dependent local vertices inside the loops and decompose Eq. ( 5.43) into the three channels (5.44) which are all determined by essentially the same two-point correlation Γ αβ;β α k (x, x ) with suitably permuted indices and arguments.We are therefore free to arbitrarily select one out of the three equivalent channels.For the s-channel the flow equation reads, where an extra factor of 1/3 in front of the 6-point diagram occurs because we evenly distribute the local 6-point contributions to the s, t and u-channels, cf.Eq. (5.44).As full 4 and 6-point functions on the r.h.s. of the flow equation in (5.45) would invalidate this s, t and u-channel splitting, they are therefore beyond the truncation scheme used here.At the present order Q = 6 this scheme requires the flow of the 4-point function only to be one-loop exact, and we may therefore replace the full 4-point and 6-point vertex functions by effective local coupling constants ν k and µ k , respectively.For the 4-point functions on the r.h.s. of (5.45) this amounts to inserting, with effective local coupling constants ν αββ α k (see Eq. (5.48) below.)Moreover, at our truncation order Q = 6 the 6-point function is determined by the scaledependent local (classical) vertex V (6) k (0) = µ k already, cf.Eq. (5.34).Inserting (5.46) into (5.45), the r.h.s. of the flow equation (5.45) therefore becomes, where summation over repeated indices is implied again.The effective local coupling constants ν k are given by summing up all the (non-local) contributions of the 4-point function, (5.48) We can then insert the ansatz (5.44) into (5.48) and do a Fourier transform, where one can readily see that the ν's are just the sum over the the three channels of (5.44) evaluated at zero momentum, consistent with the 1-loop expansion of the flow equation.The relation (5.49) together with the flow equation (5.47) constitute our final system of one-loop exact flow equations for the full 4-point function Γ αββ α k in which the scale-dependent constant vertices are calculated self-consistently.
Having the general flow equations for the 4-point functions of type (a), (b) and (c) at hand, we need to briefly discuss one minor additional subtlety.Recall that the anomalous vertex (c) is first generated at one-loop order.It hence has a structure that is highly non-local in spacetime.The local approximation (5.48) is therefore not suited in this case, the anomalous vertex has no constant contribution at tree level.It is therefore in fact more accurate to set the effective coupling constant for this anomalous vertex (c) to zero, i.e. ν cqcq k = 0, on the r.h.s. of Eq. (5.47) at our truncation order, and only employ (5.48) for the classical and quantum vertices (a) and (b).We emphasize that we nevertheless still solve (5.47) for all of the three vertex functions (a), (b) and (c).Setting ν cqcq k = 0 to zero only effectively removes all diagrams from the r.h.s. of the flow equations (5.47) that would represent a constant contribution from the anomalous vertex.As a result, the anomalous vertex is not fed back at all into the flow equations for the 4-point functions.One would have to go beyond the one-loop expansion on the r.h.s. of the flow equations (5.47) in order to do this in a consistent way, see [36].In our combined vertex and loop-expansion scheme that would imply to increase the truncation order to at least Q = 8.However, the calculated anomalous 4-point vertex function will enter the flow of the two-point function, and can thus not be neglected entirely, although it does not re-enter the flow equations for the 4-point functions themselves.
We conclude this Subsection by introducing a few further notations that will be convenient in the following Subsections.For the two-point vertex functions we define the shorthand notations , and V an k (x, x ) ≡ Γ cq;cq k (x, x ) , to emphasize the distinction between the classical (a), quantum (b) and anomalous (c) vertex functions.We furthermore introduce the diagrammatic notation to denote the two-point correlation Γ αβ;β α k (x, x ) in each channel of the full 4-point vertex function, to emphasize its one-loop structure.We will use this notation in the next Subsection to represent the two-loop order in the flow equation of the two-point function Γ αα k (x, x ).

2-Point-Function
The exact flow equation of the 2-point function at the minimum φ 0 = 0 has the formal structure of a tadpole diagram, except that it contains the full 4-point vertex function [36,38], According to our truncation scheme at order Q = 6, we need to solve the flow equation for the two-point function in a two-loop exact way, which we achieve by inserting our one-loop exact 4-point functions (5.44) from the previous section without further approximation into the r.h.s. of the flow equation (5.52).With the notation just introduced in (5.51) above, this yields The s and u-channel contributions (the second and the third diagram in (5.53)) are responsible for generating a dynamic frequency dependence in the two-point function, and therefore in the spectral function, likewise.This is why we had to go to at least the order Q = 6 in the first place, to see effects in the spectral function that are beyond a constant mass/frequency shift.These contributions generate a dynamic frequency dependence since they explicitly depend on the external frequency (and momentum in higher dimensions) that flows through the diagram.In contrast, the t-channel contribution (first diagram in (5.53)) is proportional to ∼ δ(x − x ) and thus only contributes a constant shift to the bare mass/frequency ω 0 .Translating the diagrams into formal expressions then finally leads to here for the advanced A and Keldysh K components of the 2-point function.Together with the flow equations for the 4-point vertex functions (see (5.49) and (5.47)), and the scale-dependent 6-point vertex (see Appendix B 3) they now represent a closed system of flow equations.We solve this coupled system of flow equations numerically with a method outlined in Appendix B 1, and then calculate the spectral function in the IR (k → 0) via (5.8).As we have explained in the beginning, this system of flow equations is then fully consistent with our expansion order Q = 6 in the sense that the flows of the two-point functions are two-loop exact, that of the 4-point functions are one-loop exact with self-consistent tree-level contributions, and that of the 6-point function, being a scale-dependent constant, is self-consistent at tree level.
At this point we note for completeness that the classical-statistical limit is readily implemented in the real-time FRG flows by simply deleting the quantum φ c φ q φ q φ q vertex from the microscopic action (5.27) together with replacing the quantum distribution function by the classical Rayleigh-Jeans distribution, which amounts to replacing coth(βω/2) → 2T /ω for ω T , e.g.see [17].

UV Initial Conditions
In the UV (at k = Λ) the effective average action is given by the bare action (5.27), from which we can therefore read off the initial conditions for two, four and 6-point functions, ) and µ Λ = 0 , with the bare damping γ Λ , the inverse temperature β = 1/T (entering through the FDR), the bare coupling constant λ Λ , and the bare frequency ω 0Λ .

VI. RESULTS
For the discussion of the results, we set the bare mass/frequency ω 0 = 1 to unity, which corresponds to implicitly expressing energies in units of ω 0 and all dimensionful parameters in appropriate powers thereof (e.g. the coupling λ is here measured in units of ω 3 0 ).We furthermore use a rather small damping of γ = 0.06, in order to warrant the reliability of our benchmark solution from the exact diagonalization.
We start to discuss our results at a rather high temperature of T = 32 (in units of ω 0 ) in Figure 5.The high-temperature spectral functions in principle contain all possible energy differences E n − E m in the spectrum which are allowed by selection rules, e.g. as due to the conserved parity in our case, weighted by the appropriate factor e −βEn | n|x|m | 2 that quantifies the probability for the transition in the thermal mixed state of the canonical ensemble, cf.Eq. (2.9).At a small value of λ = 1/32 for the coupling in the quartic anharmonicity in panel (a) on the left, the system behaves nearly harmonically, and all approaches more or less coincide.The classical approximation for T ω is well justified.The main peak has a Breit-Wigner-like shape which arises from the thermal ensemble of one-step transitions |n → |n + 1 , see FIG. 5. Comparison of high-temperature spectral functions, at T = 32 with weak damping γ = 0.06, over frequency (all in units of ω0) from the various methods.The panel in (a) on the left shows results at a weak coupling of λ = 1/32, and that in (b) on the right the corresponding ones at a rather strong coupling of λ = 4.The sharp individual peaks from the quantized transition energies, cf.Fig. 1, gradually build up the broad continuum distributions observed at high temperatures in (a).These represent the classical limit in which the classical-statistical spectral functions agree with the GSA results (static/adiabatic), and all of them coincide with the solution from the exact diagonalization.Increasing the coupling at fixed temperature increases the splitting between the transition energies so that the individual peaks reemerge in (b).The main peak represents the ensemble of one-step transitions |n ↔ |n + 1 , the second one that of the three-step transitions |n ↔ |n + 3 at higher excitation energies.The time-dependent second moments σxx(t), σxp(t), and σpp(t) beyond the static approximation in the GSA produce contributions which are dismissed when extracting the spectral function from the quasi-classical method described in Sec.IV B 6.These contributions, which can be neglected in the nearly harmonic system (a), are mainly responsible for the differences between static and adiabatic GSA results in (b).The classical limit also serves to assess the truncation used in the real-time FRG calculations (performed on a frequency grid with 512 points in the interval ω ∈ [0, 15]).Whether or not explicitly employing the classical limit in the real-time FRG flow equations makes no noticeable differences here.The corresponding NLO 2-loop perturbative results are shown as dashed lines for comparison, and agree with those of Ref. [36].
Fig. 1, where the |n are distributed according to the Boltzmann weight e −βEn , cf.Eq. (2.9).Because the coupling is small, the transition energies for the different n that contribute are all close to that of the groundstate transition which itself is only a little larger than ω 0 = 1 in the harmonic case.The central frequency ω c of this rather narrow main peak at ω c 1 is therefore also close to ω 0 .All non-perturbative methods describe this main peak in perfect agreement with the exact solution.In contrast, the results from NLO 2-loop perturbation theory (taken from Ref. [36]) show spurious double peak structures.Most importantly, such a splitting of the main peak does not occur in the classical-statistical limit.Although the quartic coupling λ = 1/32 is rather small in Fig. 5 (a), this is not surprising, however, because the relevant thermal coupling λT = 1 is comparatively large: In fact, one can rescale variables in the classicalstatistical theory to trade the explicit dependence on the quartic coupling λ for a dependence only on the combination λT (which is not possible in the full quantum theory).Hence, the effective expansion parameter of the classical-statistical theory in the perturbative series is not λ, but the thermal coupling λT .Since λT = 1 in Figure 5 (a) and even λT = 128 in Figure 5 (b), we must not expect the perturbative expansion to be valid here.This once again emphasizes the need of non-perturbative real-time methods here, and it might help to appreciate the huge qualitative improvements brought about by the FRG.This is particularly reassuring for field-theory applications beyond the classical-statistical limit where we neither have exact solutions nor ab-initio results from real-time simulations.
The second peak represents the corresponding thermal ensemble of three-step transitions, |n → |n + 3 .It would be absent in the harmonic case and is therefore small because λ is.In the NLO 2-loop perturbative calculation it occurs at ω = 3ω 0 corresponding to the unperturbed energy difference in all these transitions.Diagrammatically it originates from non-local 'sunset' diagrams at 2-loop level which contain three bare propagators.Cutkosky's cutting rule then implies that the spectral function, i.e. the imaginary part of the retarded propagator, peaks at ω ≈ 3ω 0 .It is therefore neither at the correct frequency nor of the correct shape.The perturbative calculation overestimates its height and underestimates its width.This is a manifestation of the fact that the perturbative expansion is not valid here, because the effective expansion parameter of the classicalstatistical theory, λT = 1, is not small.
Beyond perturbation theory, the spread of the individual transition energies generally grows with ω and the second peak therefore tends to be wider than the main peak as well.The resonance frequency of the second peak is somewhat larger but close to 3ω c (which is still well below T = 32 here).In fact, in our present truncation to the real-time FRG flows, the second peak has to be at ω ≈ 3ω c , as discussed below.The fact that its central frequency is thus somewhat underpredicted by the realtime FRG in Fig. 5 (a) is therefore a first hint at the systematic errors due to the truncation.All other methods (except perturbation theory, of course) are in perfect agreement with the exact diagonalization result around this second peak as well.In particular, we observe no noticeable differences between the static and adiabatic GSA results.This also implies that extracting the Gaussian spectral functions quasi-classically, as described in Section IV B 6, is justified because the spectral function of the nearly harmonic system depends almost solely on the evolution of the expectation values X and P .
Results at the same temperature T = 32 but a considerably larger coupling of λ = 4 are shown in Figure 5 (b).This fairly strong coupling increases the splitting between the individual transition lines, whose widths are due to the heat-bath coupling with the same small γ = 0.06 as before, in each of the two ensembles so that the two corresponding peaks are broadened significantly and their substructure becomes clearly visible in the exact-diagonalization solution.For the same reason, these peaks have moved up in energy to about 3ω 0 and 9ω 0 , respectively, which also implies that the classical limit T ω is less well satisfied here, especially in second peak.As compared to Figure 5 (a) this second peak is more pronounced, on the other hand, because the matrix elements for the three-step transitions are larger at larger coupling.
The real-time FRG, at this order Q = 6 in the truncation outlined in Section V, has notable problems reproducing the classical limit for large coupling as also seen in Figure 5 (b).Its main peak is so broad that it almost swamps the second peak.To improve the truncation, on one hand, one has to go to higher orders Q ≥ 8 in our combined vertex and loop expansion.On the other hand, self-consistent solutions might also be important, as shown for example in Ref. [36], where the 4-point function was flowed self-consistently.Quite expectedly, this turned out to increase the analytical as well as the numerical effort tremendously.Unfortunately, it also led to a decreased numerical stability of the flow at the same time, however.As shown in Ref. [36] the relevant parameter in the loop expansion is the thermal quartic coupling λT .It was observed there already that the FRG spectral functions start to deviate from the classical-statistical result for rather small but certainly non-perturbative values of λT ≈ 4 when the 4-point function of 1-loop structure is employed.This should be compared to the comparatively huge value λT = 128 in Figure 5 (b), which exceeds this proposed range of validity by nearly two orders of magnitude.As visible in Figure 5 (a) and also below in Figure 6 (d) for comparatively small but already non-perturbative values of λT ∼ 1 − 2 the real-time FRG is still capable of accurately describing the shape of the spectral function.Hence there seems to be a limiting value λT ∼ 4 for our combined loop and vertex expansion at the present order as well.At this value e.g. the splitting of the 4-point function in s, t and u channels is no-longer sufficient, and higher loop structures have to be taken into account.Compared to the perturbative results this is a tremendous improvement, however, as already discussed in relation with Figure 5.The NLO 2-loop result is not able to describe the spectral functions reasonably well for any of the parameters that we have considered.Despite its difficulty with the necessarily very large thermal couplings λT of the classical limit, the real-time FRG is thus nevertheless very well suited to describe non-perturbative phenomena at lower temperatures where the classical-statistical calculations have to break down, eventually.
Due to the sufficiently high temperature, as in Figure 5 (a), the classical and the static GSA spectral functions still coincide, with minute differences due to the higher relevant frequencies here only in the second peak.They both interpolate the substructure of the exact solution and agree with the average shape very well.However, the GSA in the adiabatic approximation, with the non-trivial time evolution of the Gaussian widths included, now shows visible deviations from those averages in the second peak of the three-step ensemble.This is because for a coupling as strong as λ = 4 used here, our quasi-classical approach of extracting the spectral function described in Section IV B 6 becomes inconsistent, which we can infer from the fact that the classical spectral function better matches the exact one than the adiabatic GSA result does.Beyond the static approximation, the time-dependence of the Gaussian widths σ xx (t), σ xp (t), and σ pp (t) produces contributions to the spectral functions that are not contained in the quasi-classical extraction scheme based on the classical unequal-time correlators of X(t) and P (t) alone.This explicitly demonstrates the relevance of the discussion at the end of Section IV B 6.
Our results for the same strong coupling of λ = 4 (and γ = 0.06) but now at successively lower temperatures between T = 4 and T = 0.5 are summarized in Figure 6.We leave out the GSA results from the adiabatic approximation with the time-dependent widths in this summary because they do not represent improvements over the static GSA results.A dedicated comparison of the two approximation schemes for the GSA is deferred to Figure 7 and discussed below.
The two ensembles of one-step and three-step transition lines in the exact spectral functions of Figure 6 become sparser with decreasing temperature because the transitions with larger n get sequentially more suppressed.Because that removes strength form the higherfrequency side of each of the two, the corresponding peaks become narrower and their central frequencies shift to smaller values with decreasing temperature.In contradistinction to the high-temperature limit and the previous Figure 5, for very low temperatures these ensembles can not be represented by broad peaks in the first place anymore, the quantization in terms of the individual transition lines can eventually no-longer be neglected.With decreasing temperature the contributions from higher states in the one-step and three-step transitions get exponentially suppressed more and more so that the corresponding ensembles of main and second peak in the spectral functions get compressed towards their lowest transition frequencies until only a few individual transition lines remain.The (static) GSA spectral function follows the ensemble averages of main and second peak more closely than the classical one as temperature is lowered beyond the range of validity of the classical approximation (T ω).In our present truncation of the real-time FRG, the second peak for ω ≈ 3ωc stays put at about three times the central frequency ωc of the main peak (a frequency grid with 320 points on ω ∈ [0, 8] was used for T = 1, 2, 4, and one with 200 points on ω ∈ [0, 5] for T = 0.5).As in Fig. 5, the corresponding NLO 2-loop perturbative results are included with dashed lines for comparison.
While all methods reproduce the general trend of the overall infrared shift and reduction of width of the main peak in the spectral function, there are quantitative differences worth discussing: At the starting temperature of T = 4 in Figure 6 (a), classical and GSA spectral functions are still very similar and both agree well (on average) with the exact result, while the FRG solution shows problems analogous to those discussed in relation with Figure 5 (b) before.Reducing the temperature to T = 2, we see in Figure 6 (b) that the classical and GSA spectral functions start to separate from each other.The classical result tends to underestimate the central frequency of the main peak, and the GSA more closely follows its shape on average.In particular, the GSA tends to reproduce better its rather abrupt start due to the relatively sharp and well isolated quantum mechanical ground-state transition line on the low-frequency side.
When it comes to the second peak representing the three-step transition lines, the GSA result is clearly able to describe their ensemble average significantly better than the classical spectral function.This trend continues towards lower temperatures where the GSA results are able to follow these ensemble averages more closely than the classical spectral functions, and show an enhanced strength in the second peak as compared to the classical-statistical result.Consequently, all these effects become more pronounced at T = 1 in Figure 6 (c).At the lowest temperature T = 0.5 of this comparison in Figure 6 (d) it becomes obvious that the GSA in our static approximation eventually also tends to underestimate the central frequency of the second peak in the Caldeira-Leggett model.
The classical spectral function is bound to approach its mean-field value in the limit T → 0, which is here simply given by the Breit-Wigner form with width γ around the unperturbed ω 0 .It should therefore only be considered valid for high temperatures.In contrast, the exact vacuum spectral function (for T → 0) in the interacting theory still contains all possible ground state transitions |0 ↔ |1 , |3 , |5 , . . .which are inherently quantum mechanical.Because even the lowest energy difference, between ground and first excited state, increases in presence of the quantum-self-interactions which are not included in the classical limit, the latter thus necessarily fails to describe the low-temperature mass shift correctly.We also notice that it produces a main peak which is systematically too broad in comparison with the exact solution.
In the language of the closed time path and the Martin-Siggia-Rose (MSR) path integral formulation of classicalstatistical mechanics, it is missing a quantum φ c φ q φ q φ q vertex.The real-time FRG, which includes such a vertex, therefore becomes better for smaller temperatures, and the location of its main peak, with central frequency ω c , fits the quantum-mechanical solution quantitatively quite well.The strength from the higher excitations in the second peak is best reproduced by the real-time FRG as well, although its central frequency stays closer to the one in the classical spectral function, as most prominently seen in the T = 0.5 plot in Figure 6 (d).The reason for this can be understood from the truncation described in Section V.The effective mass-shift of the main peak is indeed generated by the tadpole diagram in (5.53).Responsible for the second peak, however, are non-local 'sunset' diagrams, e.g.see Ref. [36].They contain three Green functions which all include the correct effective mass-shift as explained above.The imaginary part of this diagram from Cutkosky's cutting rule then fixes the location of the second peak in the spectral function to ω ≈ 3ω c .This could best be improved upon with self-consistency and by including higher-order vertex corrections.
To further assess the influence of the static approximation in the GSA, we also ran the adiabatic GSA simulations at our lowest temperature T = 0.5.The results are shown in Figure 7.We observe that the adiabatic approximation of the GSA quantitatively yields a slightly better estimate of the central frequency in the main peak.An additional structure on its right side, as shown in the inserts, might even be interpreted as a remnant of the individual transition-line substructure.Such a process is not contained in the classical-statistical approach.
The three-step transitions at higher frequencies are captured by both of the GSA results.In the adiabatic version of the GSA a double-peak structure thereby emerges resembling the corresponding substructure in the exact diagonalization result which is a purely quantum-mechanical effect (and is never observed in the classical spectral function, cf. Figure 6 (d)).Quantitatively, however, both GSA results underestimate the exact frequencies, and appear to too broad as compared to the exact diagonalization solution.The continuous high-frequency fall-off of the spectral function from there on is in fact described best by the static GSA result which includes the quasi-classical correction factor from the colored quantum noise, cf.Section IV B 5. In contrast, the large-ω behavior of the adiabatic GSA result does not and low temperature (T = 0.5).With the adiabatic corrections included in the time evolution, the central frequency of the main peak matches the exact one slightly better than in the static approximation.An additional 'flank' on its right (shown in the insert) resembles the multi-peak substructure of the individual transitions from the exact diagonalization.Moreover, we observe a splitting of the second peak in the adiabatic approximation, which at least qualitatively resembles the first two distinct three-step transitions in the exact diagonalization.
match that of the exact solution so precisely anymore.This is another indication that the quasi-classical method to measure the spectral functions as described in Section IV B 6 eventually becomes inconsistent with the full time evolution of the Gaussian widths in the adiabatic GSA at very low temperatures (or very high frequencies), in agreement with the discussion of the classical limit of the GSA in Figure 5 above.

VII. CONCLUSION & OUTLOOK
In Sections II to V we have established four different real-time methods for calculating spectral functions.These are based on exact diagonalization, classicalstatistical field theory, the Gaussian-state approximation (GSA), and the functional renormalization group (FRG) formulated on the Keldysh closed-time path.We have compared the results from these various methods in Section VI for the spectral function of the quartic anharmonic oscillator coupled to an Ohmic heat bath at finite temperature as in the Caldeira-Leggett model.
Having established the underlying quantum mechanical structure of the spectral function for the anharmonic Caldeira-Leggett oscillator from exact diagonalization in Section II, we were able to demonstrate explicitly how the quantum mechanical system with its discrete transition lines gradually turns into a classical one when temperature is increased, as illustrated here by the changing behavior of the spectral function.The quartic anharmonicity with coupling strength λ of the system is responsible for the emergence of distinct thermal ensembles, each of which consists of individual transition lines whose separation increases with λ.With increasing temperature T , more and more transitions between states higher up in the spectrum contribute to each ensemble.It then is the ensemble averages that tend to the continuous broad peaks of the spectral function in the classical limit T ω.Although computationally expensive, such an approach could be used to calculate spectral functions for field theories or at least small many-body systems exactly, which is an interesting topic for future work.Here it provides our benchmark solution.
We have summarized the most important concepts needed for classical-statistical simulations in Section III.These can be used for ab-initio calculations of spectral functions in the classical limit of large bosonic occupation numbers.This limit is realized at high temperatures or low frequencies, i.e. for T ω, and it can therefore also be used to study the static and dynamic universal behavior near a thermal phase transition at T c of the critical infrared modes with ω T c [21][22][23][24].Here it served us as the limiting case, to verify the validity of the other approaches in the high-temperature limit.We have observed that in this limit the classical result describes the exact spectral function perfectly well as a coarse-grained version of the quantum mechanical one.However, the classical approach is not able to resolve the quantummechanical substructure from the individual transition lines, and it also underestimates the frequencies of main and second peak in the spectral function when it reaches its limitations at smaller temperatures.
In Section IV we have shown how to construct a consistent description of an external heat bath in the Gaussian state approximation, describing the effects of thermalization, dissipation, and quantum-mechanical decoherence.There, we have employed two distinct approximations, named 'static' and 'adiabatic', and discussed their effects on the spectral functions in Section VI.We have verified that the GSA is perfectly consistent with the high-temperature limit, where it reproduces the classical results as expected.In the static approximation the GSA is able to follow the ensemble averages of the exact solution towards lower temperatures than the classical spectral functions, thus extending their range of applicability by some margin.While the GSA can serve as a useful qualitative tool to study the existence and rough structure of quantum effects as (small) corrections to classical-statistical simulations towards lower temperatures, it leaves the discrete substructure in the ensemble peaks unresolved, likewise.An interesting opportunity for future studies would be to investigate alternatives to Gaussian distributions of pure states, e.g. based on the logistic function as suggested in [62], for the system particle at large anharmonicities where already the leading corrections to the classical limit are expected to be non-Gaussian.
While the GSA in the static approximation presented here might not be so well-suited for high-precision calcu-lations at strong coupling λ, it does incorporate the exact classical-statistical field theory limit near a thermal second order phase transition at a critical temperature T c , where the dynamics are dominated by critical infrared modes.However, at a finite 'distance' to the critical point or in non-equilibrium phase transitions along trajectories in the phase diagram that get close to it, quantum mechanical effects might well become important.In such a situation, as relevant e.g. in heavy-ion collision experiments searching for the QCD critical point, the GSA could serve a useful indicator for that to happen.
For our real-time FRG calculations of spectral functions in Section V we have introduced the novel concept of heat-bath regulators.These are constructed from coupling the system to a fictitious external heat bath, which is introduced in the spirit of the Caldeira-Leggett model as an ensemble of harmonic oscillators, controlled by an FRG scale k dependent spectral density J k (ω).This provides a rather intuitive picture of suppressing infrared modes by overdamping, and it is particularly wellsuited for near-equilibrium real-time calculations.The construction includes regulating real and imaginary parts of self-energies, while the causal structure of the Keldysh action is built-in.At the same time, this causal structure of such a regulator added to the Keldysh action necessarily also requires adding a Callan-Symanzik counterterm as well, in order to avoid acausal regulator singularities and to suppress long-wavelength infrared modes.With these counter terms included, however, our heatbath regulators can straightforwardly be generalized to field theories in d spatial dimensions, which is work ongoing.An analogous construction scheme might also be feasible, for example by imagining the bath to be an ensemble of Klein-Gordon fields, which brings along further subtleties such as e.g.Lorentz covariance, that are not present in our 0+1 dimensional example here.
After establishing our causal regulators, we have adopted a truncation scheme for real-time FRG calculations, which is a modification of that used in Ref. [36] in that it combines vertex and loop expansions.At our present truncation order it includes a self-consistent vertex expansion of the self-energies up to two-loop order.We have demonstrated in Section VI that such a truncation yields robust and quantitatively very reasonable results for the main peak in the spectral function at low temperatures.Noticeable problems occurred only in reproducing the classical limit at strong coupling λ, which indicates that higher-order non-local structures in the vertex functions might have to be taken account to describe this regime in parameter space.
With regulators and truncation scheme for real-time FRG calculations in place, the potential of this new direction of field-theory applications for the future includes studying dynamic critical phenomena and nonequilibrium phase transitions; work we have not touched upon here but which is also under way already.
First FRG calculations of dynamical critical exponents z for various dynamical models [39][40][41] or criti-cal spectral functions [37] have only relatively recently become available and will be systematically expanded.Moreover, the approach of building causal regulators might be analogously possible also for fermions, leading to the perspective of studying different dynamics in renormalizable chiral effective theories such as the Quark-Meson Model, as an alternative and extension to the analytically continued Euclidean FRG flows [8][9][10]12], or the chiral Parity-Doublet Model for nuclear and chirally symmetric hadronic matter [11,63] within our real-time FRG framework.
For non-equilibrium phase transitions, both the realtime FRG and classical-statistical simulations are wellsuited to study off-equilibrium phenomena such as finitetime trajectories in the vicinity of a critical point in the phase diagram [64].It is believed that the theoretical understanding of such systems might yield protocols to locate the QCD critical point through the data obtained from heavy-ion collisions (see e.g.[65][66][67][68]).Both the classical-statistical simulations, possibly extended by the static GSA, as well as the real-time FRG provide promising computational frameworks for further investigations of such non-equilibrium systems in the future.and its inital value by C 0 = C(0).
Putting everything together, and combining (A.5a), (A.5b) with (A.5c), (A.5d) then yields two decoupled second-order differential equations for the irreducible correlations of x(t) with the initial heat-bath coordinates and momenta.These differential equations describe driven harmonic motion with damping (including memory effects) and, most importantly, with an in general time-dependent frequency given by the square root of the curvature C(t), For low frequencies ω s Λ, on the time-scales relevant for the driving force, the memory integrals over the damping kernel reduce to ordinary (local in time) damping terms γ d dt G xϕs (t), γ d dt G xπs (t), cf.Eq. (4.19).An exact analytic solution to the general oscillator problem with time-dependent restoring force C(t) is unfortunately not known, at least to us.Therefore, we have to resort to an additional adiabatic approximation, assuming that the curvature of the potential fluctuates slowly about a temperature-dependent equilibrium value C 0 (T ) obtained from a mean-field prescription, To go beyond this approximation, we furthermore split the time-dependent value of the curvature into this constant equilibrium value plus a small perturbation, and treat δC(t) as a correction to the exactly solvable Eqs.(A.11a), (A.11b) for the constant C 0 ≡ C 0 (T ) (dropping the temperature dependence in the following) in time-dependent perturbation theory with inital condition δC(0) = 0 so that the equilibrium value C 0 is our initial value for C(t) at the same time.
This removes the UV divergence in the stochastic forces from both equations, (4.36b) and (4.36c) at the same time.The constant stochastic force K pp on σ pp in Eq. (4.34b) is then absorbed completely in the static approximation, and the subtracted force in (4.36b) reads This is a rather simple set of inhomogeneous first-order differential equations with the particular solution, for weak damping γ < 2 √ C 0 , when all homogeneous contributions have died out for t → ∞, σ xx (t) → F (C 0 ), σ xp (t) → 0, (A.28) σ pp (t) = σ pp (0) + σ r pp (t) → σ pp (0).

b. Inital Conditions
The vanishing constant force on the subtracted σ r pp in the last line of this static approximation is consistent with the initial condition σ pϕs (0) = 0 , cf. (A.3) for t = 0.Because of the constant subtracted force (A.24) on σ xp , however, we have implicitly assumed non-trivial initial conditions for the σ xϕs .Suitable initial conditions that are consistent with the constant (subtracted) K xp can be read off from Eq. (A.2) for t = 0, With σ pϕs (0) = 0 and σ pπs (0) = 0 the time derivative of Eq. (A.2) for t = 0 then implies that σ xπs (0) = 0 also.
In summary, the presence of the constant force on σ xp in (A.26b) requires us to slightly modify the initial covariance matrix in (4.23).As we switch on the coupling between the heat bath and the particle at t = 0 we have implicitly assumed here that this happens via switching on the off-diagonal σ xϕs via the term σ xϕs (t) = Θ(t)σ xϕs (0 + ) + Assuming that the curvature of the potential C(t) varies slowly in time compared to the relevant bath oscillator frequencies we can elevate the static approximation to an adiabatic one: Using the driven oscillator solutions in Eqs.(A.14a) and (A.14b) for the fast bath oscillator degrees of freedom, we may then simply replace C 0 by C(t) in the final expressions for the equations of motion.The stochastic forces then become time dependent as well, but they remain ultraviolet finite also for the Ohmic bath J Λ (ω) with Λ → ∞, after the time-independent subtraction of the static approximation.The adiabatic approximation then yields for the Gaussian widths, with ω j = −L + jδ for j = 0, . . ., 2(N − 1).The numerical implementation also exploits the symmetry relation V qu k (ω) = V cl k (−ω) to only calculate and store the classical and anomalous vertices.Using these approximations the coupled set of flow equations consisting of partial and ordinary integro-differential equations is reduced to a finite set of ordinary integro-differential equations for the variables Γ (2),A k,j , V cl k,j , V an k,j for j = 0, . . ., 2(N − 1), and µ k .(B.2) For simplicity, the resulting differential equations are numerically solved using an explicit Euler method.To increase the efficiency, the equations are reformulated using the RG 'time' parameter t = log(k/Λ).In general, convolution integrals of the type with Ω fixed, have to be solved numerically, which are performed using a trapezoidal rule for integration.Doing this for N grid points has a numerical complexity of O(N 2 ).This may be optimized by performing a Fast Fourier transform (FFT) to real space, where the convolution integrals are simply multiplications, and then translating the result back into momentum space using an inverse FFT, as it is done in Ref. [36].Further, the convolution requires g to be evaluated outside the region where the discrete data is available.Therefore, the function value g(ω) must be extrapolated from the discrete values g j = g(ω j ).Here it is important to roughly know the behaviour of g(ω) for |ω| → ∞.
If g(ω) → ±∞ for |ω| → ∞ (i.e. for g = Γ (2),A k ), a second order Taylor expansion at the boundaries is employed, where the derivatives are evaluated using finite differences with a single-sided 3-point form for the first derivative and a single-sided 4-point form for the second derivative.
If g(ω) → const for |ω| → ∞ (i.e. for g = V cl k , V an k ), the function value g(ω) is just set to the limit lim ω →±∞ g(ω ).The drawback is that this constant must be known a priori, which is in generally not the case, but can be motivated from the UV limit (5.56d), replacing λ Λ by the corresponding effective coupling following Eq.(5.48).This method works best if L is chosen large enough that g(ω) lim ω →±∞ g(ω ) for |ω| > L is a sensible approximation.For all numerical calculations the flow parameters are characterized by Λ = 20, k IR = 0.01.In every case the ω-grid is chosen large enough to cover the full range of all relevant excitations and fine enough to properly resolve individual peaks of width γ.

Regulator Dependence
An important test for a given truncation is how sensitive it is to the specific choice of the regulator.In the optimal case, the effective average action lim k→0 Γ k in the IR should not depend on the regulator.But for any finite truncation a spurious dependence on the regulator will certainly be introduced.To analyze this dependence, we test two possibilities of causal regulators, namely (i) a linear combination R HB,k −αk 2 of a heat bath regulator defined as (5.11) and the Callan-Symanzik regulator (in the following referred to simply as 'heat bath regulator' or in short 'HB+CS') and (ii) a Callan-Symanzik-like regulator where the damping is uniformly increased by 2k for all frequencies, i.e.R

R/A k
(ω) = −k 2 ± 2ikω (in the following referred to as 'CS+Damping regulator'.) For convenience, the heat bath regulator is reparameterized according to α = α 0 + α , where the α 0 -term eliminates the negative mass shift quadratic in k from the heat bath regulator and the α -term then adjusts the regulating mass.To be explicit, we have α 0 = 1/ √ 4π for the regulator in equation (5.17).
For the parameters λ = 4 and γ = 0.06 and two different temperatures T = 0.5, 4 a comparison between the resulting spectral functions is shown in Fig. 8.We also included the comparison with the FRG result when applying the classical limit, i.e. deleting the quantum φ c φ q φ q φ q vertex and replacing the quantum distribution function by the Rayleigh-Jeans distribution, coth(ωβ/2) → 2T /ω [17].It becomes clear that the results indeed depend on the specific choice for the regulator.
We first focus on the T = 4 case, visualized in Fig. 8  (a), where we know from the discussion in Section VI that the main peak is quantitatively well described by the classical result.Taking this classical-statistical spectral function as a benchmark for the FRG there, we see that the qualitative structure of the main peak is better described by a lower value for α .In contrast, the second bump in the spectral function at around ω ≈ 6 is better described for a higher value of α .Therefore α should We also included the result for the classical limit in the realtime FRG, which arises by deleting the quantum φ c φ q φ q φ q vertex and replacing the quantum distribution function with its Rayleigh-Jeans counterpart [17].It agrees for T = 4 very well with the full quantum result from the FRG, indicating that there we already are close to the classical limit.For T = 0.5 it differs from the quantum result and matches the classical-statistical spectral function better, as expected.(For the FRG we use an N = 320, L = 8 grid for T = 4, and an N = 200, L = 5 grid for T = 0.5.) be chosen such that both the main peak and the second peak are both described sufficiently well, resulting in an optimization problem for α . 4 For the CS+Damping regulator the main peak is slightly more sharpened, reducing the agreement with the classical-statistical spectral function in comparison with the heat bath regulator.The CS+Damping regulator also shows a second bump, but which is much less pronounced than in the case of a heat bath regulator.This 0+1 dimensional example already indicates the importance of using a heat bath regulator in practical calculations.For completeness, we note that the classical limit of the FRG agrees very well with the full quantum result, which indicates that T = 4 is indeed a temperature where the system behaves classically.
We now turn to the low temperature (T = 0.5) case, shown in Fig. 8 (b), where we observe qualitatively the same regulator-dependence as in the T = 4 case in Fig. 8 (a).The resulting spectral function obtained from the FRG shows a very weak dependence on the Callan-Symanzik counter-term parameter α .It also shows a slight broadening of the main peak and a minor enhancement of the second bump for increasing α .The CS+Damping regulator produces a sharper main peak like in the T = 4 case, and further a decrease in the height of the second bump.This suggests that the effects of the regulator noted above are stable under the variation of the temperature, such that the same conclusion in favor of the heat bath regulator applies.Turn- 4 It is worth mentioning here that for higher values of α than the ones listed in Fig. 8 the stability of the flow decreased significantly.
ing to the classical limit of the FRG, we see that the resulting spectral function lies significantly closer to the classical-statistical result, most visible through the shape of the main peak, as expected.The shape of the second bump of the classical limit in the FRG also agrees with the classical-statistical spectral function.It is however located at lower frequencies, which is an effect of the structure of the non-local sunset diagrams in the flow equation of the 2-point function (cf. the discussion in Section VI in the context of Fig. 6), and therefore a purely truncational issue.We conclude from this discussion that a heat bath regulator produces more accurate results than simpler options like the CS+Damping regulator, and it therefore should be the preferred choice for real-time calculations.

Full flow equations
For the 4-point vertex, we choose to neglect the contributions arising from anomalous vertex on the r.h.s. of equation (5.47) entirely, and to approximate both the classical and the quantum vertex by combining all (possibly non-local interaction terms) into one effective local interaction.In terms of the effective coupling constants, this choice reads The flow equation for the 6-point function is gained in principle by functionally differentiating the flow equation (5.1) six times.Since this is rather cumbersome in general, as mentioned above, we instead use the method of Taylor expanding the flow of the effective potential to sixth order as described in Subsection V C 1.As a shorthand notation we define for the effective coupling constants.
t e x i t s h a 1 _ b a s e 6 4 = " 5 s G F E 6 E I g O g X y 8 0 2 z P j C 2 y F s R yE = " > A A A D f 3 i c d V N L b x M x E H Y a H i W 8 U j h y s U h Q k a D R b q D Q S B w q c e F Y J N J W i q P K 6 5 1 N r P i x s p 0 2 k e s / y Y 1 / g s Q F O 2 k l E m C k l U b f N y / P N 1 v U g l u X ZT 8 a O 8 0 7 d + / d 3 3 3 Q e v j o 8 Z O n 7 b 1 n p 1 b P D Y M h 0 0 K b 8 4 J a E F z B 0 H E n 4 L w 2 Q G U h 4 K y Y f U 7 8 2 S U Y y 7 X 6 5 p Y 1 j C W d K F 5 d 7 a Y U W o l D J S 4 c i 0 T a S n E U e b 2 z i R U / V r b T N n L 9 J 7 n x T 5 C 4 Y C e t R A O M t N L o + + b l + W b L R n D r 8 v x H a 6 N 9 7 / 6 D h 5 u P O o + f P H 2 2 t b 3 z / N T q u W E w Y F p o c 1 5 ) for a general heat-bath kernel K(ω) = |ξ(ω)| 2 β .Now we use the definition (3.10) of the classical-statistical spectral function to arrive at a balance-type equation

6 FIG. 2 .
FIG. 2. Real and imaginary parts of R Rk (ω) in units of k 2 for the bath with the analytic cutoff in Eq. (5.15).
≡ d D x . . ..(5.30)The first line in (5.29) corresponds to the 2-loop exact 2-point function, the second line to the 1-loop exact 4point function, and the third line to the '0-loop' exact 6-point function.Their detailed structures are explained in reversed order, starting from the 6-point function, in Subsections V C 1, V C 2, and V C 3, respectively.

FIG. 6 .
FIG.6.Strong-coupling spectral functions with λ = 4 (and γ = 0.06) at successively lower temperatures starting with T = 4 in (a) down to T = 0.5 in (d).With decreasing temperature the contributions from higher states in the one-step and three-step transitions get exponentially suppressed more and more so that the corresponding ensembles of main and second peak in the spectral functions get compressed towards their lowest transition frequencies until only a few individual transition lines remain.The (static) GSA spectral function follows the ensemble averages of main and second peak more closely than the classical one as temperature is lowered beyond the range of validity of the classical approximation (T ω).In our present truncation of the real-time FRG, the second peak for ω ≈ 3ωc stays put at about three times the central frequency ωc of the main peak (a frequency grid with 320 points on ω ∈ [0, 8] was used for T = 1, 2, 4, and one with 200 points on ω ∈ [0, 5] for T = 0.5).As in Fig.5, the corresponding NLO 2-loop perturbative results are included with dashed lines for comparison.

FIG. 8 .
FIG. 8. Comparison between different choices for the Callan-Symanzik counter-term coefficient α at the UV parameters λ = 4, γ = 0.06 at two temperaturs T = 4 (a) and T = 0.5 (b) for the heat bath regulator ('HB+CS'), the CS+Damping regulator, and results from the classical-statistical simulation.We also included the result for the classical limit in the realtime FRG, which arises by deleting the quantum φ c φ q φ q φ q vertex and replacing the quantum distribution function with its Rayleigh-Jeans counterpart[17].It agrees for T = 4 very well with the full quantum result from the FRG, indicating that there we already are close to the classical limit.For T = 0.5 it differs from the quantum result and matches the classical-statistical spectral function better, as expected.(For the FRG we use an N = 320, L = 8 grid for T = 4, and an N = 200, L = 5 grid for T = 0.5.) all others = 0 (B.4c)where the effective coupling constants are the same for all permutations of the CTP indices.Equipped with this choice of the effective coupling constants, it is straight-forward to draw the diagrams in equation(5.47)  to arrive at the flow equations listed in equations (B.6a -B.6c).
c. Adiabatic Approximation