## Spectral functions and critical dynamics of the O (4) model from classical-statistical lattice simulations

S¨ oren Schlichting , Dominik Smith , Lorenz von Smekal a b b

a

Fakult¨ at f¨r Physik, Universit¨t Bielefeld, 33615 Bielefeld, Germany u a b Institut f¨r Theoretische Physik, Justus-Liebig-Universit¨t, Heinrich-Bu u a GLYPH&lt;11&gt; -Ring 16, 35392 Gießen, Germany

## Abstract

We calculate spectral functions of the relativistic O (4) model from real-time lattice simulations in classical-statistical field theory. While in the low and high temperature phase of the model, the spectral functions of longitudinal ( GLYPH&lt;27&gt; ) and transverse ( GLYPH&lt;25&gt; ) modes are well described by relativistic quasi-particle peaks, we find a highly non-trivial behavior of the spectral functions in the cross over region, where additional structures appear. Similarly, we observe a significant broadening of the quasi-particle peaks, when the amount explicit O (4) symmetry breaking is reduced. We further demonstrate that in the vicinity of the O (4) critical point, the spectral functions develop an infrared power law associated with the critical dynamics, and comment on the extraction of the dynamical critical exponent z from our simulations.

Keywords: Spectral functions, O (4) model, classical-statistical simulations, critical dynamics

PACS:

## 1. Introduction

Besides static equilibrium properties, real-time correlation functions are of great interest in a wide range of physical settings, ranging from heavy-ion collisions to condensed-matter physics, as they carry important information about the dynamical properties of classical and quantum systems. Specifically, for equilibrium systems, the corresponding spectral functions contain information about the quasi-particle spectrum of a theory and can be used to reconstruct all real-time and Euclidean correlation functions in thermal equilibrium, via the fluctuation-dissipation relation. In addition, one may also extract transport properties such as the bulk viscosity [1], the life times of resonances or particle production rates from appropriate spectral functions. In the vicinity of second-order phase transitions, one can even use spectral functions to identify the dynamic universality class of a system [2].

Non-perturbative calculations of spectral functions are tremendously di GLYPH&lt;14&gt; cult. Lattice field-theory simulations o GLYPH&lt;11&gt; er a first-principles approach, but these are typically carried out in Euclidean space-time. Subsequently an analytic continuation to real time must be performed, which is an ill-posed numerical problem as it involves computing an inverse Laplace transform from a finite set of data points of finite accuracy. Different reconstruction schemes exist, such as Maximum Entropy Methods [3-5], the Backus-Gilbert method [6] the Schlessinger Point or Resonances-via-Pad´ e method [7, 8] or Tikhonov regularization [9] but each of these comes with its own set of limitations (see Ref. [10] for a comparison). One interesting alternative to such reconstructions of spectral functions from lattice data, as done for the present model in Ref. [11], is given by functional approaches such as n GLYPH&lt;0&gt; PI [12] and Functional Renor- malization Group (FRG) methods [13-19] or Dyson-Schwinger equations (DSE) [20, 21], which can be analytically continued or formulated directly in the real-frequency domain. However, such approaches necessarily require truncations of an infinite set of evolution equations or equations of motion for n -point correlation functions, and thus greatly benefit from additional insights into the structure and dynamics of excitations.

In this work, we use classical-statistical lattice simulations in real time to compute the single-particle spectral function of a scalar field theory. Since critical phenomena in quantum field theories are governed by classical dynamics, universal properties can be computed in a corresponding classical theory [22]. Likewise, spectral functions can be approximated by products of classical fields close to a second-order phase transition. This approach is based on the fluctuation-dissipation relation or Kubo-Martin-Schwinger periodicity condition [23, 24] and becomes exact as one approaches the critical point. In the past, this method has been successfully applied to a singlecomponent scalar field theory in 2 + 1 dimensions [25] and used to verify that this theory belongs to the dynamic universality class of relaxational models with conserved density (Model C) [2] according to the classification scheme of Hohenberg and Halperin [26].

Here we focus on the relativistic isovector Lorentz-scalar field theory with internal O (4) symmetry (' O (4) model') in 3 + 1 space-time dimensions, which also exhibits a second-order phase transition. Clearly, this model is of particular relevance as an e GLYPH&lt;11&gt; ective theory for low energy QCD; in particular the chiral phase transition of QCD for two degenerate light-quark flavours is believed to be in the same O (4) universality class [11, 27-29]. Other O N ( ) models are of interest in a QCD con-

text as well, such as e.g. the O (3) model in 1 + 1 dimensions, which exhibits instanton solutions, asymptotic freedom and a trace anomaly [30-33]. Central objective of our study is to calculate and analyze the features of real-time spectral functions in the O (4) model within the classical-statistical approach. Even though strictly speaking the classical-statistical approximation is only justifiable at very high temperatures or in the vicinity of the critical point, we will also explore the behavior away from criticality, where our results can still provide qualitative insights which may serve as a valuable input to the non-perturbative functional methods mentioned above.

Starting with a brief outline of the methodology and simulation setup in Secs. 2 and 3, we proceed to the extraction of the phase diagram and analysis of the static critical behavior of the O (4) scalar-field model in Sec. 4. Simulation results for real-time spectral functions and dynamic critical behavior are presented in Sec. 5, where we discuss the behavior of the spectral functions across a crossover transition and in the vicinity of the critical point. Our conclusions are provided in Sec. 6.

## 2. Spectral functions, fluctuation-dissipation theorem and classical-statistical approximation

Consider an arbitrary bosonic Heisenberg operator ˆ ( O t ; x ) in a quantum field theory described by the Hamiltonian ˆ . H The spectral function of this operator is defined via the commutator

$$\rho ( t - t ^ { \prime }, \mathbf x - \mathbf y, T ) \, = \, i \langle [ \hat { O } ( t, \mathbf x ), \hat { O } ^ { \dagger } ( t ^ { \prime }, \mathbf y ) ] _ { - } \rangle \,, \quad ( 1 ) \quad \text{where}$$

where the expectation value in thermal equilibrium is

$$\langle \hat { O } ( t, { \mathbf x } ) \rangle \, = \, \frac { 1 } { Z } \, \text{Tr} \left ( e ^ { - \hat { H } / T } \hat { O } ( t, { \mathbf x } ) \right ) \,, \ \ Z \, = \, \text{Tr} \, e ^ { - \hat { H } / T } \,. \ \ ( 2 ) \quad \text{princi} \\ \text{(see e.e.)}$$

Besides the spectral function, which characterizes the structure of possible excitations, we can also consider the statistical two-point function, which characterizes statistical fluctuations of the fields, and is defined in the quantum theory from the anticommutator:

$$F ( t - t ^ { \prime }, \mathbf x - \mathbf y, T ) & = \frac { 1 } { 2 } \langle [ \hat { O } ( t, \mathbf x ), \hat { O } ^ { \dagger } ( t ^ { \prime }, \mathbf y ) ] _ { + } \rangle - \langle \hat { O } ( t, \mathbf x ) \rangle \langle \hat { O } ^ { \dagger } ( t ^ { \prime }, \mathbf y ) \rangle \,. \\ \mathfrak z \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \colon$$

In thermal equilibrium the statistical fluctuations F ( ) are conGLYPH&lt;1&gt; nected to the spectral function GLYPH&lt;26&gt; ( ) by the fluctuation-dissipation GLYPH&lt;1&gt; relation or Kubo-Martin-Schwinger (KMS) condition [23, 24, 34], which follows from the imaginary-time periodicity of the Euclidean propagator and is stated in Fourier space as

$$F ( \omega, p, T ) \, = \, \left ( n _ { T } ( \omega ) + \frac { 1 } { 2 } \right ) \rho ( \omega, p, T ) \,. \quad \quad ( 4 ) \, \quad \ \ F _ { \text{cl} } ( t )$$

Here, nT ( ! ) is the Bose-Einstein distribution. We furthermore denote the magnitude of the spatial momentum as p GLYPH&lt;17&gt; j p j and define the Fourier transformations by

$$F ( \omega, p, T ) \ & = \quad \int d t \, d ^ { 3 } x \, e ^ { i ( \omega t - \mathbf p x ) } F ( t, \mathbf x, T ) \,, \quad \ \ ( 5 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \, \rho ( \omega, p, T ) \ & = \quad - i \int d t \, d ^ { 3 } x \, e ^ { i ( \omega t - \mathbf p x ) } \rho ( t, \mathbf x, T ) \,. \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

In the limit of small frequencies ! GLYPH&lt;28&gt; T (or high temperatures) the Bose Einstein distribution nT ( ! ) = 1 (exp( = != T ) GLYPH&lt;0&gt; 1) is well approximated by nT ( ! ) GLYPH&lt;25&gt; T =! , which is precisely the Rayleigh-Jeans distribution of the occupation-number in a classical-statistical bosonic field theory. Since the universal properties in the vicinity of a finite-temperature phase transition are governed by infrared field modes with ! GLYPH&lt;28&gt; T , it is exactly this limit which is relevant to the study of critical dynamics. In the absence of quantum anomalies critical phenomena at a finite temperature phase transition are therefore rigorously characterized by classical dynamics and we will argue in the following that (1) is approximated with increasing precision by a product of classical field variables computed in a corresponding classical-statistical theory as one approaches a critical point.

We further note that a classical-statistical description of the dynamics also becomes applicable when statistical fluctuations GLYPH&lt;24&gt; F dominate over quantum fluctuations GLYPH&lt;24&gt; GLYPH&lt;26&gt; , as the classicalstatistical approximation (CSA) can formally be seen as a leading order expansion in F GLYPH&lt;29&gt; GLYPH&lt;26&gt; , as discussed in detail in [35]. Based on this idea, the classical-statistical description has also been applied to the study of equilibrium spectral functions in the high-temperature regime of scalar field theories [22].

In the classical limit there are no commutators, so the spectral function is given by

$$\rho _ { \text{cl} } ( t - t ^ { \prime }, \mathbf x - \mathbf y, T ) \, = \, - \langle \{ O ( t, \mathbf x ), O ^ { * } ( t ^ { \prime }, \mathbf y ) \} \rangle _ { \text{cl} } \,, \quad ( 6 )$$

where fGLYPH&lt;1&gt; ; GLYPH&lt;1&gt;g denotes the Poisson bracket, the expectation value is now computed with respect to a classical-statistical ensemble and O t ( ; x ) becomes a functional of classical fields GLYPH&lt;30&gt; ( t ; x ) and their conjugate momenta GLYPH&lt;25&gt; ( t ; x ). Even though one could in principle compute the spectral function directly using Eq. (6) (see e.g. [36]), it turns out that handling the Poisson bracket is impractical and there is a more elegant way to calculate equilibrium spectral functions in classical-statistical field theory [22].

Exploiting the fluctuation dissipation relation for ! GLYPH&lt;28&gt; T in the classical-statistical theory, Eq. (4) is approximated by

$$F ( \omega, p, T ) \, \approx \, \frac { T } { \omega } \, \rho ( \omega, p, T ) \,,$$

which in the time domain can be expressed as

$$\rho ( t, p, T ) \, \approx \, - \frac { 1 } { T } \, \partial _ { t } F ( t, p, T ) \,.$$

By using Eqs. (7), (8), and the fact that in the classical limit the statistical two-point function becomes

$$F _ { \text{cl} } ( t - t ^ { \prime }, \mathbf x - \mathbf y, T ) & = \langle O ( t, \mathbf x ) O ^ { * } ( t ^ { \prime }, \mathbf y ) \rangle _ { \text{cl} } - \langle O ( t, \mathbf x ) \rangle _ { \text{cl} } \, \langle O ^ { * } ( t ^ { \prime }, \mathbf y ) \rangle _ { \text{cl} } \,, \\ \text{rmore} \quad &. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.$$

we can construct simple expressions for di GLYPH&lt;11&gt; erent spectral functions which make use only of products of field variables, are exact for classical-statistical theories but also describe the universal critical behavior of quantum field theories in the same universality class [2].

In this work, we numerically obtain the single particle spectral function (i.e. O t ( ; x ) GLYPH&lt;17&gt; GLYPH&lt;30&gt; ( t ; x )) in the momentum domain.

We consider a real scalar field theory, so O y ( t ; x ) = O t ( ; x ). Using Eq. (8), we can write for the spectral function in real time

$$\lambda & \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \\ & \rho _ { \text{cl} } ( t, p, T ) & ( 1 0 ) & \quad \begin{pmatrix} e ^ { - \beta H } \dot { \cdot } \\ \text{funct} \\ & = - \frac { 1 } { T } \int \mathrm d ^ { 3 } x \, e ^ { - i \mathrm p x } \, \partial _ { t } F _ { \text{cl} } ( t, \mathbf x, T ) & \quad \text{config} \end{pmatrix} \\ & = - \frac { 1 } { 2 T } \int \mathrm d ^ { 3 } x \, e ^ { - i \mathrm p x } \, \langle \, \pi ( t, \mathbf x ) \, \phi ( 0, \mathbf 0 ) - \phi ( t, \mathbf x ) \, \pi ( 0, \mathbf 0 ) \, \rangle _ { \text{cl} } \,, & \quad \partial _ { t } \\ \text{where} \, \pi ( t, \mathbf v ) \, - \, \Delta \, \Delta ( t, \mathbf v ) \, \, \text{ and was used that fact that the decision} \quad \text{where}$$

where GLYPH&lt;25&gt; ( t ; x ) = @ GLYPH&lt;30&gt; t ( t ; x ), and we used the fact that the disconnected part vanishes due to h GLYPH&lt;25&gt; ( t ; x ) i = 0. We will focus for simplicity on the p = 0 component, for which the spectral function is explicitly given by

$$\rho _ { \text{cl} } ( t - t ^ { \prime }, p = 0, T ) \, = \, - \frac { V } { 2 T } \, \left \langle \Pi ( t ) \Phi ( t ^ { \prime } ) - \Phi ( t ) \Pi ( t ^ { \prime } ) \right \rangle _ { \text{cl} } \,, \ \ ( 1 1 ) \quad \text{where}$$

with V = R d 3 x and

$$\Pi ( t ) = \frac { 1 } { V } \int d ^ { 3 } x \pi ( t, \mathbf x ) \,, \quad \Phi ( t ) = \frac { 1 } { V } \int d ^ { 3 } x \phi ( t, \mathbf x ) \,. \quad ( 1 2 ) \quad \text{The st} \\ \quad \text{the Eu}$$

Since in practice the spectral functions GLYPH&lt;26&gt; cl ( t GLYPH&lt;0&gt; t 0 ; 0 ; T ) are obtained directly in the time domain, it is then straightforward to obtain the corresponding spectral functions in the frequency domain by a Fourier transform.

In the above discussion, we have used a single-component scalar field for illustration. We note that for the O (4) model, the spectral function GLYPH&lt;26&gt; ab ( t ; x ) is computed individually for the di GLYPH&lt;11&gt; erent field components GLYPH&lt;30&gt; a ( t ; x ). By introducing an explicit symmetry breaking, the O (4) symmetry is broken down to O (3), and we can distinguish between the directions parallel and perpendicular to the vacuum alignment of the order parameter h GLYPH&lt;30&gt; a ( t ; x ) i , which we will refer to as the GLYPH&lt;27&gt; (parallel) and GLYPH&lt;25&gt; (perpendicular) components.

## 3. Simulation setup

We study the classical equilibrium properties of the 3 + 1 dimensional O (4) model defined by the lattice Hamiltonian

$$H = \sum _ { i } a _ { s } ^ { d } \left \{ \frac { 1 } { 2 } \pi _ { i } ^ { a } \pi _ { i } ^ { a } - \frac { 1 } { 2 a _ { s } ^ { 2 } } \sum _ { j \sim i } \phi _ { i } ^ { a } \phi _ { j } ^ { a } + \left ( \frac { m ^ { 2 } } { 2 } + \frac { d } { a _ { s } ^ { 2 } } \right ) \phi _ { i } ^ { a } \phi _ { i } ^ { a } & \quad \text{ then } \text{$i$ to a tiu} \\ + \frac { \lambda } { 4! N } \left ( \phi _ { i } ^ { a } \phi _ { i } ^ { a } \right ) ^ { 2 } + J ^ { a } \phi _ { i } ^ { a } \right \} \,, \ \ ( 1 3 ) & \quad \text{stices.} \\ \text{where $d^{a}$ are real valued field variables associated with the cities}$$

where GLYPH&lt;30&gt; a i are real valued field variables associated with the sites of a cubic lattice with periodic boundary conditions, GLYPH&lt;25&gt; a i are conjugate momenta, a = f 1 ; : : : 4 g label the components of the fields, J a = GLYPH&lt;14&gt; a 1 J denotes an explicit symmetry breaking term, d is the number of spatial dimensions (we consider the case d = 3) and as denotes the spatial lattice spacing. We set as = 1 in the following, which implies that all dimensionful quantities are understood to be expressed in units of as from here on. The sum P j GLYPH&lt;24&gt; i runs over all nearest neighbors j of site i .

In classical thermal equilibrium with inverse temperature GLYPH&lt;12&gt; = 1 = T , the expectation value of a static observable O [ GLYPH&lt;30&gt;; GLYPH&lt;25&gt; ] is defined as

$$\langle O \rangle = \, \int D \phi D \pi \, O [ \phi, \pi ] \, e ^ { - \beta H } \,, \, \quad \quad ( 1 4 )$$

and can be computed in a straightforward way by generating an ensemble of classical field configurations with distribution e GLYPH&lt;0&gt; GLYPH&lt;12&gt; H and subsequently evaluating the observable O [ GLYPH&lt;30&gt;; GLYPH&lt;25&gt; ] as a function of the fundamental fields. In practice we generate our configurations using a Langevin prescription

$$\partial _ { t _ { L } } \phi _ { i } ^ { a } = \frac { \partial H } { \partial \pi _ { i } ^ { a } } \, \ \ \partial _ { t _ { L } } \pi _ { i } ^ { a } = - \frac { \partial H } { \partial \phi _ { i } ^ { a } } - \gamma \pi _ { i } ^ { a } + \sqrt { 2 \gamma T } \, \xi _ { i } ^ { a } \, \quad ( 1 5 )$$

where t L denotes the Langevin time, GLYPH&lt;24&gt; a i corresponds to a Gaussian white noise h GLYPH&lt;24&gt; a i GLYPH&lt;24&gt; b j i = GLYPH&lt;14&gt; ab GLYPH&lt;14&gt; i j and

$$\frac { \partial H } { \partial \pi _ { i } ^ { a } } = \pi _ { i } ^ { a } \, \quad \frac { \partial H } { \partial \phi _ { i } ^ { a } } = - ( \Delta \phi ^ { a } ) _ { i } + \left ( m ^ { 2 } + \frac { \lambda } { 6 N } \phi _ { i } ^ { b } \phi _ { i } ^ { b } \right ) \phi _ { i } ^ { a } + \delta ^ { a 1 } J \, \ ( 1 6 )$$

where the Laplacian is discretized as

$$( \Delta \phi ^ { a } ) _ { i } = \sum _ { j \sim i } ( \phi _ { j } ^ { a } - \phi _ { i } ^ { a } ) \.$$

The stochastic di GLYPH&lt;11&gt; erential equation is solved numerically using the Euler-Maruyama scheme with an update step at = as = 0 01 : and - if not stated otherwise - we employ the set of parameters m 2 = GLYPH&lt;0&gt; 1, GLYPH&lt;21&gt; = 1 and GLYPH&lt;13&gt; = 0 3. : Note that in order to assess the universal critical behavior of the model, the coupling constant GLYPH&lt;21&gt; can be tuned to an optimal value to reduce scaling corrections [37]. However, we did not pursue this in our study.

Besides the static observables it is also straightforward to compute unequal time correlation functions in the classicalstatistical field theory. This allows for a simple prescription to extract the classical-statistical spectral function in real time through Eqs. (11) and (12). Since the classical fields GLYPH&lt;30&gt;; GLYPH&lt;25&gt; obey Hamilton's equations of motion

$$\partial _ { t } \phi _ { i } ^ { a } = \frac { \partial H } { \partial \pi _ { i } ^ { a } } \, \quad \partial _ { t } \pi _ { i } ^ { a } = - \frac { \partial H } { \partial \phi _ { i } ^ { a } } \, \quad \quad ( 1 8 )$$

it is straightforward to compute the unequal time correlation function h P i GLYPH&lt;30&gt; a i ( ) t P j GLYPH&lt;25&gt; a j ( t 0 ) i entering Eqs. (11) and (12). We first generate an ensemble of initial field configurations, and then independently evolve the classical field configurations up to a time max t t ( ; 0 ) based on a leap-frog scheme with at = as = 0 05 unless stated otherwise. : By saving the evolution of the order-parameter field along the classical trajectories, we subsequently extract the correlation functions between di GLYPH&lt;11&gt; erent time slices.

## 4. Results: Static universality

Critical exponents and scaling functions of the threedimensional O (4) spin model have been studied extensively using lattice simulations [38-43]. Before we discuss our results for real-time spectral functions, we verify that we reproduce the expected static critical properties and extract the phase diagram of our field theoretical model (13) in the J GLYPH&lt;0&gt; T plane. Our basic observables for this purpose are cumulants of the ferromagnetic order parameter, which in the presence of an explicit symmetry breaking is defined as

$$\phi _ { J } = \frac { 1 } { | J | V } \left ( \sum _ { i } J _ { a } \phi _ { i } ^ { a } \right ) \,.$$

Conversely, in the absence of an explicit symmetry breaking, we employ

$$\Lambda _ { 1 } ^ { \prime } \left | \phi \right | = \frac { 1 } { V } \sqrt { \left ( \sum _ { i } \phi _ { i } ^ { a } \right ) \left ( \sum _ { j } \phi _ { j } ^ { a } \right ) } \, \quad \ \ \ ( 2 0 ) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

as a proxy for the order parameter. In the following we will also use the symbol GLYPH&lt;30&gt; to generically refer to either j GLYPH&lt;30&gt; j or GLYPH&lt;30&gt; J , when relations of the same form apply to both equally.

In the vicinity of the critical point, the leading dependence of GLYPH&lt;30&gt; on the reduced temperature Tr = T GLYPH&lt;0&gt; Tc Tc , the explicit symmetry breaking J and the linear system size L (we always consider 3dlattices with Lx = Ly = Lz ), follows from the scaling relation

$$\phi ( T, J, L ^ { - 1 } ) = s ^ { \beta } \phi _ { 0 } \, \Phi \left ( s ^ { - 1 } T _ { r }, s ^ { - \beta \delta } \frac { J } { J _ { 0 } }, s ^ { - \nu } \left ( \frac { L } { L _ { 0 } } \right ) ^ { - 1 } \right ) \,, \quad ( 2 1 ) \quad \text{$\frac{\nolimits}{\text{within}$} } \\ \text{$\text{$\text{$two$ the} } }$$

with a universal scaling function GLYPH&lt;8&gt; and non-universal amplitudes GLYPH&lt;30&gt; 0 , J 0 , L 0 . By adapting the normalization conditions GLYPH&lt;8&gt; (1 ; 0 0) ; = GLYPH&lt;8&gt; (0 ; 1 0) ; = GLYPH&lt;8&gt; (0 ; 0 1) ; = 1, the critical behavior of the order parameter is determined by 1

$$\phi ( T _ { r }, 0, 0 ) \ & = \ & | T _ { r } | ^ { \beta } \, \phi _ { 0 } \, & \text{ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \phi ( 0, J, 0 ) \ & = \ & | J | ^ { 1 / \delta } \, \phi _ { 0 } J _ { 0 } ^ { - 1 / \delta } \, & \text{ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \\ \phi ( 0, 0, L ) \ & = \ & | L | ^ { - \beta / \nu } \, \phi _ { 0 } L _ { 0 } ^ { \beta / \nu } \, & \text{ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \\ ( 2 2 ) \quad \text{ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

and finite-size scaling relations take the form

$$\phi ( T _ { r }, 0, L ) = \left ( \frac { L } { L _ { 0 } } \right ) ^ { - \beta / \nu } f _ { T } \left ( \left ( \frac { L } { L _ { 0 } } \right ) _ { \lambda } ^ { 1 / \nu } T _ { r } \right ) \,, \quad \ \ ( 2 3 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

$$\phi ( 0, J, L ) = \left ( \frac { L } { L _ { 0 } } \right ) ^ { - \beta / \nu } f _ { J } \left ( \left ( \frac { L } { L _ { 0 } } \right ) ^ { \frac { \beta } { \nu / \delta } } \frac { J } { J _ { 0 } } \right ) \,. \quad \ ( 2 4 ) \quad \chi _ { | \phi | } ( \cdot )$$

We will also consider the static susceptibilities

$$\chi _ { | \phi | } = V ( \langle | \phi | ^ { 2 } \rangle - \langle | \phi | \rangle ^ { 2 } ) \,, \ \chi _ { J } = V ( \langle \phi _ { J } ^ { 2 } \rangle - \langle \phi _ { J } \rangle ^ { 2 } ) \,, \quad ( 2 5 ) \quad \text{$\ e p e c $}$$

and the Binder cumulant of the order parameter, which is given by

$$\chi _ { 4 } = 1 - \frac { \langle | \phi | ^ { 4 } \rangle } { 3 \langle | \phi | ^ { 2 } \rangle ^ { 2 } } \,. \quad \quad \ ( 2 6 ) \quad \ \ \ \ L \, \cdot \, \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

## 4.1. Static universality - T dependence

Webegin with studying the temperature dependence at J = 0, which is summarized in Figs. 1, 2. Individual points in each figure correspond to the data obtained from simulations at the corresponding temperature values, while solid bands are obtained by performing a multi-histogram re-weighting analysis [44], using the data from the closest six temperature points. Errorbars are obtained from a jacknife analysis.

We first estimate the the critical temperature Tc at J = 0. Since the critical exponents GLYPH&lt;12&gt; = 0 380(2) and : GLYPH&lt;23&gt; = 0 7377(41) :

1 We have obtained the following values GLYPH&lt;30&gt; 0 = 4 5 : GLYPH&lt;6&gt; 0 5, : GLYPH&lt;30&gt; 0 J GLYPH&lt;0&gt; 1 =GLYPH&lt;14&gt; 0 = 2 4 : GLYPH&lt;6&gt; 0 2, : GLYPH&lt;30&gt; 0 L GLYPH&lt;12&gt;=GLYPH&lt;23&gt; 0 = 4 5 : GLYPH&lt;6&gt; 0 3 as rough estimates for the non-universal amplitudes in our : model.

have been determined very precisely from spin-model simulations [43], we use these results and exploit the third identity in Eq. (22) along with the finite-size scaling relation (23) to infer the critical temperature from the order parameter. We find by plotting L GLYPH&lt;12&gt;=GLYPH&lt;23&gt; hj GLYPH&lt;30&gt; ji for di GLYPH&lt;11&gt; erent lattice sizes L = 48 64 96 128 ; ; ; as a function of T that all curves intersect in a single point with good accuracy (Fig. 1, left), which then determines the critical temperature Tc = 17 3925(10). Subsequently, we explicitly : verify the universal finite-size scaling of our data, by plotting the same observable L GLYPH&lt;12&gt;=GLYPH&lt;23&gt; hj GLYPH&lt;30&gt; ji as a function of the rescaled reduced temperature L 1 =GLYPH&lt;23&gt; Tr (Fig. 1, right). All data points collapse onto a single universal scaling curve, indicating that for typical ranges of Tr and lattice sizes L our simulations are well within the scaling window.

Even though the O (4) universality class is strongly favored by these consistency checks, we find that our data do not constrain the critical exponents at the same level of accuracy as in the spin models. By optimizing the scaling collapse across di GLYPH&lt;11&gt; erent data sets we can, for instance, obtain the estimates Tc GLYPH&lt;25&gt; 17 395 : GLYPH&lt;6&gt; 0 02, : GLYPH&lt;12&gt;=GLYPH&lt;23&gt; GLYPH&lt;25&gt; 0 53 : GLYPH&lt;6&gt; 0 015 and 1 : =GLYPH&lt;23&gt; GLYPH&lt;25&gt; 1 38 : GLYPH&lt;6&gt; 0 06, : which are consistent with the determination of Tc above and the critical-exponent values from [43]. Here we have estimated the errors, which are always dominated by the systematic uncertainties, by sequentially excluding di GLYPH&lt;11&gt; erent lattice sizes from our analysis. We have also checked that the value of the critical Binder cumulant GLYPH&lt;31&gt; 4 ( Tc ) = 0 63 : GLYPH&lt;6&gt; 0 01 agrees well with the : values reported in the literature [40, 45].

Next, we turn to the susceptibility GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j for which the critical behavior at J = 0 is determined by

$$\chi _ { | \phi | } ( T, L = \infty ) \sim | T _ { r } | ^ { - \gamma }, \quad \chi _ { | \phi | } ( T, L ) = \left ( \frac { L } { L _ { 0 } } \right ) ^ { \gamma / \nu } g _ { T } \left ( \left ( \frac { L } { L _ { 0 } } \right ) ^ { 1 / \nu } T _ { r } \right ). \\ \mathbb { m } { \colon } \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \colon \end{array}$$

This relation is obtained from Eq. (21) by di GLYPH&lt;11&gt; erentiating with respect to J and using the hyperscaling relation GLYPH&lt;13&gt; = GLYPH&lt;12&gt; ( GLYPH&lt;14&gt; GLYPH&lt;0&gt; 1). We first study the T dependence of GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j for L = 48 64 96 128 ; ; ; (Fig. 2, left) and verify that the pseudo-critical transition temperature Tpc ( L ), corresponding to the position of the peak, moves towards our estimate of Tc with increasing system size L ! 1 . We also confirm the finite-size scaling law (27), by plotting L GLYPH&lt;0&gt; GLYPH&lt;13&gt;=GLYPH&lt;23&gt; GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j as a function of L 1 =GLYPH&lt;23&gt; Tr (with GLYPH&lt;13&gt; = 1 4531(104) : taken from [43]) and verifying that the results collapse onto a single curve (Fig. 2, right). While for temperatures T &gt; Tpc ( L ) (above the pseudo-critical transition temperature) we find good agreement between di GLYPH&lt;11&gt; erent data sets, such scaling breaks down below the pseudo-critical temperature. Since for T &lt; Tpc ( L ) the susceptibility receives additional contributions of massless Goldstone modes, one expects to find a linear scaling of the susceptibility with the volume, which has been discussed in detail in [39] and is confirmed by our data.

## 4.2. Static universality - J dependence

So far we have verified the static critical behavior in the absence of explicit symmetry breaking, using the absolute value j GLYPH&lt;30&gt; j as an approximate order parameter. We now proceed by setting the temperature T to Tc GLYPH&lt;25&gt; 17 3925 and study the depen-: dence on the external field J . Again, we first consider the order

Figure 1: Temperature dependence and critical scaling of ferromagnetic order parameter at J = 0 for L = 48 64 96 128 ( ; ; ; GLYPH&lt;12&gt; and GLYPH&lt;23&gt; taken from [43]). Left: Tc = 17 3925(10) is determined by the intersection of : L GLYPH&lt;12&gt;=GLYPH&lt;23&gt; hj GLYPH&lt;30&gt; ji for di GLYPH&lt;11&gt; erent L . Right: Collapse of data points of di GLYPH&lt;11&gt; erent L onto a universal scaling function.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-ac11fd5f-f4d6-465c-8527-790f1b94bc45.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-95ea76b0-e02e-4ab6-b669-47b6dbd6e022.jpg)

Figure 2: Temperature dependence (left) and critical scaling (right) of ferromagnetic susceptibility at J = 0 for L = 48 64 96 128 ( ; ; ; GLYPH&lt;13&gt; and GLYPH&lt;23&gt; taken from [43]). Tc as obtained from finite-size scaling of hj GLYPH&lt;30&gt; ji is marked by a vertical line (deviations of the peak of GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j from this line are due to finite-size e GLYPH&lt;11&gt; ects).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-65484f12-ee6c-4ce2-a8cf-10c3e16f5717.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-09b4a7ac-6ced-4c96-a7b0-8adbd625fcd9.jpg)

φ

Figure 3: J dependence (left) and magnetic scaling (right) of order parameter GLYPH&lt;30&gt; at T = Tc . Absolute value hj GLYPH&lt;30&gt; ji and component GLYPH&lt;30&gt; J in direction of the external field h GLYPH&lt;30&gt; J i are shown. For small L GLYPH&lt;12&gt;GLYPH&lt;14&gt;=GLYPH&lt;23&gt; J the magnetic scaling of j GLYPH&lt;30&gt; j and GLYPH&lt;30&gt; J di GLYPH&lt;11&gt; ers. As J decreases, L must be increased for hj GLYPH&lt;30&gt; ji = h GLYPH&lt;30&gt; J i to hold. Solid line in left panel shows GLYPH&lt;30&gt; = cJ 1 =GLYPH&lt;14&gt; with GLYPH&lt;14&gt; = 4 824 taken from [43] and : c chosen such that the points for the largest lattices at small J are traversed. Solid line in right panel shows universal finite size scaling function for GLYPH&lt;30&gt; J as obtained in [43] from O (4) spin model simulations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-079453b4-58ea-4a21-93df-82ad02bddc4d.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-639d6065-7c12-46de-a57b-ac039939c4dd.jpg)

Figure 4: Magnetic scaling of susceptiblity GLYPH&lt;31&gt; . (left) Susceptibilities of absolute value j GLYPH&lt;30&gt; j and component GLYPH&lt;30&gt; J in direction of external field di GLYPH&lt;11&gt; er unless L GLYPH&lt;14&gt;GLYPH&lt;12&gt;=GLYPH&lt;23&gt; J is su GLYPH&lt;14&gt; ciently large. (right) Longitudinal and transverse susceptibilities GLYPH&lt;31&gt; GLYPH&lt;25&gt;=GLYPH&lt;27&gt; exhibit universal scaling at large L GLYPH&lt;14&gt;GLYPH&lt;12&gt;=GLYPH&lt;23&gt; J , but become indistinguishable at small L GLYPH&lt;14&gt;GLYPH&lt;12&gt;=GLYPH&lt;23&gt; J . Solid lines in both panels show universal finite size scaling function for GLYPH&lt;31&gt; GLYPH&lt;27&gt; as obtained in [43] from O (4) spin model simulations.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-a9dad119-3bab-48e8-a7dd-8da2ceae1c6b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-1b2b1f01-0353-4ff2-a73e-bc243cb0c8ab.jpg)

parameter itself and verify the power-law behavior GLYPH&lt;30&gt; GLYPH&lt;24&gt; j J j 1 =GLYPH&lt;14&gt; and the magnetic scaling (24), taking GLYPH&lt;14&gt; = 4 824(9) from [43]. : Fig. 3 summarizes these results. In the presence of an explicit symmetry breaking we can distinguish between the absolute value hj GLYPH&lt;30&gt; ji and the component h GLYPH&lt;30&gt; J i in direction of the external field. In principle both exhibit identical universal properties, but as J becomes smaller the quantities di GLYPH&lt;11&gt; er unless the system size L is simultaneously increased by a su GLYPH&lt;14&gt; cient amount. What is striking is that both quantities independently show magnetic scaling, whereby the data points for L GLYPH&lt;12&gt;=GLYPH&lt;23&gt; GLYPH&lt;30&gt; collapse onto single but distinct curves (Fig. 3, right).

Similarly, we also verify the J dependent critical properties of the susceptibility, which are given by

$$\chi ( J, L = \infty ) \sim | J | ^ { - \frac { \gamma } { \beta } }, \quad \chi ( J, L ) = \left ( \frac { L } { L _ { 0 } } \right ) ^ { \gamma / \nu } g _ { J } \left ( \left ( \frac { L } { L _ { 0 } } \right ) ^ { \frac { \beta \delta } { \nu } } \frac { J } { J _ { 0 } } \right ), \quad \text{cases} \quad \text{safely} \\ \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \end{cases}$$

where the symbol GLYPH&lt;31&gt; is a generic placeholder for GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j and GLYPH&lt;31&gt; J . Just as the order parameters j GLYPH&lt;30&gt; j and GLYPH&lt;30&gt; J , the corresponding susceptibilities GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j and GLYPH&lt;31&gt; J can be distinguished and independently collapse onto distinct critical scaling functions, which are di GLYPH&lt;11&gt; erent in the finite size scaling regime (for small L GLYPH&lt;12&gt;GLYPH&lt;14&gt;=GLYPH&lt;23&gt; J ) but merge for su GLYPH&lt;14&gt; ciently large L GLYPH&lt;12&gt;GLYPH&lt;14&gt;=GLYPH&lt;23&gt; J (see Fig. 4, left).

In the presence of a non-zero explicit symmetry breaking term J , we can also distinguish between the longitudinal ( GLYPH&lt;27&gt; ) and transverse ( GLYPH&lt;25&gt; ) components GLYPH&lt;31&gt; GLYPH&lt;27&gt;=GLYPH&lt;25&gt; of the susceptibility, where GLYPH&lt;31&gt; GLYPH&lt;27&gt; = GLYPH&lt;31&gt; J (cf. Eq. (25)) and GLYPH&lt;31&gt; GLYPH&lt;25&gt; is given by

$$\chi _ { \pi } = \frac { V } { N - 1 } \left [ \langle | \phi | ^ { 2 } \rangle - \langle \phi _ { J } ^ { 2 } \rangle \right ] \,. \quad \quad ( 2 9 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Independent finite size scaling of longitudinal and transverse susceptibilities is again observed (Fig. 4, right). For su GLYPH&lt;14&gt; ciently large values of L GLYPH&lt;12&gt;GLYPH&lt;14&gt;=GLYPH&lt;23&gt; J , i.e. close to the infinite volume limit, both curves are expected to approach the scaling behavior in Eq. (28), with a universal amplitude ratio GLYPH&lt;31&gt; GLYPH&lt;27&gt; =GLYPH&lt;31&gt; GLYPH&lt;25&gt; = 1 =GLYPH&lt;14&gt; (see [42]). This is nicely confirmed by our data. Even though most of our data points are outside the infinite volume scaling regime ( L GLYPH&lt;12&gt;GLYPH&lt;14&gt;=GLYPH&lt;23&gt; J !1 ) where GLYPH&lt;31&gt; exhibits a power law dependence

(cf. Eq. (28)), we also observe that the finite size scaling regime extends to much smaller values of L GLYPH&lt;12&gt;GLYPH&lt;14&gt;=GLYPH&lt;23&gt; J . In particular, for very small values of L GLYPH&lt;12&gt;GLYPH&lt;14&gt;=GLYPH&lt;23&gt; J the two scaling curves become almost indistinguishable, as the distinction between longitudinal and transverse components becomes less and less meaningful.

When comparing our results for GLYPH&lt;30&gt; J and GLYPH&lt;31&gt; GLYPH&lt;27&gt; to the universal finite size scaling functions determined in [43] for the O (4) spin model (displayed as solid lines in Fig. 4 and the right panel of Fig. 3), good agreement is found across the entire range where the parametrization is available.

## 4.3. Static universality - Conclusion

We conclude from all of the above that both, the T and J dependent static critical properties are indeed correctly reproduced in our classical statistical simulations, and we can now safely proceed to study real-time properties. In order to set the stage for our study of real-time correlation functions, we finally sketch the phase diagram in the J T -plane. For this purpose we compute the T dependence of the susceptibility GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j on a L = 64 lattice for several di GLYPH&lt;11&gt; erent values of J and carry out an interpolation for regions in between data points. For each line of constant J we then determine the maximum of GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j , which serves as an estimate for the pseudo-critical temperature along this line, and the inflection points. Our results for the phasediagram are shown in Fig. 5, where the color coding indicates the magnitude of the susceptibility. Most importantly, the horizontal and vertical dashed green lines correspond to the values of T J ; considered in our study of spectral functions.

## 5. Results: Spectral functions

We now study spectral functions GLYPH&lt;26&gt; ( t ; t 0 ) which we can directly extract as a function of the real-time variables t ; t 0 according to the procedure discussed in Secs. 2, 3. We prepare Nconf independent initial configurations for each choice of parameters and evolve each of them up to maximum time t Max (typically GLYPH&lt;24&gt; 10 4 GLYPH&lt;0&gt; 10 5 as ), recording the evolution of P i GLYPH&lt;30&gt; a i and P i GLYPH&lt;25&gt; a i . Based on this data, we then construct the un-equal time

Figure 5: Left: Temperature dependence of magnetic susceptibility GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j for di GLYPH&lt;11&gt; erent values of J on a L = 64 lattice. Pseudocritical temperatures are estimated from location of maximum. Right: Tentative phase diagram of 3 + 1 dimensional O (4) model in J GLYPH&lt;0&gt; T plane. Blue shading denotes the interpolated magnitude of the susceptibility GLYPH&lt;31&gt; j GLYPH&lt;30&gt; j at each point while dashed black and grey lines mark its maximum and inflection points respectively. Computation of spectral functions were carried out along dashed green lines.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-f865fb9b-dd48-4350-96e2-a86a95961c5f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-150dbe1f-53ad-4e0f-a35b-f53855cb55a2.jpg)

correlation function as a function of t GLYPH&lt;0&gt; t 0 according to the righthand side of Eq. (11), evaluated separately for each configuration, while immediately averaging over di GLYPH&lt;11&gt; erent positions t + t 0 in the real-time evolution. Statistical averages and errors of the spectral function GLYPH&lt;26&gt; ( t GLYPH&lt;0&gt; t 0 ) and its Fourier transform GLYPH&lt;26&gt; ( ! ) are computed from averaging over the ensemble of typically Nconf = 32 independent configurations. If not stated otherwise, all results have been obtained on L = 64 lattices, and we have checked at the example of a few data points that except for the immediate vicinity of the critical point our results remain unchanged when going to larger lattices. We further emphasize that in order to properly distinguish the longitudinal ( GLYPH&lt;25&gt; ) and transverse ( GLYPH&lt;27&gt; ) field components in our simulations, we always have to introduce a non-zero explicit symmetry breaking J .

## 5.1. Spectral functions in the crossover regime - T dependence

Before we turn our attention to the behavior in the vicinity of the critical point, we first present results for the temperature dependence of the spectral function at a relatively large explicit symmetry breaking J = 0 05 as indicated by the vertical dashed : line in the phase-diagram (cf. Fig. 5). We note that at such large values of J the transition is a relatively smooth cross over, with a pseudo-critical transition temperature Tpc ( J = 0 05) : GLYPH&lt;25&gt; 19 5. : Wehave collected all of our results for the GLYPH&lt;27&gt; - and GLYPH&lt;25&gt; -components of the spectral functions in Figs. 6 and 7. Di GLYPH&lt;11&gt; erent rows in Figs. 6 and 7 show the spectral functions in di GLYPH&lt;11&gt; erent temperature regimes, starting from very low temperatures in the top row, to temperatures just below the pseudo-critical transition temperature Tpc in the middle row, all the way to temperatures above the pseudo-critical transition temperature Tpc in the bottom row. Di GLYPH&lt;11&gt; erent columns in Figs. 6 and 7, all show the same data for the spectral functions but plotted in di GLYPH&lt;11&gt; erent ways in order to highlight the various features more clearly.

Before we turn to our simulation results we note that in the limit T ! 0 thermal fluctuations are suppressed, and the classical-statistical description reduces to the mean-field limit, where

$$\bar { \xi }, \quad \rho _ { \pi / \sigma } ( \omega, T = 0, p = 0 ) = 2 \pi \, \text{sign} ( \omega ) \, \delta ( \omega ^ { 2 } - m _ { \pi / \sigma } ^ { 2 } ) \, \quad ( 3 0 )$$

and for GLYPH&lt;21&gt; J 2 GLYPH&lt;28&gt; GLYPH&lt;0&gt; m 6 (with m 2 &lt; 0) one has

$$m _ { \sigma } ^ { 2 } \simeq - 2 m ^ { 2 } + \sqrt { \frac { 3 } { 2 N } \frac { \lambda J ^ { 2 } } { - m ^ { 2 } } } \,, \text{ and } m _ { \pi } ^ { 2 } \simeq \sqrt { \frac { 1 } { 6 N } \frac { \lambda J ^ { 2 } } { - m ^ { 2 } } } \,.$$

With our parameters, for comparison with the corresponding peaks at the lowest temperature ( T = 1 0875) in Figs. 6 and 7, : this then amounts to m GLYPH&lt;27&gt; ' 1 425 and : m GLYPH&lt;25&gt; ' 0 101. :

Clearly, however, this trivial behavior is a result of the classical-statistical approximation, which misses all quantum e GLYPH&lt;11&gt; ects that would otherwise become important in the lowtemperature regime T GLYPH&lt;28&gt; Tpc , such as for instance the decay process GLYPH&lt;27&gt; ! 2 GLYPH&lt;25&gt; which would result in strong modifications of the vacuum spectral functions. Nevertheless, even though an entirely classical description is thus not particularly well suited for the description of the low-temperature physics, it is still interesting to investigate the behavior of the classical-statistical spectral functions in this regime.

Based on our simulation results in the low-temperature regime T &lt; 15 shown in the top row of Figs. 6 and 7, we find that at low temperatures the classical-statistical spectral functions, exhibit the expected quasi-particle behavior where towards the lowest temperature T = 1 0875 the masses of : GLYPH&lt;27&gt; and GLYPH&lt;25&gt; are already nicely seen to approach the mean-field estimates below Eq. (30), up to minute shifts and some collisional broadening due to the small but finite residual temperature. Besides these quasi-particle peaks, the GLYPH&lt;25&gt; spectral function also shows an additional cusp at higher frequencies which at very low temperatures occurs approximately for frequencies ! GLYPH&lt;24&gt; 3 m GLYPH&lt;25&gt; and should be attributed to a multi-pion excitation. By further increasing the temperature, the quasi-particle peaks remain, however the mass of GLYPH&lt;27&gt; becomes lighter as the vacuum expectation value of the GLYPH&lt;27&gt; field decreases. Even though the width of GLYPH&lt;27&gt; spectral function also increases, it turns out that except for a small enhancement at low frequency, the spectral function of

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-4dc911a5-1cb4-4c9d-a43d-dab504155371.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-07aa18f3-26e7-4106-a57f-0edc752eabdf.jpg)

ω

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-648c1768-a5ca-4940-891a-c8460bfd167e.jpg)

ω

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-af135750-e201-4551-bc4e-5e347e72e1a6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-225e7126-352f-4481-985d-b139526e8623.jpg)

ω

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-5292d2ac-518e-46d8-95fe-0ed50a830ed2.jpg)

ω

Figure 6: T dependence of GLYPH&lt;27&gt; spectral functions for J = 0 05 below (top row), around (middle row) and above (bottom row) the pseudocritical temperature : Tpc GLYPH&lt;25&gt; 19 5. Left column: Raw data only. Middle column: Data with Breit-Wigner fits (31). Error bars are of the order of the pointsize. Several intermediate points : are not displayed for better visibility but are considered for the fits. Right column: Same as middle column but with log-scale.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-78e6ee0f-236c-4a08-a60a-d9e4f3c73a37.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-c405e84b-6f05-4421-880b-27b78ec48320.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-422fe63d-ef40-47fb-9ecb-2d8218108d81.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-53ef370b-f026-4ecc-8410-db53c28b85e5.jpg)

ω

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-e08281f4-965b-4458-b2ae-04683844ed2a.jpg)

ω

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-046f1bf0-ccaf-41cc-857f-86710f7ba229.jpg)

ω

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-5c7e3b06-1e61-4216-aa6a-d41fa7f933ee.jpg)

ω

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-84fd696f-3f7f-435f-b8a2-dc11ea97b6cc.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-8482bfa5-46cf-419a-a3f7-f784a9eccb29.jpg)

ω

Figure 7: T dependence of GLYPH&lt;25&gt; spectral functions for J = 0 05 below (top row), around (middle row) and above (bottom row) the pseudocritical temperature : Tpc GLYPH&lt;25&gt; 19 5. Left column: Raw data only. Middle column: Data with Breit-Wigner fits (31) to the larger peak, where applicable. Error bars are of the order of the : pointsize. Several intermediate points are not displayed for better visibility but are considered for the fits. Right column: Same as middle column but with log scale.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-d22754be-6435-45ae-aa69-caec18b51707.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-947392c2-05a3-4a94-8ea3-23ed37e50e7f.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-f0108e4b-e8c3-459b-8689-e204df2828e9.jpg)

Figure 8: Spectral function GLYPH&lt;26&gt; GLYPH&lt;25&gt; ( ! ) for J = 0 05 in the intermediate temperature range where two peaks are visible. : Solid lines represent fits to a Breit-Wigner distribution with two resonance (32).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-389cc332-236a-4df0-ada1-4d45e729640e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-4c23e4aa-0539-490f-b12d-499bc9efc7a1.jpg)

the GLYPH&lt;27&gt; mode can still be well described in terms of a single BreitWigner resonance

$$\rho ( \omega ) = \frac { \omega \Gamma } { ( \omega ^ { 2 } - m ^ { 2 } ) ^ { 2 } + \omega ^ { 2 } \Gamma ^ { 2 } } \quad \quad ( 3 1 ) \quad \text{$\text{$\text{$speed$}$}$} \\ \text{$bv$} \, t$$

as indicated by the solid lines, representing Breit-Wigner fits of the spectral function.

temperature Tpc = 19 5 but also further into the symmetry re-: stored phase. Despite the clear presence of a second peak, we find that for temperatures T &gt; 16 the dominant peak of the GLYPH&lt;25&gt; spectral function can again be described to reasonable accuracy by the Breit-Wigner distribution in Eq. (31), as indicated by the solid lines in Fig. 7.

Conversely, the spectral function for the GLYPH&lt;25&gt; mode exhibits a much more non-trivial behavior as the frequency threshold for the scattering states lowers and the resonance becomes more pronounced as temperature increases. Beyond T = 13 05 the : spectral function GLYPH&lt;26&gt; GLYPH&lt;25&gt; features an interesting double peak structure, where the e GLYPH&lt;11&gt; ective mass and spectral weight of the lower frequency peak decrease as a function of temperature, while the upper frequency peak becomes increasingly dominant when further increasing the temperature. We further illustrate this double peak structure in Fig. 8, which shows a close up of the GLYPH&lt;25&gt; spectral function in the same temperature regime. In order to track the widths and positions of the individual peaks, we also present fits to a double Breit-Wigner distribution, featuring two distinct resonances, of the form

$$\rho ( \omega ) = c _ { 1 } \frac { \omega \Gamma _ { 1 } } { ( \omega ^ { 2 } - m _ { 1 } ^ { 2 } ) ^ { 2 } + \omega ^ { 2 } \Gamma _ { 1 } ^ { 2 } } + c _ { 2 } \frac { \omega \Gamma _ { 2 } } { ( \omega ^ { 2 } - m _ { 2 } ^ { 2 } ) ^ { 2 } + \omega ^ { 2 } \Gamma _ { 2 } ^ { 2 } } \,. \ \ ( 3 2 ) \quad \sigma \text{spef} \quad \\ \text{sum} \quad.$$

While Eq. (32) provides a good description of the peaks, it tends to overestimate the spectral weight in the low frequency tails of the spectral function.

Spectral functions just below the crossover transitions at Tpc = 19 5, are presented in the middle rows of Figs. 6 and : 7 and show a smooth continuation of the temperature dependence observed at the lower temperatures. The spectral function of the GLYPH&lt;27&gt; mode continues to show a quasi-particle peak, where the e GLYPH&lt;11&gt; ective mass continues to becomes lighter but the width decreases again as expected when the pseudo-critical transition temperature is approached. While the double peak structure in the GLYPH&lt;25&gt; spectral function is most pronounced below the crossover temperature (i.e. around T ' 15) some remnants of the low frequency peak clearly persists not only up to the pseudo-critical

Beyond the pseudo-critical temperature T &gt; 19 5, the dom-: inant features of the spectral functions for GLYPH&lt;25&gt; and GLYPH&lt;27&gt; begin to coincide, as can be seen from comparing the results in the bottom rows of Figs. 6 and 7. One finds that as the temperature is increased further beyond Tpc , both GLYPH&lt;27&gt; and GLYPH&lt;25&gt; spectral functions are increasingly well described by the Breit-Wigner ansatz in Eq. (31), with increasing mass and decay width as a function of temperature. Between T = 19 5 and : T = 22 the additional low frequency peak in the GLYPH&lt;25&gt; spectral function slowly disappears, such that at the highest temperature T = 26 1 the spectral func-: tions for GLYPH&lt;27&gt; and GLYPH&lt;25&gt; become almost degenerate up to small di GLYPH&lt;11&gt; erences at very low frequencies ! glyph[lessorsimilar] 0 1, signaling the approxi-: mate restoration of the full O (4) symmetry on the level of the spectral functions.

Our results for the temperature dependence of the GLYPH&lt;25&gt; and GLYPH&lt;27&gt; spectral functions in the crossover regime, are compactly summarized in Fig. 9, where we show the temperature dependence of the e GLYPH&lt;11&gt; ective masses m GLYPH&lt;25&gt;=GLYPH&lt;27&gt; and decay width GLYPH&lt;0&gt; GLYPH&lt;25&gt;=GLYPH&lt;27&gt; obtained from the (single and double peak) Breit-Wigner fits. While at low temperatures GLYPH&lt;25&gt; and GLYPH&lt;27&gt; spectral functions in the classical-statistical approximation show well defined quasiparticle peaks, the GLYPH&lt;27&gt; mass rapidly decreases with increasing temperature, and the in-medium decay widths of GLYPH&lt;25&gt; and GLYPH&lt;27&gt; increase significantly below the pseudo-critical transition temperature Tpc . The scattering states in the GLYPH&lt;25&gt; spectral function lead to the development of an additional resonance peak as Tpc is approached. Interestingly, it is this emerging second peak which appears to develop further into the resonance peak that eventually becomes degenerate with the GLYPH&lt;27&gt; mode as temperature is increased further beyond Tpc . The original low-temperature quasi-particle peak on the other hand slowly melts and disap-

Figure 9: Dependence of (single-peak and double-peak) Breit-Wigner fit parameters for GLYPH&lt;25&gt; and GLYPH&lt;27&gt; spectral functions on the temperature T around the crossover transition at finite explicit symmetry breaking J = 0 05. :

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-31f98857-f728-48a0-ad16-70597707923e.jpg)

pears around the pseudo-critical temperature. The two distinct peaks at m GLYPH&lt;25&gt;; 1 and m GLYPH&lt;25&gt;; 2 around T ' 15 in fact show signs of an interesting avoided-crossing behavior which has not been observed in the corresponding solutions of analytically continued FRG flow equations, for example, so far.

## 5.2. Spectral functions in the crossover regime - J dependence

So far we have investigated the temperature dependence of the GLYPH&lt;27&gt; and GLYPH&lt;25&gt; spectral functions in the vicinity of the crossover transition, at a fixed relatively large explicit symmetry breaking J = 0 05. : Since we always have to keep a non-vanishinig explicit symmetry breaking in order to distinguish between GLYPH&lt;25&gt; and GLYPH&lt;27&gt; components, we will now fix the temperature close to Tc at T = 17 4 and decrease the explicit symmetry breaking : J by successive factors of two to approach the critical point, as indicated by the horizontal line in the phase diagram in Fig. 5.

Our results for the J dependence of spectral functions close to Tc are summarized in Fig. 10, where top and bottom rows show the GLYPH&lt;25&gt; and GLYPH&lt;27&gt; spectral functions at di GLYPH&lt;11&gt; erent explicit symmetry breaking. Starting from J = 0 05 employed in our tem-: perature scan of the crossover transition, we find that lowering the explicit symmetry breaking J results in a rapid decrease of the e GLYPH&lt;11&gt; ective mass of GLYPH&lt;27&gt; and GLYPH&lt;25&gt; along with a simultaneous increase of the decay width. E GLYPH&lt;11&gt; ectively the combination of these two phenomena leads to a melting of the quasi-particle peaks, in both GLYPH&lt;27&gt; and GLYPH&lt;25&gt; spectral functions, as can be seen from Fig. 11, where we present the J dependence of the Breit-Wigner resonance parameters. However, one should caution, that already at J = 0 05 : GLYPH&lt;2&gt; 2 GLYPH&lt;0&gt; 1 the GLYPH&lt;25&gt; spectral function develops an additional enhancement at low frequencies, which is no longer fully captured by the Breit-Wigner fits. Even though initially the GLYPH&lt;27&gt; spectral function can still be reasonably well described in terms of a single resonance, we find that below J = 0 05 : GLYPH&lt;2&gt; 2 GLYPH&lt;0&gt; 5 the description in terms of Breit-Wigner distribution becomes increasingly inaccurate also for the GLYPH&lt;27&gt; spectral function, as both spectral functions start to feature a strong enhancement at low frequency, which is no longer captured by a simple quasi-particle peak.

Eventually, the amount of explicit symmetry breaking is no longer large enough to guarantee the alignment of the order parameter in our finite volume system, such that for very small values of J the spectral functions of GLYPH&lt;25&gt; and GLYPH&lt;27&gt; e GLYPH&lt;11&gt; ectively become degenerate. Even though this is a finite volume artifact, it is also clear that extending the study to larger and larger lattices will only shift the problem towards smaller and smaller values of J , as in any finite system the alignment of the order parameter with the symmetry breaking axis will only be guaranteed above a certain amount of explicit symmetry breaking. We illustrate this problem in Fig. 12, where we compare the results for the spectral functions GLYPH&lt;26&gt; GLYPH&lt;25&gt; and GLYPH&lt;26&gt; GLYPH&lt;27&gt; obtained on L = 64 and L = 128 lattices.

One finds that for J = 0 05 : GLYPH&lt;2&gt; 2 GLYPH&lt;0&gt; 2 = 0 0125 the results ob-: tained on L = 64 and L = 128 lattices are in good agreement with each other, indicating the absence of finite-volume e GLYPH&lt;11&gt; ects. Decreasing the amount of explicit symmetry breaking J further, to approach the critical point, the distinction between GLYPH&lt;25&gt; and GLYPH&lt;27&gt; becomes less and less prominent, as the simulations develop a significant volume dependence in the vicinity of the critical point. In particular, for the L = 64 lattice, one clearly observes that at some point the relevant infrared cut-o GLYPH&lt;11&gt; is no longer set by J but rather by the finite system size, leading to J independent results for the spectral function for J glyph[lessorsimilar] 0 05 : GLYPH&lt;2&gt; 2 GLYPH&lt;0&gt; 8 = 0 000195. While the results obtained on : L = 128 lattices, continue to show an increase of the low frequency enhancement with decreasing J , a clear distinction between GLYPH&lt;25&gt; and GLYPH&lt;27&gt; modes ceases to exist for J glyph[lessorsimilar] 0 05 : GLYPH&lt;2&gt; 2 GLYPH&lt;0&gt; 9 = 0 000098 due to finite : volume e GLYPH&lt;11&gt; ects.

While previous studies of the critical dynamics of a Z 2 symmetric GLYPH&lt;30&gt; 4 theory in 2 + 1 dimensions dealt with this problem by performing simulations on extremely large lattices [2], it is worth pointing out that the problem is substantially more severe for the breaking of a continuous symmetry, where the orientation of order parameter field can rotate continuously over the course of the simulation, and we will therefore have to explore di GLYPH&lt;11&gt; erent strategies to study the dynamical critical behavior in the limit J ! 0 and T ! Tc .

## 5.3. A glance at critical dynamics

We now focus on the behavior of the spectral function in the vicinity of the critical point, realized by setting T GLYPH&lt;25&gt; Tc and J ! 0 in our simulations. Since in the vicinity of the critical point J ! 0 our finite volume simulations do not allow us to distinguish between GLYPH&lt;25&gt; and GLYPH&lt;27&gt; modes, we set J = 0 directly and investigate the behavior of the combined spectral function

$$\rho \equiv \frac { 1 } { N } \, \text{tr} \, \rho = \frac { 1 } { 4 } \rho _ { \sigma } + \frac { 3 } { 4 } \rho _ { \pi } \,.$$

In the limit T ! Tc , for ! ! 0 and p ! 0 the spectral function is expected to exhibit a scaling behavior of the form [2]

$$\rho ( s ^ { z } \omega, s p, s ^ { \frac { 1 } { \nu } } T _ { r } ) \, = \, s ^ { - ( 2 - \eta ) } \rho ( \omega, p, T _ { r } ) \,, \quad \ \ ( 3 4 )$$

where GLYPH&lt;17&gt; is the anomalous dimension, GLYPH&lt;23&gt; is the correlation-length exponent and z denotes the dynamic scaling exponent which can be used to classify the time-dependent critical behavior of

Figure 10: Spectral functions GLYPH&lt;26&gt; GLYPH&lt;25&gt; and GLYPH&lt;26&gt; GLYPH&lt;27&gt; for di GLYPH&lt;11&gt; erent amounts of explicit symmetry breaking J at a nearly critical temperature T = 17 4. Data is shown in log-lin : representation in the left panel and lin-lin representation in the right panel. Curves in the right panel show fits to a Breit-Wigner distribution (31).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-f6d77724-94d9-4444-9c55-f3cec42218e0.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-0ab71842-3981-4dba-b00f-6a66bef236b6.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-610d6c27-ce33-45d4-b2d9-36f34f014443.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-2a421b95-fe68-4e43-87d7-014f133b902f.jpg)

Figure 11: Dependence of Breit-Wigner fit parameters for GLYPH&lt;25&gt; and GLYPH&lt;27&gt; spectral functions on the explicit symmetry breaking J at a nearly critical temperature T = 17 4. :

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-c23b00ca-a3ca-489c-86a0-a17e263f3444.jpg)

a system. Taking GLYPH&lt;23&gt; = 0 7377(41) and : GLYPH&lt;13&gt; = 1 4531(104) from : [43] and using GLYPH&lt;13&gt;=GLYPH&lt;23&gt; = 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; we obtain GLYPH&lt;17&gt; = 0 03022. : Based on the analyis of Halperin and Hohenberg [26], the dynamic universality class is determined by the conserved quantities of the system along with their coupling to the order parameter.

We first note that for the Hamiltonian dynamics of the relativistic O (4) model considered in our study, we have a conserved energy

$$\epsilon \ = \ \int d ^ { 3 } x \left [ \frac { 1 } { 2 } \pi _ { a } ( x ) \pi _ { a } ( x ) + \frac { 1 } { 2 } ( \nabla \phi _ { a } ( x ) ) ( \nabla \phi _ { a } ( x ) ) \quad ( 3 5 ) \quad & \text{lead t} \\ + \frac { m ^ { 2 } } { 2 } \phi _ { a } ( x ) \phi _ { a } ( x ) + \frac { \lambda } { 4! N } ( \phi _ { a } ( x ) \phi _ { a } ( x ) ) ^ { 2 } \right ] \, \quad & \text{lead t} \\ \quad. \quad \. \quad \.$$

the conversed momentum

$$\Pi _ { i } & = \, \int d ^ { 3 } x \, \pi _ { a } ( x ) \nabla _ { i } \phi _ { a } ( x ) & ( 3 6 ) \quad \text{tem.th}$$

and a conserved O (4) current

$$j ^ { a b } & = \epsilon ^ { a b c d } \, \int d ^ { 3 } x \, \phi _ { c } ( x ) \pi _ { d } ( x ) & ( 3 7 ) \quad \text{tem si}$$

Since the Poisson brackets of the N -component order parameter GLYPH&lt;30&gt; a = 1 V R d x 3 GLYPH&lt;30&gt; a ( x ) with the conserved quantities GLYPH&lt;15&gt; and j ab are non-vanishing, the dynamics of these modes can a GLYPH&lt;11&gt; ect the critical dynamics of the order parameter, and the relativistic O (4) model with Hamiltonian dynamics does not belong to one of the standard dynamic universality classes according to the classification scheme of Halperin and Hohenberg [26]. However, it has been argued [46] that for negative values of the specific heat exponent GLYPH&lt;11&gt; &lt; 0 (which is the case for the O (4)) the coupling to the conserved energy is irrelevant. Since for J = 0 the order parameter becomes an N = 4 component field (i.e. carrying information about the orientation as well as the magnitude), the analysis of Wilczek and Rajagopal [47] suggests that the critical dynamics of the relativistic O (4) scalar theory follows an extension of model G , where the dynamical critical exponent z = d 2 can be determined from a renormalization group analysis.

Simulation results for 1 N tr GLYPH&lt;26&gt; are shown in Fig. 13 where we present results in the frequency and time domain obtained for various di GLYPH&lt;11&gt; erent lattice sizes between L = 48 and L = 256. We find that in the vicinity of the critical point, a fine time step of the numerical integrator is needed to correctly reproduce the late-time behavior of the spectral function; we have therefore decreased the time step in our numerical integration by a factor of four to GLYPH&lt;1&gt; t = 0 00125 and checked explicitly for our : L = 96 data that reducing the time step by an additional factor of four does not a GLYPH&lt;11&gt; ect the results.

Despite the strong finite-size dependence in our simulations, we also observe first indications of the emergence of an infrared power law dependence of the spectral function in the low frequency domain. Based on the scaling relation in Eq. (34) the critical spectral function is expected approach the following scaling behavior

$$\rho ( \omega, p = 0 ) \infty \, \omega ^ { \frac { 2 - \eta } { z } } \, \text{ \quad \ \ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{\ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text{ } \text$$

for an infinite system. Di GLYPH&lt;11&gt; erent curves in the left panel of Fig. 13 indicate power law fits, employing the values z = 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; and z = 3 2 for the dynamical critical exponent. = We find that our results for the spectral function at zero spatial momentum favor the value z = 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; of the so called 'conventional theory' of dynamic critical phenomena, which emerges when the critical divergencies of kinetic coe GLYPH&lt;14&gt; cients is not taken into account in the scaling analysis [26]. In fact this observation may be reinforced further by looking at the behavior of the spectral function in the time domain, where values of z &lt; 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; would lead to an increase of the spectral function function GLYPH&lt;26&gt; ( ) t / t 1 GLYPH&lt;0&gt; 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; z as a function at late times, which is clearly not observed in our simulations. Conversely, a value of z = 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; (or z &gt; 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; ) would lead to a logarithmic time dependence (or power law decay) of the spectral function GLYPH&lt;26&gt; ( ), which appears to be more consistent t with our results.

It is also evident from Fig. 13, that for any finite-size system the critical behavior of the spectral function GLYPH&lt;26&gt; ( t ! 1 ) is suppressed by the exponential decay with the auto-correlation time / exp( GLYPH&lt;0&gt; t =GLYPH&lt;24&gt; t ( L )). While for any finite system the autocorrelation time GLYPH&lt;24&gt; t ( L ) is finite, it diverges with increasing system size as GLYPH&lt;24&gt; t ( L ) / L z , corresponding to the well known phenomenon of critical slowing down, and the typical way to extract the dynamical critical exponent z in Monte-Carlo simulations. Our results for the analysis of the auto correlation time are compactly summarized in Fig. 14, where we we present fits to the late time exponential behavior of the spectral function GLYPH&lt;26&gt; ( ) t / exp( GLYPH&lt;0&gt; t =GLYPH&lt;24&gt; t ( L )) along with the results for GLYPH&lt;24&gt; t ( L ) shown in the inset. While for small volumes L &lt; 128 the scaling of the auto-correlation time GLYPH&lt;24&gt; t ( L ) appears to be consistent with the z = 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; predicted by the conventional theory, the behavior of GLYPH&lt;24&gt; t ( L ) for large volumes hints at a weaker divergence of the auto-correlation time consistent with z = 3 2 on the larger lat-= tices, as indicated by the solid and dashed curves in the inset of Fig. 14.

Since the frequency dependence of the critical spectral function GLYPH&lt;26&gt; ( !; p = 0) / ! 2 GLYPH&lt;0&gt; GLYPH&lt;17&gt; z , and the finite-size scaling of the autocorrelation time GLYPH&lt;24&gt; ( ) t / L z lead to di GLYPH&lt;11&gt; erent extractions of the dynamic critical exponent z , we are unable to determine the

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-06a09c34-4bc3-4b2a-844f-d5b4ba3d4f8c.jpg)

ω

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-98c15597-0fee-43ba-9939-bf06885c17dc.jpg)

ω

Figure 12: Spectral functions GLYPH&lt;26&gt; GLYPH&lt;25&gt; (left) and GLYPH&lt;26&gt; GLYPH&lt;27&gt; (right) at T GLYPH&lt;25&gt; Tc for di GLYPH&lt;11&gt; erent amounts of explicit symmetry breaking J and di GLYPH&lt;11&gt; erent lattice sizes L = 64 and L = 128. In the vicinity of the critical point ( J ! 0), fluctuations of the order parameter prohibit the distinction between GLYPH&lt;25&gt; and GLYPH&lt;27&gt; modes. Lines represent raw data, points shown to guide the eyes and illustrate errorbars.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-63b8e6e0-0721-42d1-a0b5-f423639551f3.jpg)

Figure 13: Behavior of the spectral function close to the critical point ( T GLYPH&lt;25&gt; Tc , J = 0) for di GLYPH&lt;11&gt; erent lattice sizes L in the frequency (left) and time domain (right). Lines represent raw data and points illustrate error estimates and guide the eyes in right figure.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-0fbc00dc-e49c-4b24-bba4-8d2f53c69551.jpg)

dynamical critical behavior precisely from our current simulations. One possible explanation of the observed discrepancies could be due to the fact that we have set the spatial momentum p = 0 prior to taking the limit ! ! 0 (or t ! 1 ), which may or may not a GLYPH&lt;11&gt; ect the critical scaling of the spectral function. In any case, it would be interesting to investigate the critical dynamics in more detail as a function of p and ! at non-vanishing Tr and J , to further elucidate on the structure of excitations in the vicinity of the critical point. However, this will require significant computational resources and is well beyond the scope of our present work.

and GLYPH&lt;26&gt; GLYPH&lt;27&gt; . While at very low temperatures, the classical-statistical approximation is inadequate and e GLYPH&lt;11&gt; ectively reduces to a meanfield approximation, we argued that a classical-statistical description becomes accurate in the vicinity of a second order phase transition and demonstrated some intriguing features of the spectral functions close to the crossover transition and in the vicinity of the O (4) critical point.

## 6. Conclusions

We have performed a detailed study of classical-statistical spectral functions in the relativistic O (4) model. While the static critical behavior is naturally reproduced correctly within the classical-statistical lattice approach, the focus of our study has been on the behavior of the real-time spectral functions GLYPH&lt;26&gt; GLYPH&lt;25&gt;

In the broader context of non-perturbative calculations of real-time spectral functions, the results from classical-statistical simulations reported in this paper may provide additional guidance to alternative theoretical approaches, based, e.g. on functional methods or analytic continuation of Euclidean correlation functions, where prior information on the structure of excitations is required to devise suitable ans¨ atze or e GLYPH&lt;14&gt; cient truncation schemes. Specifically, there is an interesting possibility to benchmark the quality of results obtained within functional approaches, based on a direct comparison of the results obtained in the classical-statistical limit. This is work in progress and will be reported elsewhere.

Figure 14: Extraction of the auto-correlation time GLYPH&lt;24&gt; t ( L ) from the spectral function close to the critical point ( T GLYPH&lt;25&gt; Tc , J = 0). Di GLYPH&lt;11&gt; erent curves show results for the spectral function GLYPH&lt;26&gt; ( t ; p = 0) for di GLYPH&lt;11&gt; erent lattice sizes L as a function of the scaled time t ( L L = 0) 3 2 = with L 0 = 128.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00912v1-figure-c10ec446-32d6-4268-8c91-127da39c1466.jpg)

## Acknowledgements

This research was supported by the Helmholtz International Center (HIC) for FAIR within the LOEWE initiative of the State of Hesse, and by the Deutsche Forschungsgemeinschaft (DFG) through the grant CRC-TR 211 'Strong-interaction matter under extreme conditions.'

## References

- [1] G. D. Moore, O. Saremi, Journal of High Energy Physics 2008 (09) (2008) 015. doi:10.1088/1126-6708/2008/09/015 .
- [2] J. Berges, S. Schlichting, D. Sexty, Nucl. Phys. B832 (2010) 228-240. arXiv:0912.3135 doi:10.1016/j.nuclphysb.2010.02.007 , .
- [3] M. Jarrell, J. Gubernatis, Physics Reports 269 (3) (1996) 133 - 195. doi: 10.1016/0370-1573(95)00074-7 .
- [4] M. Asakawa, Y. Nakahara, T. Hatsuda, Progress in Particle and Nuclear Physics 46 (2) (2001) 459 - 508. doi:10.1016/S0146-6410(01) 00150-8 .
- [5] Y. Burnier, A. Rothkopf, Phys. Rev. Lett. 111 (2013) 182003. doi:10. 1103/PhysRevLett.111.182003 .
- [6] B. B. Brandt, A. Francis, H. B. Meyer, D. Robaina, Phys. Rev. D 92 (2015) 094510. doi:10.1103/PhysRevD.92.094510 .
- [7] L. Schlessinger, Phys. Rev. 167 (1968) 1411-1423. doi:10.1103/ PhysRev.167.1411 .
- [8] R.-A. Tripolt, I. Haritan, J. Wambach, N. Moiseyev, Physics Letters B 774 (2017) 411 - 416. doi:10.1016/j.physletb.2017.10.001 .
- [9] D. Dudal, O. Oliveira, P. J. Silva, Phys. Rev. D 89 (2014) 014010. doi: 10.1103/PhysRevD.89.014010 .
- [10] R.-A. Tripolt, P. Gubler, M. Ulybyshev, L. Von Smekal, Comput. Phys. Commun. 237 (2019) 129-142. arXiv:1801.10348 doi:10.1016/j. , cpc.2018.11.012 .
- [11] J. Engels, O. Vogt, Nucl. Phys. B832 (2010) 538-566. arXiv:0911. 1939 doi:10.1016/j.nuclphysb.2010.02.006 , .
- [12] D. Roeder, J. Ruppert, D. H. Rischke, Nucl. Phys. A775 (2006) 127151. arXiv:hep-ph/0503042 , doi:10.1016/j.nuclphysa.2006. 05.007 .
- [13] K. Kamikado, N. Strodtho GLYPH&lt;11&gt; , L. von Smekal, J. Wambach, Eur. Phys. J. C74 (3) (2014) 2806. arXiv:1302.6199 , doi:10.1140/epjc/ s10052-014-2806-6 .
- [14] R.-A. Tripolt, N. Strodtho GLYPH&lt;11&gt; , L. von Smekal, J. Wambach, Phys. Rev. D89 (3) (2014) 034010. arXiv:1311.0630 doi:10.1103/PhysRevD. , 89.034010 .
- [15] R.-A. Tripolt, L. von Smekal, J. Wambach, Phys. Rev. D90 (7) (2014) 074031. arXiv:1408.3512 doi:10.1103/PhysRevD.90.074031 , .
- [16] D. Mesterhazy, J. H. Stockemer, Y. Tanizaki, Phys. Rev. D92 (7) (2015) 076001. arXiv:1504.07268 doi:10.1103/PhysRevD.92.076001 , .
- [17] J. M. Pawlowski, N. Strodtho GLYPH&lt;11&gt; , Phys. Rev. D92 (9) (2015) 094009. arXiv:1508.01160 doi:10.1103/PhysRevD.92.094009 , .
- [18] N. Strodtho GLYPH&lt;11&gt; , Phys. Rev. D95 (7) (2017) 076002. arXiv:1611.05036 , doi:10.1103/PhysRevD.95.076002 .
- [19] J. M. Pawlowski, N. Strodtho GLYPH&lt;11&gt; , N. Wink, Phys. Rev. D98 (7) (2018) 074008. arXiv:1711.07444 doi:10.1103/PhysRevD.98.074008 , .
- [20] J. A. Mueller, C. S. Fischer, D. Nickel, Eur. Phys. J. C70 (2010) 1037-1049. arXiv:1009.3762 , doi:10.1140/epjc/ s10052-010-1499-8 .
- [21] C. S. Fischer, J. M. Pawlowski, A. Rothkopf, C. A. Welzbacher, Phys. Rev. D98 (1) (2018) 014009. arXiv:1705.03207 , doi:10.1103/ PhysRevD.98.014009 .
- [22] G. Aarts, J. Smit, Nucl. Phys. B511 (1998) 451-478. arXiv:hep-ph/ 9707342 doi:10.1016/S0550-3213(97)00723-2 , .
- [23] R. Kubo, Journal of the Physical Society of Japan 12 (6) (1957) 570-586. doi:10.1143/JPSJ.12.570 .
- [24] P. C. Martin, J. Schwinger, Phys. Rev. 115 (1959) 1342-1373. doi: 10.1103/PhysRev.115.1342 .
- [25] G. Aarts, Phys. Lett. B518 (2001) 315-322. arXiv:hep-ph/0108125 , doi:10.1016/S0370-2693(01)01081-4 .
- [26] P. C. Hohenberg, B. I. Halperin, Rev. Mod. Phys. 49 (1977) 435-479. doi:10.1103/RevModPhys.49.435 .
- [27] R. D. Pisarski, F. Wilczek, Phys. Rev. D 29 (1984) 338-341. doi:10. 1103/PhysRevD.29.338 .
- [28] F. Wilczek, International Journal of Modern Physics A 07 (16) (1992) 3911-3925. doi:10.1142/S0217751X92001757 .
- [29] K. Rajagopal, F. Wilczek, Nuclear Physics B 399 (2) (1993) 395 - 425. doi:10.1016/0550-3213(93)90502-G .
- [30] A. Polyakov, Physics Letters B 59 (1) (1975) 79 - 81. doi:10.1016/ 0370-2693(75)90161-6 .
- [31] V. Novikov, M. Shifman, A. Vainshtein, V. Zakharov, Physics Reports 116 (3) (1984) 103 - 171. doi:10.1016/0370-1573(84)90021-8 .
- [32] J. O. Andersen, D. Boer, H. J. Warringa, Phys. Rev. D69 (2004) 076006. arXiv:hep-ph/0309091 doi:10.1103/PhysRevD.69.076006 , .
- [33] E. Seel, D. Smith, S. Lottini, F. Giacosa, JHEP 07 (2013) 010. arXiv: 1209.4243 doi:10.1007/JHEP07(2013)010 , .
- [34] G. Parisi, Statistical Field Theory, Addison-Wesley (Frontiers in Physics 66), 1988.
- [35] J. Berges, AIP Conf. Proc. 739 (2005) 3-62. arXiv:hep-ph/0409233 , doi:10.1063/1.1843591 .
- [36] A. Pieiro Orioli, J. Berges, Phys. Rev. Lett. 122 (15) (2019) 150401. arXiv:1810.12392 doi:10.1103/PhysRevLett.122.150401 , .
- [37] M. Hasenbusch, T. Torok, J. Phys. A32 (1999) 6361. arXiv:cond-mat/ 9904408 doi:10.1088/0305-4470/32/36/301 , .
- [38] K. Kanaya, S. Kaya, Phys. Rev. D 51 (1995) 2404-2410. doi:10.1103/ PhysRevD.51.2404 .
- [39] J. Engels, T. Mendes, Nucl. Phys. B572 (2000) 289-304. arXiv: hep-lat/9911028 doi:10.1016/S0550-3213(00)00046-8 , .
- [40] A. Pelissetto, E. Vicari, Phys. Rept. 368 (2002) 549-727. arXiv: cond-mat/0012164 doi:10.1016/S0370-1573(02)00219-3 , .
- [41] J. Engels, S. Holtmann, T. Mendes, T. Schulze, Physics Letters B 514 (3) (2001) 299 - 308. doi:10.1016/S0370-2693(01)00798-5 .
- [42] J. Engels, L. Fromme, M. Seniuch, Nucl. Phys. B675 (2003) 533554. arXiv:hep-lat/0307032 doi:10.1016/j.nuclphysb.2003. , 09.060 .
- [43] J. Engels, F. Karsch, Phys. Rev. D90 (1) (2014) 014501. arXiv:1402. 5302 doi:10.1103/PhysRevD.90.014501 , .
- [44] A. M. Ferrenberg, R. H. Swendsen, Phys. Rev. Lett. 61 (1988) 26352638. doi:10.1103/PhysRevLett.61.2635 .
- [45] P. Springer, Finite Volume E GLYPH&lt;11&gt; ects and the Binder Cumulant in O N ( )-Model, M.Sc. Thesis, Technische Universi¨ at M¨nchen, 2013. u
- [46] B. I. Halperin, P. C. Hohenberg, S.-k. Ma, Phys. Rev. B10 (1974) 139153. doi:10.1103/PhysRevB.10.139 .
- [47] K. Rajagopal, F. Wilczek, Nucl. Phys. B399 (1993) 395-425. arXiv: hep-ph/9210253 doi:10.1016/0550-3213(93)90502-G , .
Sören Schlichting Dominik Smith Lorenz von Smekal
Spectral functions and critical dynamics of the $O(4)$ model from classical-statistical lattice simulations
We calculate spectral functions of the relativistic $O(4)$ model from
real-time lattice simulations in classical-statistical field theory. While in
the low and high temperature phase of the model, the spectral functions of
longitudinal $(\sigma)$ and transverse $(\pi)$ modes are well described by
relativistic quasi-particle peaks, we find a highly non-trivial behavior of the
spectral functions in the cross over region, where additional structures
appear. Similarly, we observe a significant broadening of the quasi-particle
peaks, when the amount explicit $O(4)$ symmetry breaking is reduced. We further
demonstrate that in the vicinity of the $O(4)$ critical point, the spectral
functions develop an infrared power law associated with the critical dynamics,
and comment on the extraction of the dynamical critical exponent $z$ from our
simulations.