text:
# Quantum Monte Carlo simulation of thin magnetic films

## I Introduction

The driving force behind much of the research on thin magnetic films is their application in data storage devices. Magnetic thin films also display many remarkable physical phenomena, such as the reorientation effect,Bland and Heinrich (1994) in which the axis of magnetization changes as a function of film thickness, temperature and applied fields. Many theoretical methods, such as ab-initio calculations,Hjortstam et al. (1997) mean-field theories,Moschel and Usadel (1994); Jensen and Bennemann (1998) classical Monte Carlo simulations,Erickson and Mills (1991a); Iglesias et al. (1999) Green’s function methodsJensen et al. (1999); Ecker et al. (1999); Fröbrich et al. (2000a, b); Fröbrich et al. (2002); Erickson and Mills (1991b) and Schwinger bosonsTimm and Jensen (2000) have been applied to thin-film systems. Much of this work has focused on using the Heisenberg model to study the reorientation effect. The ground state and the lowest (one-magnon) excitations are known for the two-dimensional ferromagnetic Heisenberg model, but there is no closed form analytic solution at finite temperatures. Since results obtained using different approximate methods differ a great deal from each other, there is a need for an unbiased check of the various methods used; this is achieved by quantum Monte Carlo (QMC) calculations. In Ref. Ecker et al., 1999, for example, it has been shown by comparing with QMC resultsTimm et al. (1998); Henelius et al. (2000) that the Tyablicov Tyablicov (1959) decoupling (random phase approximation: RPA) is a very good approximation for the magnetization of a spin $`S=1/2`$ monolayer in an external magnetic field (perpendicular to the film plane). The main purpose of the present paper is to show the feasibility of large-scale quantum Monte Carlo (QMC) calculations, free of systematic errors, in the study of thin magnetic films. In order to achieve this goal we have included higher spins, local anisotropies, and a transverse magnetic field in the operator-loop formulation of the stochastic series expansion (SSE) QMC method.Sandvik (1999)

Two often-cited trends in quantum Monte Carlo development are the emergence of methods free of systematic errors,Sandvik and Kurkijärvi (1991); Prokof’ev et al. (1996); Beard and Wiese (1996); Rombouts et al. (1999); Henelius et al. (2000) and the development of highly efficient loop-cluster algorithms.Evertz et al. (1993); Kawashima et al. (1994); Prokof’ev et al. (1996); Beard and Wiese (1996); Sandvik (1999); Syljuåsen and Sandvik  In this work we particularly want to emphasize the general applicability of the SSE operator-loop method, which makes it possible to use the same algorithm to study a wide variety of Hamiltonians. Whereas previously it was necessary to rewrite large sections of the computer code when changing the model, one can now use the same code (compiled only once) to simulate a wide range of different systems. The user of the program no longer necessarily needs detailed knowledge of the algorithm and code to be able to conduct a thorough study of many quantum spin models in any dimension with, for example, non-zero magnetic fields, anisotropies and varying spin size.

A brief description of the QMC method is given in Sec. II. Thereafter we discuss the applicability of the method and the introduction of general spin size and a transverse field in the SSE operator-loop algorithm. In Sec. III we compare some examples of the QMC simulations to approximate theoretical approaches. Finally, we comment on possible future applications of QMC in the context of thin magnetic films.

## II SSE cluster-loop algorithm

There are excellent descriptions of the SSE loop algorithm,Sandvik (1999); Syljuåsen and Sandvik  so we only give a brief summary here in order to introduce the general framework of the method. We will, however, try to describe the main features of the method pictorially. The focus is on the new aspects that arise when introducing arbitrary spin size and a transverse field.

We consider a lattice spin model described by a Hamiltonian $`H`$. The SSE method relies on a Taylor expansion of the partition function $`Z`$:

$$Z=\underset{\alpha }{}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(\beta )^n}{n!}\alpha |H^n|\alpha ,$$
(1)
where $`|\alpha `$ are basis states in which the matrix element above can be evaluated, and $`\beta `$ is the inverse temperature.

To describe the updating procedure, we write the Hamiltonian as a sum over all $`M`$ bonds representing interacting spins in the system

$$H=\underset{b=1}{\overset{M}{}}H_b.$$
(2)
The bond operator $`H_b`$ can be decomposed into its diagonal and off-diagonal parts,

$$H_b=H_{D,b}+H_{O,b},$$
(3)
where subscript $`D`$ denotes a diagonal operator and $`O`$ an off-diagonal operator. For a ferromagnetic Heisenberg model these two operators are of the form

$$H_{D,b}=S_{i(b)}^zS_{j(b)}^z$$
(4)
and

$$H_{O,b}=S_{i(b)}^+S_{j(b)}^{}+S_{i(b)}^{}S_{j(b)}^+,$$
(5)
where $`i(b)`$ and $`j(b)`$ denote the two spins connected by bond $`b`$. If we introduce a cutoff order $`L`$ in the Taylor expansion (which, when done properly, does not cause any systematic errorsSyljuåsen and Sandvik ), and include additional unit operators $`I`$, the expansion can be rewritten in the form

$$H=\underset{\alpha }{}\underset{S_L}{}\frac{\beta ^n(Ln)!}{L!}\alpha |S_L|\alpha ,$$
(6)
where $`S_L`$ is an operator string

$$S_L=\underset{p=1}{\overset{L}{}}H_p,$$
(7)
with $`H_p\{H_{D,b},H_{O,b},I\}`$. The Monte Carlo procedure must sample the space of all states $`|\alpha `$ and all operator sequences $`S_L`$.

We next consider the SSE space in more detail. Denoting a propagated state by

$$|\alpha (p)=\underset{i=1}{\overset{p}{}}H_i|\alpha ,$$
(8)
the matrix element in Eq. (6) can be written as a product of elements of the form $`\alpha (p)|H_b|\alpha (p1)`$, where the bond-operator $`H_b`$ only acts on two spins. For the spin-1/2 Heisenberg model the only elements that can appear are $`,|H_b|,`$, $`,|H_b|,`$, $`,|H_b|,`$, and the spin-reversed versions of the same set. From now on we refer to these matrix elements, consisting of four spin states and a bond operator, as “vertices”. The matrix element in Eq. (6) can thus be viewed as a list of vertices of the above kind. In Fig. 1 we depict such a list graphically.

In the operator-loop algorithm, two basic updates ensure that the complete SSE space is sampled. The diagonal update attempts to exchange diagonal operators $`H_{D,b}`$ and unit operators $`I`$. The probability for inserting a diagonal operator (exchanging it for a unit operator) at position $`p`$ in the operator sequence is

$$P_{\text{insert}}=\frac{M\beta \alpha (p)|H_{D,b}|\alpha (p)}{Ln},$$
(9)
while the probability for removing a diagonal operator is

$$P_{\text{remove}}=\frac{Ln+1}{M\beta \alpha (p)|H_{D,b}|\alpha (p)}.$$
(10)
In a diagonal update, one exchange attempt is made for each diagonal and unit operator. A typical outcome of a diagonal update is shown in Fig. 2.

The second type of update is a global operator-loop update, which leaves unit operators unaffected. The idea of the loop move is to form and flip a closed loop of spins in the vertex list. In the process both the affected vertices and states are changed. The operator-loop update together with the above diagonal update ensure that the complete SSE configuration space is sampled.

The loop move is easy to formulate. In Fig. 3 we illustrate a loop move graphically. We have removed all the unit operators $`I`$ from the operator sequence in Fig. 2. Furthermore we only show the spins that are members of a vertex, and we turn the configuration into a linked list by connecting the same spin at different vertices by vertical lines, which we refer to as “links”. Next we describe the loop move in more detail.

First a random spin, belonging to a vertex, is selected and flipped, see Fig. 3(a). We refer to this spin as the “entrance” spin to the first vertex. For the loop to proceed, we select an “exit” spin, also belonging to the same vertex, see Fig. 3(b). The probability of choosing a given exit spin is proportional to the vertex matrix element, that results from flipping the entrance and exit legs. This choice of exit spin ensures detailed balance.Syljuåsen and Sandvik The exit spin is linked to another spin in the linked vertex list, and this spin is chosen as the entrance spin to the next vertex, see Fig. 3(c). In this manner a loop of spins on the space-time lattice is formed and flipped. Note that when the first spin is flipped a discontinuity is introduced into the vertex list. A link-discontinuity appears when two spins linked together are not in the same state. When, for example, the initial entrance and exit spins (assuming they are different) are flipped, two discontinuities are introduced in the linked list, see Fig. 3(b), where a link-discontinuity is marked by a short horizontal bar. One of the discontinuities gets propagated by the loop until the loop passes through the initial spin a second time, when the discontinuities “annihilate” each other and the loop closes, see Fig. 3(d). A practical and sufficient criterion for closing the loop is that the exit spin (after having been flipped) is in the same state as the next entrance spin (before flipping it). When this is the case, all discontinuities have been removed and the loop can be closed.

### II.1 Applicability of the algorithm

The SSE loop algorithm is extremely general and can be efficiently applied to any model with a positive definite partition function (if this is not the case, one encounters the usual difficulties due to the sign problemHenelius and Sandvik (2000)). This means that the algorithm can be applied to ferromagnetic and antiferromagnetic (on bipartite lattices) Heisenberg models of any quantum spin size in any dimension, including magnetic fields and local as well as exchange anisotropies. The interactions need not be short-ranged, as long as they do not cause frustration (which leads to the sign problem).

From a programming point of view, one need only calculate the energies of all possible vertices (which typically is a simple task) and generate a table of exit probabilities as a function of four variables: the initial vertex type, the entrance spin, the entrance spin state and the exit spin.

As an illustration we have written a code that works for a general $`d`$-dimensional Heisenberg model of the form

$`H=`$ $`{\displaystyle \underset{ij}{}}\left[J_{ij}^zS_i^zS_j^z+J_{ij}^t(S_i^+S_j^{}+S_i^{}S_j^+)\right]`$ (11)
$`{\displaystyle \underset{i}{}}\left[K_2(S_i^z)^2+K_4(S_i^z)^4\right]`$
$`{\displaystyle \underset{i}{}}\left[B^zS_i^z+B^xS_i^x\right],`$
where the exchange interaction is given by $`J_{ij}`$ (negative for antiferromagnets and positive for ferromagnets). $`𝐒`$ denotes the quantum spin and $`K_2`$ and $`K_4`$ are the second and fourth-order local uniaxial anisotropies. $`𝐁`$ is an external magnetic field. In this study we limit ourselves to a nearest neighbour interaction, but in principle one can include long-range interactions. In the following we describe the modifications to the loop algorithm that allow the introduction of general spin size and a transverse field.

### II.2 Inclusion of arbitrary spin

For the sake of simplicity, the SSE loop-operator algorithm was originally describedSandvik (1999) for a spin-1/2 system, but the generalization to a spin-$`S`$ system is straightforward. For an arbitrary spin system there are more allowed vertices, but the probability for inserting and extracting diagonal operators is still given by Eq. (9) and Eq. (10) respectively. The only necessary change in the loop-update is the choice of the initial spin state. For the spin-1/2 model the initial entrance spin was simply flipped, but for the spin-$`S`$ model one can randomly choose among the $`2S`$ spin states that differ from the initial spin state. All other aspects of the algorithm remain the same.

### II.3 Inclusion of a transverse field

Previously a transverse field has been included in the SSE algorithm using local updates.Henelius et al. (2000) In the present work we include the transverse field directly in the global loop update. There are very likely many different ways to formulate a loop-update which includes a transverse magnetic field. Here we will not make an exhaustive study, but rather present one possible algorithm, which we found particularly simple to implement. It is probably not the most efficient update, and we are currently making a more detailed study of more efficient algorithms for including the transverse field.

Including a transverse magnetic field is more involved because the total $`z`$ component of the magnetization is no longer a conserved quantity, $`[H,_iS_i^z]0`$, due to the presence of single lowering and raising operators ($`2S^x=S^++S^{}`$) in the Hamiltonian. As a consequence additional “flip vertices”, with different magnetization in the initial and final states (such as $`,|H_b|,`$), are also allowed in the vertex-list. The presence of flip vertices implies that, in general, one can exit the same exit spin with several different spin states. The increase in the number of exit states means that the exit probability is a function of five variables instead of four: original vertex, entrance spin, entrance state, exit spin and exit state. For a model that conserves the total $`z`$ component of the magnetization the exit state is a function of the other four variables.

To include a magnetic field, we also have to reconsider how to start, and end, the operator-loop. With no magnetic field we introduce an initial link-discontinuity by flipping the first entrance spin. A second link-discontinuity appears when the first exit spin is flipped (unless the entrance and exit spins coincide, in which case the loop closes immediately). The second discontinuity is propagated as the loop progresses until the two discontinuities annihilate each other and the loop closes. With flip vertices allowed, the number of link-discontinuities can change by the introduction or removal of single flip vertices. If one therefore chooses to start, and end, the loop in the same way as previously described, the loop can close without having corrected for the initial discontinuity. In such a case it would be necessary to return to the starting point of the loop and continue execution until the last discontinuity is removed.

Here we chose a slightly different approach, which can be implemented by a very minor change in the original formulation. Instead of picking an initial entrance spin and changing its state we leave its state unchanged. Without a magnetic field the loop would close immediately. With a magnetic field this possibility still exists, but the loop may also start with a flip vertex and thereafter proceed in a “normal” fashion until another flip vertex is inserted, or removed, when the loop is closed according to the usual criterion (the exit spin and the following entrance spin are in the same state). We have compared high-precision QMC data for small systems with exact diagonalization, to ensure that detailed balance is indeed satisfied, but we do not provide a strict proof here.

Fig. 4 depicts an example of a loop move for the case of a non-zero transverse field. In Fig. 4(a) an initial entrance spin is selected and left in its original state. An exit spin is selected and flipped (Fig. 4(b)), resulting in a flip vertex. The next entrance spin can be seen in Fig. 4(c), with a corresponding exit spin in Fig. 4(d). Note that the exit state in Fig. 4(d) is unchanged, resulting in a second flip vertex, leading to termination of the loop. The loop does not close on itself in this case, but starts and ends at two different spins.

The transverse magnetization is particularly easy to calculate within the SSE formulation, since it simply is equal to the average number of flip vertices, $`N_{\text{flip}}`$, in the operator sequence,Henelius et al. (2000)

$$M_x=\frac{1}{\beta }N_{\text{flip}}.$$
(12)
To conclude this section we note that only two changes need to be made to add a transverse field. First of all the initial entrance spin state should be left unchanged. Secondly, the exit probability is now also a function of exit state, since there are, in general, several possible exit states for a given exit spin. Next we apply this algorithm to a two-dimensional ferromagnetic system.

## III Application to thin magnetic films

The simplest effective model for a thin magnetic film is a ferromagnetic Heisenberg monolayer. However, the Mermin-Wagner theoremMermin and Wagner (1966) tells us that this model cannot have a finite critical temperature in two dimensions. The continuous symmetry can be explicitly broken by an anisotropy in order to induce a finite critical temperature. Typical experimentally observed values of the local anisotropy in 3$`d`$ transition metal films are about two orders of magnitude smaller than the exchange coupling, and one might therefore expect a very small critical temperature. However, it turns out that the critical temperature contains a logarithmic dependence on the anisotropy,Bander and Mills (1988); Erickson and Mills (1991a); Timm and Jensen (2000) so that very small anisotropies induce a critical temperature of the order of the coupling constant. In this section we concentrate on results for the magnetization as a function of temperature. The model is given by Eq. (11), where the double sum is over all nearest neighbours on a two-dimensional square lattice. Since the main purpose of the present paper is to show the feasibility of the quantum Monte Carlo method, and not to explore the whole parameter range of the model as defined in Eq. (11), we limit ourselves to some illustrative examples. From now on we will assume an isotropic exchange interaction ($`J_{ij}^t=J_{ij}^z/2`$), zero vertical magnetic field ($`B^z=0`$) and only second-order anisotropy ($`K_4=0`$). First we discuss how the finite-size effects were treated. Thereafter we compare the QMC results with approximate theoretical approaches. Finally, we also show a case where a transverse field drives the magnetization in the plane.

### III.1 Finite-size effects

The first problem one encounters in a finite-size system with spin-inversion symmetry is that the magnetization should vanish because opposite spin orientations occur equally likely. One can circumvent this problem by calculating either the absolute value of the magnetization, or the magnetization squared (and take the root afterwards). Both approaches are equivalent in the thermodynamic limit, and, since the absolute value of the magnetization is found to converge faster, we show only the former. The magnetization curves for a particular value of the anisotropy, $`K_2/J=0.01`$, and different system sizes are shown in Fig. 5. The finite-size effects increase closer to the Curie temperature $`T_c`$, as expected for a second-order phase transition. For low temperatures, however, one can see that it is possible to extract data which have converged (within statistical error) with respect to system size. In the next section we restrict ourselves to results for the magnetization that show no discernible finite-size effects.

The critical temperature can be determined directly. We have used the Binder ratioBinder (1981) to extract the value of $`T_c`$. By plotting the ratio of two moments of the magnetization,

$$Q=\frac{M_z^4^{\frac{1}{4}}}{M_z^2^{\frac{1}{2}}},$$
(13)
the finite-size effects around $`T_c`$ should largely cancel, and curves for different system sizes are expected to intersect at $`T_c`$. In Fig. 6 we show an example for $`K_2/J=1.5`$, where one can clearly see how data for different system sizes intersect at one point ($`T_c/J`$).

In transition-metal thin films, the second-order anisotropy is believed to be of the order $`K_2/J=0.01`$. The relatively small energy scale of the anisotropy induces large finite-size effects. In order to see the effects of a small energy scale one needs to examine large systems at low temperatures. In Fig. 7 we show the Binder ratio as a function of temperature and system size for $`K_2/J=0.01`$. It seems that results for linear system sizes $`N=4,8`$ and 16 have converged around $`T_c/J1.05`$, but as the system size is increased strong corrections appear and push $`T_c/J`$ down to around $`T_c/J0.9`$. This is a clear case where it would be dangerous to draw conclusions from a study of small system sizes. The results for system size $`N=256`$ have statistical errors which are slightly larger than the symbol size, while all other statistical errors are much smaller. To determine $`T_c`$ very accurately for such a small anisotropy one would need accurate data for even larger system sizes. However, the current precision is enough for the comparisons with results of approximate theoretical methods that we present in the next section.

### III.2 Magnetization

We compare the QMC data with Schwinger Boson mean-fieldTimm and Jensen (2000) and many-body Green’s functions calculations.Fröbrich et al. (2000b); Fröbrich et al. (2002) Both methods include spin wave excitations approximately and represent significant improvements over simple mean field theories in which magnon excitations are neglected completely. Of these two methods the Schwinger boson theories are numerically less demanding and much easier to extend to arbitrary spin than the RPA approach.Timm and Jensen (2000)

In the Schwinger boson theory the Heisenberg model is mapped onto an equivalent bosonic system. This can be done by using the SU(2) symmetry in spin space of the Heisenberg model. The SU(2) model can then be generalized to an SU($`N`$) model, containing $`N`$ bosons per site. In the limit $`N\mathrm{}`$ mean-field theory becomes exact, and in this section we label mean-field SU($`N\mathrm{}`$) results by “SU($`N`$)”. Using the local equivalence between the SU(2) and O(3) groups, this can be repeated for an O(3) model, and we label mean-field O($`N\mathrm{}`$) results by “O($`N`$)”.

The many-body Green’s function calculations for the magnetization are done by a procedure where one works at the second level of the hierarchy of equations of motion for the Green’s functions.Fröbrich et al. (2002) This allows an exact treatment of the terms stemming from the single-ion anisotropy, whereas the exchange interaction terms are treated by the TyablicovTyablicov (1959), or random-phase approximation (RPA), and from now on Green’s function results will be labeled by RPA. This procedure is an improvement over the Anderson-Callen decouplingFröbrich et al. (2000b), in which the single-ion anisotropy terms are decoupled at the level of the lowest-order equation of motion, which is a good approximation only for small anisotropies.Fröbrich et al. (2002) The results in this subsection are calculated using the exact treatment of the single-ion anisotropy. Owing to problems with numerical stability, the calculations for the reorientation of the magnetization discussed in the next subsection are done with the Anderson-Callen decoupling.

We compare magnetization curves for three values of the anisotropy covering three orders of magnitude, $`K_2/J=0.01`$, 0.2 and 1.5. The smallest value is of greatest relevance for experiments, but it is of interest to see how well the analytic methods work outside this region as well. In Fig. 8 we see that, for $`K_2/J=0.01`$ the RPA calculations give a rather accurate magnetization curve, while the O($`N`$) and SU($`N`$) theories give a low and high estimate, respectively. In the low-temperature limit both the SU($`N`$) theory and the RPA calculation recover the correct spin-wave result. Note also by how much a simple mean-field (MF) theory overestimates the magnetization ($`T_c^{\mathrm{MF}}2.7T_c^{\mathrm{QMC}}!`$). As is well known, mean-field theory also totally fails at low temperatures (exponential instead of power law behaviour), due to the neglect of spin waves. In Fig. 9 we see that for $`K_2/J=0.2`$ the RPA still gives very accurate values at low temperatures, while overestimating the magnetization at higher temperatures. The O($`N`$) theory happens to give a good estimate of $`T_c`$, while the SU($`N`$) again overestimates the magnetization. In Fig. 9 the QMC points are connected by straight lines, except for the line between $`T_c`$ and the highest temperature below $`T_c`$. This line is a fit to the Ising-like critical behavior

$$M_z(TT_c)^{\frac{1}{8}}.$$
(14)
The good fit indicates that the QMC data have come close enough to $`T_c`$ for the critical behavior to set in. For the largest anisotropy, $`K_2/J=1.5`$ (Fig. 10), the RPA calculation again gives the most accurate estimate, but all approximate curves result in too high a magnetization. For high temperatures the SU($`N`$) theory yields a larger magnetization than simple mean-field theory. In Fig. 10 we have also included the quantum Monte Carlo data raised to the eighth power. The resulting straight line is further evidence of how close to the critical point the simulation has come.

The large deviations of the approximate theoretical treatments from QMC, in particular close to $`T_c`$ and for large anisotropies, indicate that neither the Green’s function approach nor the Schwinger boson method treats spin wave interactions in a satisfactory way. It is interesting to note that that RPA and SU($`N`$) give a (mean-field) exponent of $`1/2`$, whereas O($`N`$) gives $`1/3`$.

In this subsection we have compared QMC data with approximate theoretical results for the magnetization curve. Except for the value of $`T_c`$, the QMC data shown above have converged in system size within error bars that are not discernible in the figures. The largest system size used is $`256\times 256`$ spins. Error bars are shown for the Binder estimate of $`T_c`$. As can be seen the error increases with decreasing anisotropy, which is to be expected.

### III.3 Reorientation

We consider two examples for the reorientation of the magnetization in a transverse field. For an anisotropic spin-2 model we have calculated the vertical and transverse components of the magnetization, $`M_z`$ and $`M_x`$, as a function of a transverse magnetic field at a fixed temperature.

The anisotropy favors an out-of-plane magnetization, while the transverse field wants the magnetization to be in the plane. This competition results in quite interesting phase diagrams, where the order of the transition in general depends on the order of the anisotropy.Timm and Jensen (2000) Here we concentrate on two examples, to demonstrate the applicability of the QMC approach.

In Fig. 11 we show the magnetization curves for a small anisotropy, $`K_2/J=0.01`$, at a temperature $`T/J=1`$. Primarily we show QMC results for system size $`128\times 128`$, but in addition we show results for $`64\times 64`$ in the region where they differ from each other. For small and intermediate fields the results have converged, while we again see increasing finite-size effects closer to the critical field. The trend is, however, clear. The transverse magnetization does appear to increase linearly with the transverse field, up to the highest fields for which the results have converged.

The SU($`N`$) and O($`N`$) calculations both show a linear dependence on the transverse magnetization, $`M_x`$, similar to QMC, whereas the vertical component is over- or underestimated, respectively. The SU($`N`$) results for $`M_z`$ are good at weak field, though. The mean-field results for $`M_x`$ are also linear but mean-field theory consistently overestimates $`M_z`$. The RPA calculation with the approximate decoupling of the anisotropy termsFröbrich et al. (2000b) follows the QMC curves at small values of the transverse field, but the reorientation occurs at a considerably smaller critical field than in QMC and in the other approximations. Also, in the RPA the magnetization is not a unique function of field close to the transition. The same system has been studied for a larger anisotropy, $`K_2/J=0.20`$, in Fig. 12. Here the SU($`N`$) and O($`N`$) solutions are mean-field-like and agree very well with the QMC solution, except close to the transition, whereas the RPA behaves as in Fig. 11. Close to the transition the SU($`N`$), O($`N`$), and mean-field solutions overestimate the vertical magnetization. The QMC results display a small deviation from the linear increase in the in-plane magnetization, which is not reproduced by the Schwinger boson or mean-field methods.

## IV Conclusion

This work shows the feasibility of using large-scale QMC calculations to examine microscopic thin film models. The QMC approach can be used both as a sanity check on approximate theoretical treatments and as a method in its own right. Our results indicate that in the absence of a transverse field, RPA with an exact treatment of the anisotropy termsFröbrich et al. (2002) appears to be more accurate than the Schwinger boson calculation. Both methods, however, represent an improvement over simple mean-field theory for small anisotropies but QMC reveals weaknesses in the approximate theories at large anisotropies and close to the Curie temperature. In the presence of a transverse field the results obtained from Schwinger boson methods are quite mean-field-like, which turns out to be appropriate for most of the field range. The reason that the RPA is worse in this case is probably due to the Anderson-Callen decoupling of the anisotropy terms in RPA.Fröbrich et al. (2000b) This was necessary because the more accurate treatmentFröbrich et al. (2002) led to numerical difficulties when applied to the reorientation problem. To make more definite conclusions about the merits of different approximate methods, we would have to investigate a much larger parameter space, including higher anisotropies, general spin, temperature, and transverse field and the extension to several layers. We leave such an investigation as work for the future.

It is possible to extend the QMC calculations to include several layers with arbitrary inter-layer coupling, as well as, for example, anisotropies in the exchange coupling. Unfortunately, the dipole interaction introduces frustration and thereby the sign problem, and is therefore currently out of reach for QMC studies. On the other hand, long-range ferromagnetic couplings are not a problem. The implementation of the dipole coupling also leads to problems in the Schwinger boson theory but is possible in a Green’s function description.Jensen et al. (1999); Fröbrich et al. (2000b) A very exciting recent development in QMC is the introduction of directed loop moves,Syljuåsen and Sandvik  which, according to our initial calculations, can reduce the autocorrelation time by one order of magnitude in spin-1 systems, and therefore make it possible to reach substantially larger system sizes. However, it is not clear that the directed loops are as easy to implement for a general model as the method used in this work.

We are grateful to A. Sandvik for fruitful discussions. P.H. acknowledges support by the Swedish Research Council and Ella och Georg Ehrnrooths stiftelse.


format: text

text:
# Real-time pion propagation in finite-temperature QCD

## I Introduction

Properties of hadrons at high temperatures and densities are of great interest from both experimental and theoretical perspectives. One motivation for studying temperature effects on hadrons comes from the suggestion that some features of the dilepton spectrum observed in heavy-ion collisions can be explained by the modification of masses and widths of mesons by the thermal medium dilepton . Nevertheless, reliable information on the temperature modification of hadronic properties is still lacking. Lattice simulations, which rely on the imaginary-time formulation of quantum field theory, have serious difficulties with real-time quantities.111Some progress, however, may have been made recently lattice-dilepton . The absence of Lorentz invariance at finite temperature implies that there is no direct relationship between real-time characteristics of hadrons (for example, the so-called “pole masses,” which are supposedly the positions of poles in propagators) and quantities that can be extracted from Euclidean propagators (e.g., the “screening masses,” which characterize the exponential falloff of static Euclidean correlators). Thus, as a rule, lattice measurements of correlation functions at finite temperature cannot be used to draw conclusions about real-time propagation of hadrons.

The aim of this paper is to demonstrate that pions present an exception to this rule. We shall argue that it is possible to determine the dispersion relation of soft pions (more precisely, its real part) at all temperatures below the chiral phase transition, knowing only equal-time (or static) correlation functions, which, in principle, can be determined on the lattice. It should be emphasized that we do not assume the temperature $`T`$ to be small compared to the chiral phase transition temperature $`T_c`$: we must have $`T<T_c`$, but $`T/T_c`$ is allowed to be of order 1. Our results are valid in the nontrivial regime where neither perturbative QCD nor chiral perturbation theory are reliable. Moreover, our method, with minimal modification, can be applied to any field theory with a broken symmetry, at temperatures below symmetry restoration.

That the dispersion relation of a mode can be expressed fully in terms of static correlation functions is nontrivial, but by no means unprecedented. We recall the sound waves, whose velocity is $`u=(p/ϵ)^{1/2}`$, where $`p`$ and $`ϵ`$ are the pressure and the energy density, respectively. The sound speed, while being a real-time quantity (the pole in the correlator of the energy density $`T^{00}`$), can be determined solely from thermodynamics. The relation between the speed of sound and the thermodynamic functions is an exact consequence of the existence of the hydrodynamic description. A less familiar example is the variational Feynman-Bijl formula which relates the phonon spectrum in superfluid helium to the static density-density correlation function BijlFeynman . The example most closely related to our problem, however, is that of spin waves in antiferromagnets HHPR69 : the velocity of spin waves at any temperature below phase transition is equal to the ratio of the stiffness and the magnetic susceptibility in a direction perpendicular to magnetization. Both quantities can be defined from the static response of the system to external fields. The only difference between QCD and antiferromagnets is that in the former case the symmetry is $`\mathrm{SU}(2)_V\times \mathrm{SU}(2)_A\mathrm{O}(4)`$, which is spontaneously broken to $`\mathrm{SU}(2)_V\mathrm{O}(3)`$, while in ferromagnets $`\mathrm{O}(3)`$ is broken down to $`\mathrm{O}(2)`$ PisarskiWilczek ; RajagopalWilczek .

The paper is constructed as follows. In Sec. II we summarize the findings of the paper. In Sec. III we use simple, but nonrigorous, arguments relying on an effective Lagrangian to understand these results. In Sec. IV the results are derived in a more rigorous way from a set of assumptions about the real-time correlation functions, which comes from hydrodynamics. In Sec. V we show that our result holds for the simplest field-theoretical model of a scalar field theory with broken symmetry. In the Appendix we give a simple derivation of the known result about the dynamical critical exponent $`z`$, and derive the critical scaling of a diffusion coefficient.

## II Summary of results

We claim that, in QCD with two light flavors, at temperatures below the chiral phase transition, the real part of the dispersion relation of sufficiently soft pions is given by the following equation:

$$\omega _𝒑^2=u^2(𝒑^2+m^2).$$
(1)
In this paper the following terminology is used: $`u`$ is the pion velocity (although, strictly speaking, it is the pion velocity only when $`m=0`$), and $`m`$ is the pion screening mass (we shall show that it is the same screening mass as defined on the lattice). The energy of a pion at $`𝒑=\mathrm{𝟎}`$, $`m_p=um`$ is called the pion pole mass.

At finite temperature, the meaning of soft pions may need some clarification. We shall understand Eq. (1) as a statement that the correlators of operators carrying pion quantum numbers have a pole at the frequency with the real part determined by Eq. (1).

The pion velocity $`u`$ is the ratio of two statically measurable quantities, the temperature-dependent pion decay constant $`f^2`$ and the axial isospin susceptibility $`\chi _{I5}`$,

$$u^2=\frac{f^2}{\chi _{I5}}.$$
(2)
The axial isospin susceptibility $`\chi _{I5}`$ can be defined as the second derivative of the pressure with respect to the axial isospin chemical potential (see Sec. III). Equivalently, it can be defined via the static Euclidean correlator of the axial isospin charge densities,

$$\delta ^{ab}\chi _{I5}=\underset{0}{\overset{1/T}{}}𝑑\tau 𝑑𝒙A_0^a(\tau ,𝒙)A_0^b(0,\mathrm{𝟎}),$$
(3)
where

$$A_0^a\overline{\psi }\gamma ^0\gamma ^5\frac{\tau ^a}{2}\psi ,$$
(4)
$`\psi `$ is the quark field, $`a,b=1,2,3`$, $`\tau ^a`$ are isospin Pauli matrices, $`\mathrm{Tr}\tau ^a\tau ^b=2\delta ^{ab}`$, and $`\mathrm{}`$ denotes thermal averaging, which can be taken by evaluating a Euclidean Feynman path integral with appropriate boundary conditions. The quantity $`\chi _{I5}`$ has been considered previously Hansen .

The pion decay constant $`f`$ and the screening mass $`m`$ can be determined from the static Euclidean pion correlators at small momenta, which are predicted to have the form

$$\underset{0}{\overset{1/T}{}}𝑑\tau 𝑑𝒙e^{i𝒒𝒙}\phi ^a(\tau ,𝒙)\phi ^b(0,\mathrm{𝟎})=\frac{1}{f^2}\frac{\delta ^{ab}}{𝒒^2+m^2},$$
(5)
where the scalar field $`\phi ^a`$ is defined as

$$\phi ^a\frac{i\overline{\psi }\gamma ^5\tau ^a\psi }{\overline{\psi }\psi }.$$
(6)
Equation (5) is supposed to be valid when $`|𝒒|m_\sigma `$, where $`m_\sigma `$ is the (temperature-dependent) screening mass of the $`\sigma `$ meson.

The way $`m`$ enters Eq. (5) explains why we term it the pion screening mass. As far as $`f`$ is concerned, the conventional definition of the pion decay constant, $`f_\pi p_\mu =0|A_\mu |\pi (p)`$, cannot be used at finite temperature, since neither the vacuum state $`|0`$ nor the one-pion state $`|\pi `$ allows generalization in thermal media. Both $`\sqrt{\chi _{I5}}`$ and $`f`$ approach $`f_\pi `$ as $`T0`$, so both can be viewed as the generalization of $`f_\pi `$ to finite temperature. We, however, reserve the name “pion decay constant” for the quantity $`f`$ defined in Eq. (5).222See Sec. IV.5 for a discussion related to this point. In contrast to the susceptibility $`\chi _{I5}`$, $`f`$ cannot be defined at temperatures above critical.

We shall also show that the Gell-Mann–Oakes–Renner (GOR) relation can be generalized to finite temperature, and does in fact become two separate relations for the screening and pole masses of the pions,

$$f^2m^2=\chi _{I5}m_p^2=m_q\overline{\psi }\psi .$$
(7)
For simplicity, in this paper we assume $`m_u=m_d=m_q`$.

Equation (2) is the direct analogue of a similar equation for the velocity of spin waves in a quantum antiferromagnet at temperatures below the phase transition HHPR69 . The counterpart of $`f`$ is the stiffness (denoted as $`\rho _s`$ in Ref. HHPR69 ), while $`\chi _{I5}`$ is similar to the magnetic susceptibility in a direction perpendicular to magnetization.

An important point we wish to emphasize is that the relations (1)–(7) are exact in the chiral limit, i.e., the limit when $`m_\pi `$ (vacuum pion mass) and $`𝒑`$ are infinitesimally small, at any temperature in the interval from 0 to $`T_c`$. Although the methods we use to derive them may look similar to chiral perturbation theory, unlike the chiral perturbation theory, we do not consider perturbations around $`m_\pi =0`$, $`T=0`$, treating both $`m_\pi `$ and $`T`$ as small parameters. Rather, we consider perturbations around $`m_\pi =0`$ at a fixed $`T`$, the latter not assumed to be small. Therefore, unlike in chiral perturbation theory, we are not able to calculate the temperature dependence of parameters such as $`\chi `$, $`f`$, $`u`$, or $`\overline{\psi }\psi `$. However, we are able to show that their temperature dependence must be such that relations (1)–(7) hold.

Given these equations, exact in the chiral limit, it is legitimate to ask how useful these relations can be at the physical value of the pion mass $`m_\pi =140`$ MeV. To answer this question quantitatively, a (perhaps lattice) calculation of the $`T`$ and $`m_\pi `$ dependence is needed, which is beyond the scope of this paper. However, we shall attempt to give a semiquantitative answer to this question.

It is easier to begin with $`T=0`$. In this case formula (1) becomes trivial: the dispersion relation is exact, simply by Lorentz invariance, for any $`𝒑`$ and $`m=m_\pi `$. It is a straightforward exercise to show, using PCAC (partial conservation of axial vector current), that to leading order, i.e., $`𝒪(m_\pi ^2)`$, the correlator in Eq. (3) whose value at $`m_\pi =0`$ is $`\chi _{I5}`$, does not depend on $`m_\pi `$. Thus one can expect that measuring the correlator (3) even at the physical $`m_\pi `$ (or higher, as is typical in a lattice calculation due to the high price of simulating with light quark masses) is not significantly different from $`\chi _{I5}`$. In order to extract the value of $`f`$ from the relation (5) reliably, one requires that higher mass states in the channel with pion quantum numbers do not significantly contaminate the exponential falloff of the correlator in coordinate space. Since the masses of such states are typically above 1 GeV, the exponential tail contribution from the light state of mass $`m_\pi =140`$ MeV should be easy to separate.

At nonzero $`T`$ the region of applicability of Eq. (2) is limited by terms of higher powers in $`𝒑`$. For simplicity, in the exact chiral limit $`m_\pi =0`$ the dispersion relation takes the form $`\omega =|𝒑|\frac{i}{2}D^{}𝒑^2+\mathrm{}`$, where $`D^{}`$ is a temperature-dependent parameter (see Sec. IV and the Appendix). The condition on momenta which is required to neglect nonlinearity in the dispersion relation is $`|𝒑|1/D^{}`$. Calculation of the diffusion constant $`D^{}`$ is a challenging task. Taking an estimate from Ref. Smilga , $`D^{}=CT^3/f_\pi ^4`$, with a numerically rather small $`C0.1`$, we can conclude that even for $`Tf_\pi `$ at momenta of order $`|𝒑|100`$ MeV the nonlinearity in the dispersion relation is still small.

As $`T=T_c`$ the dispersion relation is essentially nonlinear, which manifests itself in the divergence of $`D^{}`$ as $`TT_c`$ (see the Appendix). As $`TT_c`$ the maximum momentum at which the dispersion relation can be considered linear (in the chiral limit) decreases and vanishes at $`T_c`$. Although the power with which the width of the linearity window shrinks to zero can be determined (by extension of arguments given in the Appendix), the pre-exponent is unknown, and a nonperturbative calculation is, in principle, required to determine quantitatively the size of the nonlinearity in the dispersion relation at a given $`T`$ and $`|𝒑|`$. Such a calculation is beyond the scope of the paper.

Equations (1)–(7) were used in Ref. critical\_pions  to extract information about the critical behavior of the pion velocity and masses near the critical temperature. A brief derivation of Eqs. (1)–(7) was also sketched in Ref. critical\_pions . We present a more extended version of this derivation in Sec. III.

In Sec. IV we provide a new derivation of the relations (1)–(7) using the operator approach based on hydrodynamic equations. For this purpose we shall need expressions for $`\chi _{I5}`$ and $`f`$ in terms of equal-time rather than static (zero frequency) correlators. These relations follow from definitions (3) and (5) in the exact chiral limit ($`m_\pi =0`$) or when the temperature is high enough ($`Tm_\pi `$):

$$\delta ^{ab}\chi _{I5}=\frac{1}{T}𝑑𝒙A_0^a(t,𝒙)A_0^b(t,\mathrm{𝟎}).$$
(8)
In the chiral limit, the equality of the static correlator \[Eq. (3)\] and the equal-time correlator \[Eq. (8)\] of $`A_0^a`$ is a consequence of the conservation of the axial isospin charge. Similarly,

$$\frac{1}{T}𝑑𝒙e^{i𝒒𝒙}\phi ^a(t,𝒙)\phi ^b(t,\mathrm{𝟎})=\frac{1}{f^2}\frac{\delta ^{ab}}{𝒒^2+m^2}.$$
(9)
Note that, in contrast to Eqs. (3) and (5), Eqs. (8) and (9) require an additional condition $`Tm_\pi `$. This condition is, however, not so dramatic, if one recalls that the relavant energy scale of order $`T`$ is the lowest nonzero Matsubara frequency, $`2\pi T`$. Thus, in practice, one requires $`2\pi Tm_\pi `$, which is satisfied reasonably well for $`T`$ of order $`T_c160`$ MeV and the physical pion mass.

However, a much stronger condition is required for the applicability of the hydrodynamic description: $`|𝒑|`$ (and $`m_\pi `$) must be smaller than the typical collision rate $`\tau ^1T^5/f_\pi ^4`$ (according to Ref. Goity ). At small $`Tf_\pi `$ this condition is much stronger than the conditions that are required for the effective Lagrangian derivation to hold (typically $`|𝒑|,m_\pi m_\sigma `$ — see Sec. III). However, both derivations apply in the required regime infinitesimally close to the chiral limit. The fact that the hydrodynamic derivation has a smaller validity range, however, does not mean that the result is not valid outside this range. As an example, consider the case of $`T=0`$, where the results (1), (2) hold trivially, while hydrodynamics does not apply at all. A less trivial example is presented in Sec. V, where explicit calculation verifies Eq. (2). This weak-coupling calculation, however, does not require any notion of hydrodynamics.

The advantage of the hydrodynamic approach is that it allows us to consider properly the effects of the dissipation, which are inherently beyond the Lagrangian approach. Before proceeding to the new operator derivation (Sec. IV) we review the derivation based on the effective Lagrangian approach.

## III Effective Lagrangian approach

### III.1 Axial isospin susceptibility and $`f_\pi `$ at $`T=0`$

Since an important role in our analysis is played by the axial isospin susceptibility, we first consider this quantity at zero temperature and its relation to the pion decay constant. The quark part of the QCD Lagrangian at finite axial isospin chemical potential $`\mu _{I5}`$ is given by

$$_{\mathrm{quark}}=i\overline{\psi }\gamma ^\mu D_\mu \psi m_q\overline{\psi }\psi +\mu _{I5}\overline{\psi }\gamma _0\gamma _5\frac{\tau _3}{2}\psi ,$$
(10)
where $`D_\mu `$ is the color covariant derivative. This chemical potential $`\mu _{I5}`$ is coupled to the axial isospin charge $`A_0^3`$ \[cf. Eq. (4)\], which generates the $`\mathrm{SU}(2)_A`$ part of the $`\mathrm{SU}(2)_V\times \mathrm{SU}(2)_A`$ chiral symmetry.

The response of the QCD vacuum to $`\mu _{I5}`$ can be found from the effective chiral Lagrangian. The latter, to lowest order of momenta, masses, and chemical potential, is completely fixed by the chiral symmetries and is given by

$$_{\mathrm{eff}}=\frac{f_\pi ^2}{4}\mathrm{Tr}^\nu \mathrm{\Sigma }_\nu \mathrm{\Sigma }^{}+\frac{f_\pi ^2m_\pi ^2}{2}\mathrm{Re}\mathrm{Tr}\mathrm{\Sigma },$$
(11)
where $`\mathrm{\Sigma }`$ is an $`\mathrm{SU}(2)`$ matrix whose phases describe the pions, $`\mathrm{\Sigma }=e^{i\tau ^a\pi ^a/f_\pi }`$, and $``$ denotes the covariant derivative, which is defined as

$$_0\mathrm{\Sigma }=_0\mathrm{\Sigma }\frac{i}{2}\mu _{I5}(\tau _3\mathrm{\Sigma }+\mathrm{\Sigma }\tau _3),_i\mathrm{\Sigma }=_i\mathrm{\Sigma },i=1,2,3.$$
(12)
At $`\mu _{I5}=0`$, the Lagrangian (11) is the standard chiral Lagrangian with two phenomenologically determined constants $`f_\pi `$ and $`m_\pi `$. The way $`\mu _{I5}`$ enters the effective description is completely fixed by symmetries to lowest order. This can be seen by promoting the $`\mathrm{SU}(2)_A`$ symmetry to a local symmetry and treating $`\mu _{I5}`$ as the time component of the $`\mathrm{SU}(2)_A`$ vector potential mu2 .

The Lagrangian (11) and its derivation is analogous to the case of the effective Lagrangian at finite (vector) isospin chemical potential $`\mu _I`$ considered in Ref. mui . A significant difference between the cases studied here and in Ref. mui  is that the QCD vacuum breaks the $`\mathrm{SU}(2)_A`$ (axial isospin) symmetry spontaneously, but remains symmetric under the $`\mathrm{SU}(2)_V`$ (vector isospin) symmetry. It is important to note, however, that the $`\mathrm{SU}(2)_A`$ is a symmetry of the Lagrangian (at $`m_q=0`$), as good as the $`\mathrm{SU}(2)_V`$. In particular, the axial isospin current $`A_\mu ^a`$ is conserved in the chiral limit. Thus, it is entirely legitimate to consider the theory at finite $`\mu _{I5}`$ and use symmetry arguments to fix the $`\mu _{I5}`$ dependence of the effective Lagrangian.

The vacuum energy density depends nontrivially on $`\mu _{I5}`$ already for arbitrarily small $`\mu _{I5}`$. (This is in contrast to the case of the isospin chemical potential $`\mu _I`$, where $`\mu _I`$ needs to be larger than a threshold equal to $`m_\pi `$ in order to change the ground state.) The isospin axial susceptibility, at $`\mu _{I5}=0`$, is easy to determine using the effective Lagrangian (11):

$$\chi _{I5}\frac{^2_{\mathrm{vac}}}{\mu _{I5}^2}|_{\mu _{I5}=0}=\frac{^2_{\mathrm{eff}}}{\mu _{I5}^2}|_{\mathrm{\Sigma }=1}=f_\pi ^2.$$
(13)
This result is of potential importance for lattice QCD calculations because, in principle, it allows one to determine $`f_\pi `$ directly by measuring the axial isospin susceptibility. To our knowledge, this has not been done on the lattice for temperatures below the chiral phase transition. The isospin susceptibility (as well as the isoscalar, i.e., baryon number susceptibility) measurement has been done using staggered fermions in gottlieb . In the formulation used in Ref. gottlieb , introducing the axial isospin chemical potential would require replacing $`\mathrm{exp}(\mu a)`$ on the links (where $`a`$ is the lattice spacing) with $`\mathrm{exp}(\mu a\zeta _x)`$, where $`\zeta _x`$ is the usual staggered factor $`\zeta _x(1)^{x_1+x_2+x_3+x_4}`$ — the representation of $`\gamma _5`$ in the staggered fermion action.

### III.2 Pion velocity

To obtain Eq. (2), we expand the previous discussion to nonzero temperature. We first presume that the dynamics of the pions is described by some effective Lagrangian $`_{\mathrm{eff}}`$. Strictly speaking, this is not correct since dissipative effects cannot be included in the effective Lagrangian. We can expect to recover the correct answers if, in the infrared, the pion thermal width is negligible compared to its energy. This has been seen in explicit calculations at low $`T`$, where chiral perturbation theory can be used Shuryak ; Goity . We furthermore assume that this Lagrangian is local and can be expanded in powers of momenta. To lowest order, the Lagrangian is fixed by symmetries up to three coefficients, $`f_t`$, $`f_s`$, and $`m`$,

$$_{\mathrm{eff}}=\frac{f_t^2}{4}\mathrm{Tr}_0\mathrm{\Sigma }_0\mathrm{\Sigma }^{}\frac{f_s^2}{4}\mathrm{Tr}_i\mathrm{\Sigma }_i\mathrm{\Sigma }^{}+\frac{m^2f_s^2}{2}\mathrm{Re}\mathrm{Tr}\mathrm{\Sigma }.$$
(14)
Due to the lack of Lorentz invariance, $`f_t^2`$ and $`f_s^2`$ are independent parameters. The covariant derivative $`_0`$ is the same as defined in Eq. (12). The dispersion relation following from this Lagrangian has the form (1) where

$$u=\frac{f_s}{f_t}.$$
(15)
At zero temperature $`f_t=f_s=f_\pi `$, and the pion velocity $`u`$ is equal to the speed of light. We now show that, at finite temperature, all three parameters $`f_t`$, $`f_s`$, and $`m`$ can be determined from equal-time (or static) correlation functions.

Repeating the same argument as in Sec. III.1, we can show that $`f_t`$ is related to the axial isospin susceptibility,

$$\chi _{I5}=f_t^2,$$
(16)
but now $`\chi _{I5}`$ is defined as the susceptibility at finite temperature. Note that, as a susceptibility with respect to a conserved charge, $`\chi _{I5}`$ is free of ultraviolet divergences.

We now need a prescription to compute $`f_s`$ and $`m`$ from static correlation functions. In order to make the connection, we generalize the mass term in Eq. (10),

$$_{\mathrm{quark}}=i\overline{\psi }\gamma ^\mu D_\mu \psi (\overline{\psi }_LM\psi _R+\mathrm{H}.\mathrm{c}.)+\mu _{I5}A_0^3,$$
(17)
and regard $`M`$ as an external field, $`M=M(x)`$. The Lagrangian (17) possesses the symmetry $`\psi _LL\psi _L`$, $`\psi _RR\psi _R`$, $`MLMR^{}`$, where $`L,R\mathrm{SU}(2)`$. If $`M(x)`$ is a slowly varying function of $`x`$, its effect can be captured in the effective chiral Lagrangian. The requirement that the effective Lagrangian preserves this symmetry fixes the form of its mass term,

$$_{\mathrm{eff}}=\frac{f_t^2}{4}\mathrm{Tr}_0\mathrm{\Sigma }_0\mathrm{\Sigma }^{}\frac{f_s^2}{4}\mathrm{Tr}_i\mathrm{\Sigma }_i\mathrm{\Sigma }^{}\frac{1}{2}\overline{\psi }\psi \mathrm{Re}\mathrm{Tr}M^{}\mathrm{\Sigma }.$$
(18)
We shall limit ourselves to a particular ansatz of the external field $`M(x)`$:

$$M(x)=m_qe^{i\alpha ^a(x)\tau ^a}.$$
(19)
The second derivative of the partition function with respect to $`\alpha ^a`$ can be computed in both the microscopic theory (17) and in the effective theory (18). In the microscopic theory, we find

$$\frac{\delta ^2\mathrm{ln}𝒵}{\delta \alpha ^a(x)\delta \alpha ^b(0)}=m_q^2\overline{\psi }\psi ^2\phi ^a(x)\phi ^b(0),$$
(20)
where $`\phi ^a`$ is defined in Eq. (6). On the other hand, from the the effective Lagrangian (11) we find

$$\frac{\delta ^2\mathrm{ln}𝒵}{\delta \alpha ^a(x)\delta \alpha ^b(0)}=m_q^2\overline{\psi }\psi ^2\varphi ^a(x)\varphi ^b(0),\text{ where }\varphi ^a(x)\mathrm{Re}\mathrm{Tr}i\tau ^a\mathrm{\Sigma }(x)/2.$$
(21)
Comparing Eqs. (20) and (21), we see that correlation functions of $`\phi ^a(x)`$ defined in the microscopic theory and of $`\varphi ^a(x)`$ defined in the effective theory are equal. This equality should hold for small momenta when the effective theory is applicable.

The correlation function of $`\varphi ^a(x)`$, on the other hand, can be calculated by expanding the effective Lagrangian to second order in $`\varphi ^a`$. The Matsubara propagator of $`\varphi ^a`$ is

$$\underset{0}{\overset{1/T}{}}𝑑\tau 𝑑𝒙e^{iqx}\varphi ^a(x)\varphi ^b(0)=\frac{\delta ^{ab}}{f_t^2q_0^2+f_s^2(𝒒^2+m^2)},q_0=2\pi Tn.$$
(22)
This propagator, for $`q_0=0`$ and small $`𝒒`$, should be equal to the propagator of $`\phi ^a`$. This establishes Eq. (5), with the identification $`f=f_s`$. Together with Eqs (15) and (16), this is our result for the pion dispersion relation. Furthermore, it is natural to assume that, for small $`𝒒`$, the dynamics of $`\phi ^a`$ is slow, so at high enough temperature ($`Tm_\pi `$) one can regard $`\phi ^a(\tau ,𝒙)`$ as independent of $`\tau `$. In this case the left hand side of Eq. (22) is equal to that of Eq. (9).

It is instructive to write Eq. (9) in coordinate space,

$$\phi ^a(0,𝒙)\phi ^b(0,\mathrm{𝟎})=\frac{T}{4\pi f_s^2}\frac{e^{m|𝒙|}}{|𝒙|}(|𝒙|T^1,m_\sigma ^1).$$
(23)
We see that measuring, in the microscopic theory, the large-distance equal-time correlation function of the $`\phi ^a`$ defined in Eq. (6) we can extract the screening mass $`m`$ and the decay constant $`f`$. Combined with the determination of the susceptibility $`\chi _{I5}`$ and Eq. (16), the dispersion relation of soft pions is now completely known.

Below we shall provide a more systematic proof of these relationships using an operator approach, with crucial inputs from the hydrodynamic theory. We shall also demonstrate the validity of the relation (2) in an explicit lowest order perturbative calculation in the linear sigma model.

## IV Hydrodynamic (operator) approach

From a modern perspective, hydrodynamics is an effective theory operating at sufficiently large distance and time scales. (By “sufficiently large” normally we mean scales larger than the mean free path, or the relaxation time.) As such, it is the most suitable framework to discuss low-energy degrees of freedom (like pions) at finite temperatures. This modern point of view, as opposed to the view of hydrodynamics as a purely phenomenological description, has been in existence for a long time hydro\_corr . From this philosophy, it is not surprising that the same hydrodynamic theory describes systems with very different microscopic dynamics, and that systems with different symmetries (like normal fluids and superfluids) correspond to different hydrodynamic theories. Similarly to the effective Lagrangian, hydrodynamic equations can also be viewed as a particular way to satisfy Ward identities at finite temperatures Yaffe . In our case, hydrodynamics is nontrivial due to the chiral symmetry breaking Sonhydro . We will derive the constraints placed by hydrodynamics on the dynamics in our problem, and show how our results on the pion dispersion relation follow from there. Our treatment is similar, to some degree, to that of Ref. HHPR69 .

### IV.1 Basic assumptions of hydrodynamics

We shall assume that, as one goes sufficiently far into the infrared, the dynamics of any interacting finite-temperature system can be described in terms of a finite number of fields, which will be called hydrodynamic variables. To be relevant in the infrared, the fluctuations of these fields (either thermal fluctuations, or those due to external sources) should relax arbitrarily slowly. This requirement eliminates most of the degrees of freedom, which typically relax during a relaxation time determined by the microscopic dynamics. However, the following fields are obvious hydrodynamic variables:

* The densities of conserved quantities, including the energy density $`T^{00}`$, the momentum density $`T^{0i}`$, and the densities of conserved global charges (i.e., the zeroth components of conserved currents). These fields cannot relax quickly because of the conservation laws. A configuration where charges fluctuate over a length scale $`L`$ much larger than the mean free path can relax only by diffusion, which takes place over a time proportional to $`L^2`$. The relaxation time diverges with the wavelength of the perturbation.
* The phases of the condensates which break global symmetries. At zero temperature the fluctuations of such phases correspond to Goldstone bosons, whose energy can be arbitrarily small. At finite temperature below symmetry restoration, one should also expect the long-wavelength fluctuations of the condensate phases to relax slowly.
* Near the critical temperatures of second-order phase transitions, the order parameters themselves (not just the phases) should be considered hydrodynamic variables Hohenberg . For example, it is believed that, in QCD with two massless flavors, the chiral phase transition is of the second order, where the $`\sigma `$ meson becomes degenerate with the pions. Near $`T_c`$, hence, $`\sigma `$ should be included in the hydrodynamic description. In contrast to the fields in the categories (i) and (ii), the rate of relaxation of order parameters is controlled by the closeness of $`T`$ to $`T_c`$, but not by the wavelength of perturbations.333We shall not consider the possibility of Abelian gauge fields in the Coulomb phase (as opposed to the Higgs phase) present in magnetohydrodynamics.

If one makes a further assumption that the fields listed in (i)–(iii) exhaust all slowly relaxing ones, then the set of hydrodynamic degrees of freedom is completely fixed once all symmetries of the theory and the pattern of symmetry breaking at the given temperature are known. In this paper we shall limit ourselves to temperatures far away from any second order phase transition, so order parameters are excluded from hydrodynamics. This set of hydrodynamic variables always contains the energy and momentum densities $`T^{00}`$ and $`T^{0i}`$. For QCD below the chiral phase transition, one has, in addition, the densities of baryon, isospin, and axial isospin charges, and the phases of the chiral condensate.

After identifying the hydrodynamic degrees of freedom, one can proceed in various ways. One can ask about the the equations of motion that the hydrodynamic variables obey. For our case, it is a rather nontrivial task, because of the multitude of the fields involved. In Ref. Sonhydro , the Poisson-bracket technique is used to derive the dissipationless hydrodynamic equations. One finds a system of fully nonlinear coupled differential equations, generalizing the equations of fluid dynamics and equations of motion of the nonlinear sigma model. However, the connection of this procedure to the fundamental (microscopic) field theory has not been made, and the physical meaning of several temperature-dependent parameters appearing in the final equations (called $`f_t`$, $`f_s`$, and $`f_v`$ in Ref. Sonhydro ) is not at all clear.

Alternatively, one can ask the question: what are the constraints that hydrodynamics places on the correlation functions of the hydrodynamic variables? It is clear from the discussion above that the real-time correlators of the hydrodynamic variables \[i.e., $`𝒪(t,𝒙)𝒪(0,\mathrm{𝟎})`$, where $`𝒪`$ is the variable\] are long range, i.e., have power-law (but not exponential) falloff, at least in the timelike regime $`t\mathrm{}`$, $`𝒙=\text{fixed}`$.444The hydrodynamic correlators are not required to have power-law decay in the regime $`t=\text{fixed}`$, $`𝒙\mathrm{}`$. It is less trivial to decide about the correlators of the operators not belonging to the set of hydrodynamic variables. To this end, the following additional assumption is made.

* All local operators can be expressed as local functions of the hydrodynamic operators and their spatial derivatives, up to corrections which have short-ranged correlations that go to zero exponentially when either temporal or spatial separation goes to infinity ($`t`$ or $`𝒙\mathrm{}`$). (These short-range parts correspond to the “noises” in the hydrodynamic theory.)

Physically, these assumptions mean that the dynamics in the infrared can be described in terms of hydrodynamic variables only, and is equivalent to the assumption of local thermodynamic equilibrium: the values of all variables are determined by specifying a few. The assumptions (i)–(iv) form the starting point of our construction of the hydrodynamic equations.

### IV.2 Linearized hydrodynamics for soft pions

To be less abstract, let us consider QCD in the chiral limit, at a finite temperature below the chiral phase transitions, and zero chemical potentials. The hydrodynamic operators are the energy density $`T^{00}`$, the momentum density $`T^{0i}`$, the baryon density $`\overline{\psi }\gamma ^0\psi `$, the densities of vector and axial isospin charges,

$$V_0^a=\overline{\psi }\gamma ^0\frac{\tau ^a}{2}\psi ,A_0^a=\overline{\psi }\gamma ^0\gamma ^5\frac{\tau ^a}{2}\psi ,$$
(24)
and the pion field, defined as Eq. (6),

$$\phi ^a=\frac{i\overline{\psi }\gamma ^5\tau ^a\psi }{\overline{\psi }\psi }.$$
(25)
As we shall see below, to linear order the dynamics of $`A_0^a`$ and $`\phi ^a`$ is decoupled from other modes.

Let us consider the equation for $`A_0^a`$. Infinitesimally close to the chiral limit, we can derive from the QCD Lagrangian the familiar PCAC relation:

$$_\mu A^{a\mu }=m_q\overline{\psi }\psi \phi ^a.$$
(26)
The left hand side contains, in addition to $`A_0^a`$, the spatial components of the axial current $`A^{ai}`$, which are not hydrodynamic variables. According to the assumption (iv), we can express $`A^{ai}`$ as local functions of the hydrodynamic operators and their spatial derivatives, plus a short-ranged part. If we work to leading order in the power of the fields (which means the linear order), the only ones suitable are $`\phi ^a`$ and $`A_0^a`$ which are parity odd and isovectors. The spatial index in $`A^{ai}`$ forces one to have at least one spatial derivative. Thus, to leading orders in fields and derivatives, the only terms consistent with symmetries are

$$A^{ai}=f^2_i\phi ^aD_iA_0^a\xi ^{ai},$$
(27)
where $`f^2`$ and $`D`$ are coefficients depending on the temperature, and $`\xi _i^a`$ is the short-range (“noise”) part of $`A^{ai}`$. (At this step, we have not yet related $`f^2`$ to the static correlation functions. We will do it later on.) Equation (26) now takes the form

$$_0A_0^a=f^2^2\phi ^a+m_q\overline{\psi }\psi \phi ^a+D^2A_0^a+_i\xi ^{ai}.$$
(28)
The parameter $`D`$ can be interpreted as the diffusion coefficient for $`A_0^a`$; however, Eq. (30a) is more complex than a diffusion equation. Now let us discuss the equation for $`\phi ^a`$. The time derivative of $`\phi ^a`$, not being a hydrodynamic operator, can be expanded as

$$_0\phi ^a=\frac{1}{\chi }A_0^a\kappa _1^{}\phi ^a+\kappa _2^2\phi ^a+\eta ^a,$$
(29)
where $`\chi `$, $`\kappa _1^{}`$, and $`\kappa _2`$ are again coefficients dependent on temperature, and $`\eta `$ is a short-range noise. We have kept the term with the second derivative of $`\phi ^a`$ as well as the term with no derivative. The reason for doing so is that $`\kappa _1^{}`$ is suppressed by the quark masses. Indeed, in the chiral limit the state with $`\phi ^a0`$, $`A_0^a=0`$ can be another vacuum which stays unchanged with time. This is consistent with Eq. (30b) only if $`\kappa _1^{}=0`$ in this limit, thus at small quark masses $`\kappa _1^{}`$ is small. For this reason we keep the $`^2\phi ^a`$ term.

Introducing the parameter $`m^2`$ defined so that $`f^2m^2=m_q\overline{\psi }\psi `$ (which we still have to relate to static correlators), and $`\kappa _1`$ so that $`\kappa _1^{}=\kappa _1m^2`$, Eqs. (28) and (29) can be written as

$`_0A_0^a`$ $`=`$ $`f^2(^2m^2)\phi ^a+D^2A_0^a+_i\xi ^{ai},`$ (30a)
$`_0\phi ^a`$ $`=`$ $`{\displaystyle \frac{1}{\chi }}A_0^a+(\kappa _2^2\kappa _1m^2)\phi ^a+\eta ^a.`$ (30b)
Equations (30) are the linearized hydrodynamic equations governing the evolution of $`A_0^a`$ and $`\phi ^a`$. The correlation functions of $`A_0^a`$ and $`\phi ^a`$ can be found if one knows the correlators of the “noises” $`\xi ^{ai}`$ and $`\eta ^a`$. By construction, these fields have only short-range correlations, so, if one is interested only in the dynamics at large distance and time scales, these correlations can be replaced by delta functions. Isospin symmetry and rotational invariance require the correlators to be of the following forms:

$`\xi ^{ai}(x)\xi ^{bj}(y)`$ $`=`$ $`F_\xi \delta ^{ab}\delta ^{ij}\delta ^4(xy),`$ (31a)
$`\eta ^a(x)\eta ^b(y)`$ $`=`$ $`F_\eta \delta ^{ab}\delta ^4(xy),`$ (31b)
$`\xi ^{ai}(x)\eta ^b(y)`$ $`=`$ $`0.`$ (31c)
up to corrections proportional to derivatives of $`\delta ^4(xy)`$ which will be neglected since they are of higher order in momentum. Equations (30) and (31) completely determine the hydrodynamics of soft pions, to the linearized order.

### IV.3 Relation to static correlators

Our next task is to relate, as much as possible, the parameters appearing in Eqs. (30) and (31) with the equal-time correlators of $`A_0^a`$ and $`\phi ^a`$.

First, we multiply Eq. (30b) by $`A^b`$, taken at the same time moment, and integrate over space. One finds

$$\begin{array}{cc}\hfill 𝑑𝒙\dot{\phi }^a(t,𝒙)A_0^b(t,\mathrm{𝟎})=& d𝒙[\frac{1}{\chi }A_0^a(t,𝒙)A_0^b(t,\mathrm{𝟎})\kappa _1m^2\phi ^a(t,𝒙)A_0^b(t,\mathrm{𝟎})\hfill \\ & +\eta ^a(t,𝒙)A_0^b(t,\mathrm{𝟎})].\hfill \end{array}$$
(32)
In the right hand side, the second term is proportional to $`m^2`$, which is small in the chiral limit, and hence can be neglected. The last term in the integrand will be shown later to vanish,

$$\eta ^a(t,𝒙)A_0^b(t,\mathrm{𝟎})=0.$$
(33)
Equation (33) is trivial if understood as $`\eta ^a(t+\epsilon ,𝒙)A_0^b(t,\mathrm{𝟎})=0`$: $`A_0^b`$ cannot be correlated with the noise $`\eta `$ in the future. What is somewhat less trivial (and will be checked a posteriori) is that $`\eta ^a(t\epsilon ,𝒙)A_0^b(t,\mathrm{𝟎})`$ also vanishes: the equal-time correlator of $`\eta `$ and $`A_0`$ does not depend on how the equal-time limit is taken. Therefore, only the first term survives; and by definition of $`\chi _{I5}`$ it is equal to

$$𝑑𝒙\frac{1}{\chi }A_0^a(t,𝒙)A_0^b(t,\mathrm{𝟎})=\frac{T\chi _{I5}}{\chi }.$$
(34)
On the other hand, the left hand side of Eq. (32) can be computed explicitly. To this end, we write

$$\dot{\phi }^a(𝒙)A_0^b(\mathrm{𝟎})=i\mathrm{Tr}\{e^{\beta H}[H,\phi ^a(𝒙)]A_0^b(\mathrm{𝟎})\}.$$
(35)
(we drop the time variable $`t`$ which is an argument of all operators). We now make use of the following expansion:

$$[e^{\beta H},\phi ^a]=\beta e^{\beta H}[H,\phi ^a]\frac{\beta ^2}{2}e^{\beta H}[H,[H,\phi ^a]]+\mathrm{}.$$
(36)
The expansion parameter here is $`\beta q_0q_0/T`$ where $`q_0`$ is the frequency of variation of $`\phi ^a`$. Since we are dealing with the low-frequency modes in $`\phi ^a`$, we can ignore all terms beyond the first in Eq. (36). Equation (35) now becomes

$$\dot{\phi }^a(𝒙)A_0^b(\mathrm{𝟎})=iT\mathrm{Tr}\{[e^{\beta H},\phi ^a(𝒙)]A_0^b(\mathrm{𝟎})\}=iT[\phi ^a(𝒙),A_0^b(\mathrm{𝟎})]=T\delta ^{ab}\delta ^3(𝒙).$$
(37)
In the last transformation, we make use of the commutation relation

$$[\phi ^a(𝒙),A_0^b(\mathrm{𝟎})]=i\delta ^{ab}\delta ^3(𝒙)\frac{\overline{\psi }\psi }{\overline{\psi }\psi }.$$
(38)
Comparing to Eq. (34), we find

$$\chi =\chi _{I5}.$$
(39)
Thus, we show that $`\chi `$ is the axial isospin susceptibility of the system, which is a static quantity. The proof we have just presented is similar to that of the equipartition theorem in statistical mechanics.

Analogously, we can show that

$$𝑑𝒙e^{i𝒒𝒙}\dot{A}_0^a(𝒙)\phi ^b(\mathrm{𝟎})=T\delta ^{ab}.$$
(40)
On the other hand, from Eq. (30a)

$$\begin{array}{cc}\hfill 𝑑𝒙e^{i𝒒𝒙}\dot{A}_0^a(𝒙)\phi ^b(\mathrm{𝟎})=& f^2(𝒒^2+m^2)𝑑𝒙e^{i𝒒𝒙}\phi ^a(𝒙)\phi ^b(\mathrm{𝟎})\hfill \\ & D𝒒^2𝑑𝒙e^{i𝒒𝒙}A_0^a(𝒙)\phi ^b(\mathrm{𝟎}),\hfill \end{array}$$
(41)
where we dropped the term proportional to

$$_i\xi _i^a(t,𝒙)\phi ^b(t,\mathrm{𝟎})$$
(42)
as it will be shown to vanish in the same manner as $`\eta ^aA_0^b`$ in Eq. (32). Moreover, $`Dq^2A_0^a\phi ^b`$ can be neglected compared to $`\dot{A}_0^a\phi ^bq_0A_0^a\phi ^b`$, since we expect the pion to have a linear dispersion relation and, for small enough momenta, $`q_0Dq^2`$. Equating Eqs. (40) and (41), we find that equal-time correlators of the pion field must have the form of a Yukawa potential,

$$𝑑𝒙e^{i𝒒𝒙}\phi ^a(𝒙)\phi ^b(\mathrm{𝟎})=\frac{T}{f^2}\frac{\delta ^{ab}}{𝒒^2+m^2}.$$
(43)
Thus, the parameters $`f^2`$ and $`m^2`$ appearing in Eq. (30a) are the same ones defined in Eq. (9) via the equal-time correlator of $`\phi ^a`$.

The parameters $`D`$, $`\kappa _1`$, $`\kappa _2`$, $`F_\xi `$, and $`F_\eta `$ cannot be expressed individually in terms of the equal-time correlators, but some relations between them will be derived below.

### IV.4 Hydrodynamic correlation functions

From Eqs. (30) and (31) one can easily compute the real-time correlators of $`A_0^a`$ and $`\phi ^a`$:

$`{\displaystyle d^4xe^{iqx}A_0^a(x)A_0^b(0)}`$ $`=`$ $`\delta ^{ab}{\displaystyle \frac{q_0^2𝒒^2F_\xi +\chi _{I5}^2\omega _𝒒^4F_\eta }{[(q_0\omega _𝒒)^2+\frac{1}{4}\mathrm{\Gamma }_𝒒^2][(q_0+\omega _𝒒)^2+\frac{1}{4}\mathrm{\Gamma }_𝒒^2]}},`$ (44a)
$`{\displaystyle d^4xe^{iqx}\phi ^a(x)\phi ^b(0)}`$ $`=`$ $`\delta ^{ab}{\displaystyle \frac{\chi _{I5}^2𝒒^2F_\xi +q_0^2F_\eta }{[(q_0\omega _𝒒)^2+\frac{1}{4}\mathrm{\Gamma }_𝒒^2][(q_0+\omega _𝒒)^2+\frac{1}{4}\mathrm{\Gamma }_𝒒^2]}},`$ (44b)
where

$`\omega _𝒒^2`$ $`=`$ $`{\displaystyle \frac{f^2}{\chi _{I5}}}(𝒒^2+m^2),`$ (45a)
$`\mathrm{\Gamma }_𝒒`$ $`=`$ $`\kappa _1m^2+(D+\kappa _2)𝒒^2.`$ (45b)
The correlation functions (44) peak around $`q_0\pm \omega _𝒒`$, and the width of the peaks $`\mathrm{\Gamma }_𝒒`$ is much smaller than $`\omega _𝒒`$ at small $`𝒒`$. The correlators have poles corresponding to the pion collective excitations with the dispersion relation $`q_0=\omega _𝒒\frac{i}{2}\mathrm{\Gamma }_𝒒`$. We note, moreover, that all parameters entering the real part of the dispersion relation (45a) can be determined from the static correlation functions.

It is now possible to explicitly check that the noise-field correlators appearing in Eqs. (33) and (42) vanish. Indeed, by using Eqs. (30) and (31) we find

$`\eta (t,𝒙)A_0(t^{},\mathrm{𝟎})`$ $`=`$ $`\delta ^{ab}{\displaystyle \frac{d^4q}{(2\pi )^4}e^{iq_0(tt^{})+i𝒒𝒙}\frac{\chi _{I5}\omega _𝒒^2F_\eta }{q_0^2\omega _𝒒^2+iq_0\mathrm{\Gamma }_𝒒}},`$ (46a)
$`_i\xi _i(t,𝒙)\phi (t^{},\mathrm{𝟎})`$ $`=`$ $`\delta ^{ab}{\displaystyle \frac{d^4q}{(2\pi )^4}e^{iq_0(tt^{})+i𝒒𝒙}\frac{\chi _{I5}^1𝒒^2F_\xi }{q_0^2\omega _𝒒^2+iq_0\mathrm{\Gamma }_𝒒}}.`$ (46b)
The integrands in Eqs. (46) have two poles, both located in the upper half plane: $`q_0=\pm \omega _𝒒+\frac{i}{2}\mathrm{\Gamma }_q`$. When $`t>t^{}`$, when taking integrals over $`q_0`$, one can close the contour in the lower half plane, so the integrals are obviously zero, which corresponds to our previous remark that the fields cannot correlate with the noises in the future. If one takes $`tt^{}`$ from below, the integrals are also zero since they are equal to the sums of the residues of the integrands (which are zero because these functions behave as $`q_0^2`$ at large $`q_0`$.) Therefore, we have checked that the equal-time noise-field correlators (33) and (42) indeed vanish.

One can find the relations between the amplitude of the noise correlators, $`F_\xi `$ and $`F_\eta `$, and the parameters characterizing the damping in Eq. (45b). One integrates Eqs. (44) over $`q_0`$ and obtains the equal-time correlation functions

$`{\displaystyle 𝑑𝒙e^{i𝒒𝒙}A_0^a(t,𝒙)A_0^b(t,\mathrm{𝟎})}`$ $`=`$ $`\delta ^{ab}{\displaystyle \frac{𝒒^2F_\xi +\chi _{I5}^2\omega _𝒒^2F_\eta }{2\mathrm{\Gamma }_𝒒}},`$ (47a)
$`{\displaystyle 𝑑𝒙e^{i𝒒𝒙}\phi ^a(t,𝒙)\phi ^b(t,\mathrm{𝟎})}`$ $`=`$ $`\delta ^{ab}{\displaystyle \frac{𝒒^2F_\xi +\chi _{I5}^2\omega _𝒒^2F_\eta }{2\chi _{I5}^2\omega _𝒒^2\mathrm{\Gamma }_𝒒}}.`$ (47b)
Comparing these correlators with Eqs. (8) and (9), taking into account Eqs. (45), one finds

$`F_\xi `$ $`=`$ $`2T\chi _{I5}(D+\kappa _2\kappa _1),`$ (48a)
$`F_\eta `$ $`=`$ $`{\displaystyle \frac{2T\kappa _1}{f^2}}.`$ (48b)
Our results for the hydrodynamic correlations functions are given by Eqs. (44) and (48). For completeness, we give here also the result for the cross correlator of $`\phi ^a`$ and $`A_0^a`$:

$$d^4xe^{iqx}\phi ^a(x)A_0^b(0)=\frac{2iTq_0\mathrm{\Gamma }_𝒒}{[(q_0\omega _𝒒)^2+\frac{1}{4}\mathrm{\Gamma }_𝒒^2][(q_0+\omega _𝒒)^2+\frac{1}{4}\mathrm{\Gamma }_𝒒^2]}.$$
(49)
It is instructive to compare the correlator of the axial isospin charge density $`A_0^a`$ \[Eq. (44a)\] with that of the vector isospin charge density $`V_0^a`$. The dynamics of the latter is completely diffusive and is given by the equations

$`_0V_0^aD_I^2V_0^a`$ $`=`$ $`_i\zeta ^{ai},`$ (50a)
$`\zeta ^{ai}(x)\zeta ^{bj}(0)`$ $`=`$ $`2TD_I\chi _I\delta ^{ab}\delta ^{ij}\delta ^4(x),`$ (50b)
where $`D_I`$ is the diffusion constant for the isospin charge, and $`\chi _I`$ is the (vector) isospin susceptibility. The correlator of $`V_0^a`$ is

$$d^4xe^{iqx}V_0^a(x)V_0^b(0)=\frac{2TD_I\chi _I𝒒^2}{q_0^2+D_I^2𝒒^4}.$$
(51)
The pole of the correlator is located at purely imaginary frequencies $`q_0=\pm iD_I𝒒^2`$, as it should be for a purely diffusive mode.

Using Eqs. (51) and (47) we can understand the limitations of the Lagrangian approach and the role of the dissipative processes. For example, one could try to apply the method of Sec. III to determine the vector isospin susceptibility $`\chi _I`$, introducing an isospin chemical potential $`\mu _I`$. At zero temperature such a method was used in Ref. mui . As expected, the isospin susceptibility $`\chi _I`$ vanishes at $`T=0`$ since it takes finite energy to excite isospin degrees of freedom (pions) and change the isospin density. However, at nonzero temperature, the naive application of the effective Lagrangian method would predict that $`\chi _I=0`$ also.555We thank T. Schäfer for pointing this out to us. On the other hand, we should expect $`\chi _I0`$ at finite temperature, even at very small $`T`$, because of the presence of isospin-carrying pions in the thermal medium. Looking at the correlator (51), one sees that the equal-time correlator defining $`\chi _I`$, being the integral of Eq. (51) over $`q_0`$ in the limit $`𝒒0`$, receives the main contribution from very slow (diffusive) modes: $`q_0D_I𝒒^2`$. However, the effective Lagrangian (14) only describes faster (propagating) modes of the $`\phi `$ field with $`q_0=\omega _𝒒D_I𝒒^2`$. The slow diffusive modes which contribute to $`\chi _I`$ are not present in the Lagrangian (14). For the axial isospin susceptibility $`\chi _{I5}`$, the situation is different, for the integral of Eq. (51) over frequencies in the limit $`𝒒0`$ is concentrated entirely near values $`q_0=\pm \omega _𝒒`$. Thus, unlike $`\chi _I`$, the value of $`\chi _{I5}`$ at finite $`T`$ can be found correctly using the effective Lagrangian approach.

### IV.5 Comparison to previous results

As a by-product of our analysis, we obtain a generalization of the Gell-Mann–Oakes–Renner relation to finite temperature:

$$f^2m^2=m_q\overline{\psi }\psi .$$
(52)
This equation has the same form as at zero temperature. At finite temperature, $`m`$ should be understood as the pion screening mass, and the exact meaning of the temperature-dependent “pion decay constant” $`f^2`$ is given by Eq. (43). One can also write the GOR relation in an alternative form,

$$\chi _{I5}m_p^2=m_q\overline{\psi }\psi ,$$
(53)
where $`m_p=um`$ is the “pole mass” of pions.

That the velocity of pions at finite temperature is different from the velocity of light has been seen in the second order of chiral perturbation theory (i.e., in order $`T^4/f_\pi ^4`$) Pisarski ; Toublan . The authors of Ref. Pisarski  also introduced two pion decay constants $`f^t`$ and $`f^s`$, which correspond to our $`\sqrt{\chi _{I5}}`$ and $`f`$, and checked the validity of the GOR relation (53). However, these constants were defined only in chiral perturbation theory, and only at small $`T`$; there has not been any attempt to give a precise definition of the constants $`f^t`$ and $`f^s`$ at temperatures comparable to $`T_c`$. We have, in contrast, given a precise meaning to the constants $`\chi _{I5}`$ and $`f`$ in terms of equal-time correlation functions that can be measured on the lattice. The GOR relation now contains only well-defined quantities. Although Eq. (53) cannot be checked on the lattice, since it contains the pion pole mass, the version (52) can be verified numerically since all quantities entering it are statically measurable.

We also note that Eqs. (30), without the noise and dissipation terms, can be obtained by linearizing the hydrodynamic equations obtained in Ref. Sonhydro  by the Poisson-bracket technique. One can then identify the parameters $`f_t`$ and $`f_s`$ in Ref. Sonhydro  as $`f_t^2=\chi _{I5}`$, $`f_s^2=f^2`$ (the parameter $`f_v^2`$ of Ref. Sonhydro ) is equal to the vector isospin susceptibility). We did not try to reproduce the full nonlinear hydrodynamic equation of Ref. Sonhydro  in the present approach.

## V An explicit example: a linear sigma model

It is instructive to explicitly verify our relation (2) between the velocity of the Goldstone bosons, the axial isospin susceptibility, and the temperature-dependent decay constant in a model where weak-coupling calculations are possible. The simplest theory with spontaneous breaking of a continuous symmetry is the linear sigma model. We shall show that Eq. (2) indeed holds in this model to leading order of perturbation theory

It is important to emphasize that the sigma model considered in this section is only meant to serve as an example of a theory where we can explicitly check our results by perturbative calculations. These models do not describe QCD at $`T`$ of order $`T_c`$ — a theory with no apparent small parameter. We claim, however, that the connection between the dispersion relation of the pions and static correlation functions is model independent, and can be derived from only a few general assumptions stated in Sec. IV.

We start with the following Lagrangian

$$=\frac{1}{2}^\mu \varphi ^a_\mu \varphi ^a+\frac{\mu ^2}{2}\varphi ^a\varphi ^a\frac{\lambda }{4}(\varphi ^a\varphi ^a)^2,$$
(54)
where $`a=1,2,\mathrm{},N`$. We assume $`\lambda 1`$, so perturbation theory can be applied. When $`\mu ^2>0`$, the $`O(N)`$ symmetry is spontaneously broken at zero temperature. We choose the vacuum to align in the $`\varphi ^N`$ direction. At zero temperature, and at tree level, the vacuum expectation value of $`\varphi ^N`$ is

$$\varphi ^Nv_0=\frac{\mu ^2}{\lambda }.$$
(55)
At finite temperature the expectation value of $`\varphi ^N`$ is different from $`v_0`$. Denoting it as $`v`$ and replacing $`\varphi ^N=v+\sigma `$ in the Lagrangian (54), we obtain

$$\begin{array}{cc}\hfill & =\frac{1}{2}(_\mu \sigma )^2+\frac{1}{2}(_\mu 𝝅)^2(\lambda v^2\mu ^2)v\sigma \frac{1}{2}(3\lambda v^2\mu ^2)\sigma ^2\frac{1}{2}(\lambda v^2\mu ^2)𝝅^2\hfill \\ & \lambda v\sigma ^3\lambda v\sigma 𝝅^2\frac{\lambda }{4}\sigma ^4\frac{\lambda }{4}(𝝅^2)^2\frac{\lambda }{2}\sigma ^2𝝅^2,\hfill \end{array}$$
(56)
where $`𝝅=(\varphi ^1,\varphi ^2,\mathrm{},\varphi ^{N1})`$. From Eq. (56) the Feynman rules can be easily written down.

The value of $`v`$ is determined by the condition that $`\sigma =0`$. For computation, we will use the Matsubara (Euclidean) formalism, in which this condition reads, to one-loop order,

$$(\lambda v^2\mu ^2)v+\lambda vT\underset{p_0}{}\frac{d𝒑}{(2\pi )^3}\left(\frac{3}{P^2+m_\sigma ^2}+\frac{N1}{P^2}\right)=0,$$
(57)
In this Section, $`P^2=p_0^2+p^2`$, $`p=|𝒑|`$. The sum-integral in Eq. (57) is ultraviolet divergent. This temperature-independent divergence can be absorbed into the redefinition of $`\mu ^2`$. The temperature dependence of $`v`$ comes from the thermal part of the Lagrangian. We shall be interested in temperatures of order $`v_0`$, so $`Tm_\sigma `$. In this case,

$$\lambda v^2\mu ^2+\frac{N+2}{12}\lambda T^2=0,$$
(58)
or

$$v^2(T)=v_0^2\left(1\frac{T^2}{T_c^2}\right)\text{ with }T_c^2=\frac{12}{N+2}\frac{\mu ^2}{\lambda }.$$
(59)
Equation (59) is valid everywhere except for a narrow Ginzburg region near $`T_c`$.

The masses of $`\sigma `$ and $`\pi `$ are computed, e.g., in Ref. Kapusta . At the one-loop level $`m_\sigma `$ receives contributions from four diagrams

(where $`\sigma `$ propagators are denoted by solid lines, and $`\pi `$ propagators are drawn as dashed lines). The last two bubble diagrams are negligible in the finite-temperature regime we consider, so

$$m_\sigma ^2=3\lambda v^2\mu ^2+\lambda vT\underset{p_0}{}\frac{d𝒑}{(2\pi )^3}\left(\frac{3}{P^2+m_\sigma ^2}+\frac{N1}{P^2}\right)=2\lambda v^2,$$
(60)
and thus decreases with temperature. The Goldstone boson receives corrections to its mass from three one-loop graphs,

and remains massless since

$$\begin{array}{cc}\hfill m_\pi ^2& =\lambda v^2\mu ^2\mu ^2+\lambda vT\underset{p_0}{}\frac{d𝒑}{(2\pi )^3}\left(\frac{1}{P^2+m_\sigma ^2}+\frac{N1}{P^2}\right)\hfill \\ & 4\lambda ^2v^2T\underset{p_0}{}\frac{d𝒑}{(2\pi )^3}\frac{1}{P^2(P^2+m_\sigma ^2)}=0.\hfill \end{array}$$
(61)
To find the dispersion relation of the pion near $`q=0`$, one has to expand the pion self-energy $`\mathrm{\Sigma }(q)`$ in powers of $`q_0`$ and $`𝒒`$. The only one-loop diagram that depends on $`q`$ is the bubble diagram, which can be evaluated either by using the Schwinger-Keldysh real-time formalism, or by doing calculations in the Matsubara formalism and then performing an analytic continuation from imaginary to real $`q_0`$. By either method one finds

$`\mathrm{\Sigma }(q)`$ $`=`$ $`4\lambda ^2v^2{\displaystyle }{\displaystyle \frac{d𝒑}{(2\pi )^32\omega _{𝒑+𝒒}2\mathrm{\Omega }_𝒑}}\{[f(\omega _{𝒑+𝒒})+f(\mathrm{\Omega }_𝒑)+1]({\displaystyle \frac{1}{\mathrm{\Omega }_𝒑+\omega _{𝒑+𝒒}q_0}}+{\displaystyle \frac{1}{\mathrm{\Omega }_𝒑+\omega _{𝒑+𝒒}+q_0}})`$ (62)
$`+[f(\omega _{𝒑+𝒒})f(\mathrm{\Omega }_𝒑)]({\displaystyle \frac{1}{\mathrm{\Omega }_𝒑\omega _{𝒑+𝒒}q_0}}+{\displaystyle \frac{1}{\mathrm{\Omega }_𝒑\omega _{𝒑+𝒒}+q_0}})\},`$
where $`\omega _𝒑`$ and $`\mathrm{\Omega }_𝒑`$ are the energies of the $`\pi `$ and $`\sigma `$ particles with momentum $`𝒑`$: $`\omega _𝒑=p`$, $`\mathrm{\Omega }_𝒑=(p^2+m_\sigma ^2)^{1/2}`$; and $`f(\omega )`$ is the Bose-Einstein distribution function $`f(\omega )=(e^{\beta \omega }1)^1`$.

We need to compute the coefficients of $`q_0^2`$ and $`q^2`$ in the expansion of $`\mathrm{\Sigma }(Q)`$. To compute the coefficient of $`q^2`$, one can set $`q_0=0`$ and Eq. (62) becomes

$$\mathrm{\Sigma }(0,𝒒)=4\lambda ^2v^2\frac{d𝒑}{(2\pi )^3}\frac{1}{\mathrm{\Omega }_𝒑^2\omega _{𝒑+𝒒}^2}\left[\frac{f(\mathrm{\Omega }_𝒑)+\frac{1}{2}}{\mathrm{\Omega }_𝒑}\frac{f(\omega _{𝒑+𝒒})+\frac{1}{2}}{\omega _{𝒑+𝒒}}\right].$$
(63)
The right hand side of Eq. (63) can be expanded in powers of $`q`$. The constant is compensated by the other diagrams. The $`q^2`$ term gives rise to an integral which is dominated in the infrared region $`pm_\sigma `$, and is of order

$$\mathrm{\Sigma }(0,𝒒)\frac{\lambda ^2v^2T}{m_\sigma ^3}𝒒^2=O(\lambda ^{1/2})𝒒^2+\mathrm{}.$$
(64)
In the last equation we assumed $`Tv_0`$, and $`m_\sigma ^2\lambda T`$, which is valid when $`T`$ is not very close to $`T_c`$.

Now let us put $`𝒒=0`$ and expand in $`q_0`$. Since the constant term is canceled out by other diagrams, one has to look only at terms of order $`q_0^2`$. One finds

$$\mathrm{\Sigma }(q_0,\mathrm{𝟎})=8\lambda ^2v^2q_0^2\frac{d𝒑}{(2\pi )^32\omega _𝒑2\mathrm{\Omega }_𝒑}\left[\frac{f(\omega _𝒑)+f(\mathrm{\Omega }_𝒑)+1}{(\mathrm{\Omega }_𝒑+\omega _𝒑)^3}+\frac{f(\omega _𝒑)f(\mathrm{\Omega }_𝒑)}{(\mathrm{\Omega }_𝒑\omega _𝒑)^3}\right].$$
(65)
The first term in the square brackets gives rise to an integral that is dominated by the infrared, i.e., by $`pm_\sigma `$, and so is completely analogous to the coefficient of $`p^2`$ above. The coefficient of $`q_0^2`$ coming from the first term is hence of order $`O(\lambda ^{1/2})`$. In contrast, the integral of the second term in the square brackets is dominated by $`pT`$. For such $`p`$ one can write, approximately,

$$\mathrm{\Omega }_𝒑\omega _𝒑=\frac{m_\sigma ^2}{2p},f(\omega _𝒑)f(\mathrm{\Omega }_𝒑)=\frac{m_\sigma ^2}{2p}\frac{n(p)}{p},$$
(66)
which gives

$$\mathrm{\Sigma }(q_0,\mathrm{𝟎})=\frac{8\lambda ^2v^2}{m_\sigma ^4}q_0^2\frac{d𝒑}{(2\pi )^3}\frac{f(p)}{p}=\frac{4\lambda ^2v^2}{3m_\sigma ^4}T^2q_0^2.$$
(67)
Substituting $`m_\sigma ^2=2\lambda v^2`$, we finally find

$$\mathrm{\Sigma }(q_0,𝒒)=\frac{T^2}{3v^2}q_0^2.$$
(68)
Since $`T`$ and $`v`$ are both of order $`v_0`$, the coefficient in front of $`q_0^2`$ is of order $`𝒪(\lambda ^0)`$, in contrast to that of $`𝒒^2`$. The velocity of pions in our model is

$$u^2=\left(1+\frac{T^2}{3v^2}\right)^1,Tm_\sigma .$$
(69)
It is instructive to compare this result to the one obtained in the framework of the chiral perturbation theory: $`u=1𝒪(T^4/f_\pi ^4)`$ Pisarski ; Toublan . One should bear in mind that the regime in which the chiral perturbation theory result applies corresponds to $`Tm_\sigma `$. In our linear sigma model we consider a different regime: $`Tvm_\sigma \sqrt{\lambda }v`$. In the chiral perturbation theory $`T/f_\pi `$ serves as an expansion parameter. In our weak-coupling calculation the expansion is in $`\lambda `$ while the $`T`$ dependence is included to all orders in $`T/v`$.

Now we need to check that this coincides with $`f^2/\chi _{I5}`$. Recall that our definition of $`f^2`$ is as follows: if we define $`\phi ^a=\pi ^a/v`$, then

$$𝑑𝒙e^{i𝒒𝒙}\phi ^a(t,𝒙)\phi ^b(t,\mathrm{𝟎})=\frac{T}{f^2}\frac{\delta ^{ab}}{𝒒^2}.$$
(70)
The left hand side is proportional to the $`\pi `$ Matsubara propagator, summed over $`q_0`$,

$$\frac{T}{v^2}\underset{q_0}{}\frac{1}{q_0^2+𝒒^2}.$$
(71)
When $`|𝒒|T`$ the dominant term in the sum is the one with $`q_0=0`$, therefore,

$$f^2=v^2.$$
(72)
Taking into account one-loop graphs, as we have seen, will change $`f^2`$ only by an amount of order $`O(\lambda ^{1/2})`$.

Now to compute the susceptibility $`\chi _{I5}`$ we need to turn on a chemical potential coupled to a broken charge. There are $`N1`$ broken generators. Let us consider the one that transforms $`\sigma `$ and $`\pi _1`$ into each other. The change of the Lagrangian when the corresponding chemical potential is turned on is

$$\delta =\mu (\pi _1_0\sigma \sigma _0\pi _1)+\frac{\mu ^2}{2}[(v+\sigma )^2+\pi _1^2].$$
(73)
The susceptibility can be computed in Matsubara formalism. There are two contributions: one from the $`\mu ^2`$ term in $`\delta `$, the other from the bubble graph:

$$\begin{array}{cc}\hfill \chi _{I5}& =v^2+2T\underset{p_0}{}\frac{d𝒑}{(2\pi )^3}\frac{1}{p_0^2+p^2}4T\underset{p_0}{}\frac{d𝒑}{(2\pi )^3}\frac{p_0^2}{(p_0^2+p^2)^2}\hfill \\ & =v^2+\frac{1}{2T}\frac{d𝒑}{(2\pi )^3}\mathrm{sinh}^2\frac{p}{2T}=v^2+\frac{T^2}{3}.\hfill \end{array}$$
(74)
The interpretation of this formula is rather direct: $`v^2`$ is the contribution of the condensate, and $`T^2/3`$ is the contribution from the free gas of $`\sigma `$ and $`\pi _1`$. The square of the velocity of pions, Eq. (69), is equal to the ratio of $`f^2`$ in Eq. (72) and $`\chi _{I5}`$ in Eq. (74), which is what we need to verify. We did not, however, attempted to turn on an explicit symmetry breaking and verify, e.g., the GOR relation at finite temperature. Such a calculation should be straightforward.

## VI Conclusion

Our goal has been to demonstrate that the dispersion relation of pions can be expressed in terms of quantities obtainable from equal-time or static correlation functions. The precise relation is given in Sec. II. Our result enables one to find the real part of the pion dispersion relation on the lattice. However, it does not enable one to compute the imaginary part, which characterizes the damping of pion modes.

Nowhere in our treatment did we assume any condition on the temperature (except $`T<T_c`$), as long as we stay infinitesimally close to the chiral limit. Thus, if $`m_q`$ is infinitesimally small, our result applies to all temperatures smaller than the temperature of the chiral phase transition $`T_c`$. However, at any fixed (small) value of $`m_q`$, our results do not apply at temperatures too close to critical, where the pion screening mass becomes of the same order as the screening mass of the sigma meson. The width of this region near $`T_c`$ shrinks to zero as a power of $`m_q`$ (more precisely, as $`m_q^{1/\beta \delta }`$). Our treatment must fail there because, as explained in Sec. IV, the sigma boson also needs to be included into the hydrodynamic theory. However, the scaling of different quantities in this temperature region can be determined using scaling and universality arguments, as discussed in Ref. critical\_pions . In the Appendix of the present paper we derive some additional interesting scaling properties omitted in Ref. critical\_pions .

Finally, assuming $`m_q`$ to be very small, as $`TT_c`$, the chiral condensate $`\overline{\psi }\psi 0`$. As shown in Ref. critical\_pions , this implies that also $`f0`$ (although with a slightly different critical exponent). On the other hand, the axial isospin susceptibility $`\chi _{I5}`$ becomes degenerate with the (vector) isospin susceptibility at $`T_c`$, where both remains finite. One concludes from Eq. (2) that the pion velocity $`u`$ tends to zero as one approaches the critical temperature. In fact, it can be shown that $`u`$ approaches zero faster than the divergence of the screening mass $`m`$, so the pion pole mass $`m_p=um`$ goes to zero as $`TT_c`$ critical\_pions .

###### Acknowledgements.

We thank Daniel Boyanovsky, Sangyong Jeon, Guy Moore, Thomas Schäfer, Edward Shuryak and Larry Yaffe for stimulating discussions. The authors thank RIKEN, Brookhaven National Laboratory, and the U.S. Department of Energy (DE-AC02-98CH10886) for providing the facilities essential for the completion of this work. D.T.S. is supported, in part, by DOE Grant No. DOE-ER-41132 and by the Alfred P. Sloan Foundation. M.A.S. is supported by a DOE OJI grant.

## Appendix A Dynamical scaling at the chiral transition

In this appendix we provide a simple derivation of the dynamical critical exponent $`z`$ characterizing critical slowing down near the QCD chiral phase transition. It is based on a scaling argument similar to the one used in Ref. Hohenberg  to find $`z`$ in antiferromagnets.

For simplicity, here we consider the chiral limit $`m_q=0`$. Scaling and universality arguments presented in Ref. critical\_pions  predict that the velocity $`u`$ vanishes as $`T`$ approaches $`T_c`$ from below. The quantity $`u^2`$ scales as $`u^2t^{(d2)\nu }`$, where $`t(T_cT)/T_c`$. Since the inverse correlation length of the order parameter $`\overline{\psi }\psi `$, i.e., the static screening mass of the sigma particle $`m_\sigma `$, scales as $`m_\sigma t^\nu `$, we conclude that for $`t1`$

$$u^2m_\sigma ^{d2}.$$
(75)
The fact that $`u0`$ at $`T_c`$ means that the dispersion relation ceases to be linear. Moreover, the effect of damping also becomes important; in other words, we expect the real and imaginary parts of $`\omega `$ to become comparable. In such a situation on cannot refer to $`\omega `$ as a quasiparticle energy. Rather, it is a characteristic frequency, or inverse relaxation time, of a mode of a given wave number $`𝒑`$. Scaling hypothesis dictates that the relation between $`𝒑`$ and $`\omega `$ should be homogeneous, i.e., $`\omega |𝒑|^z`$. The dynamical scaling exponent $`z`$ is, in a generic case, new exponent independent of the static exponents, e.g., $`\nu `$ and $`\eta `$. However, in the case of QCD (similar to the case of an antiferromagnet RajagopalWilczek ; Hohenberg ), since the dispersion relation of pions is given in terms of static quantities only, it turns out that $`z`$, as one would expect, can be derived from static scaling laws only.

To determine $`z`$ we observe that the dispersion relation $`\omega |𝒑|^z`$ applies at scales $`m_\sigma |𝒑|T`$. At softer scales, $`|𝒑|m_\sigma `$, the dispersion relation is still linear, $`\omega =u|𝒑|`$. Requiring that the two expressions for $`\omega `$ match at $`|𝒑|m_\sigma `$ we find

$$m_\sigma ^zum_\sigma m_\sigma ^{d/2},$$
(76)
where to obtain the last scaling relation we used Eq. (75). Therefore, at the QCD chiral phase transition the dynamical critical exponent is given by

$$z=\frac{d}{2}$$
(77)
(see also an alternative derivation in Refs. RajagopalWilczek ; Hohenberg ).666 This result agrees with RajagopalWilczek ; Hohenberg , but disagrees with Boyanovsky .

A similar argument can be applied to determine the scaling of the diffusion coefficient $`D^{}=D+\kappa _2`$, characterizing the pion damping in the chiral limit \[Eq. (45b)\]. The small momentum expansion of $`\omega `$ is $`u|𝒑|\frac{i}{2}D^{}𝒑^2+\mathrm{}`$ . The scaling hypothesis dictates that all terms in this expansion become of the same order in magnitude when $`|𝒑|m_\sigma `$ (otherwise, another scale, e.g., $`u/D^{}`$, appears in addition to the inverse correlation length $`m_\sigma `$). This requires

$$D^{}m_\sigma ^{(d4)/2}t^{\nu (d4)/2}.$$
(78)
Thus $`D^{}`$ diverges near the chiral phase transition in QCD ($`d=3`$).