Dynamic Universality Class of Model H with Frustrated Diffusion: $\epsilon$-Expansion

We study a variation of the dynamic universality class of model H in a spatial dimension of $d=4-\epsilon$, by frustrating charge diffusion and momentum density fluctuations along $d_T=1$ or $d_T=2$ dimensions, while keeping the same dynamics of model H in the other $d_L=d-d_T$ dimensions. The case of $d_T=2$ describes the QCD critical point in a background magnetic field. We find that these models belong to a different dynamical universality class due to extended conservation laws compared to the model H, although the static universality class remains the same as the 3-dimensional Ising model. We compute the dynamic critical exponents of these models in first order of $\epsilon$-expansion to find that $x_\lambda\approx 0.847\,\epsilon$, $x_{\bar\eta}\approx 0.153\,\epsilon$, and $z=4-x_\lambda\approx 3.15$ when $\epsilon=1$ and $d_T=2$. For $d_T=1$ the results are numerically similar to the model H values: $z\approx 3.08$.


Introduction
In this work we study a particular variation of dynamic universality class of model H [1,2] in the -expansion [3,4]. Recall that the model H has two dynamical fields: a conserved charge density and the shear component of momentum density fluctuations. Their dynamics in long wave-lengths near the critical regime is mainly governed by diffusion dynamics due to their conservation laws. What gives a non-trivial IR fixed point is a non-linear coupling between these two diffusion dynamics [5,6]. This coupling constant at the fixed point is of order , and a perturbation theory in terms of -expansion is possible. The results for dynamical critical exponents such as x λ , x η and z are available up to 2 [1]. When = 1, they are x λ ≈ 1, x η ≈ 0 and z ≈ 3.
The model we study is a simple twist of the model H: we will frustrate the charge conductivity and the momentum density fluctuations along d T = 1 or d T = 2 spatial dimensions (denoted as x T ), while keeping them in the rest d L = 4 − d T − dimensions (denoted as x ). There are at least two motivations to consider these: 1) d T = 1 case can describe the layered system where the charge transport and the momentum flow across the layers are frustrated or damped with a finite relaxation time, so that these modes decouple near the critical regime, 2) d T = 2 case describes a system in a background magnetic field F 12 = 0, where momentum fluctuations in (1,2) directions are damped with a finite relaxation time and are no longer a hydrodynamic variable [7]. In strong magnetic field limit, a charge conductivity along (1,2) directions should be suppressed as well due to small cyclotron orbits. The d T = 2 case should describe the dynamic universality class for the QCD critical point [8,9,10] in a strong background magnetic field.
Let us explain the case 2) in some more detail [7]. The hydrodynamic equations with a background magnetic field are where J µ is the charge current, T µν = wu µ u ν + pη µν is the energy-momentum, u µ is the fluid velocity field and λ is the conductivity. For small transverse velocity v 1 , v 2 , the first equation gives the current J i = λF 12 ij v j (i, j = 1, 2). In the homogeneous limit (k → 0), the second equation gives so that v 1,2 has a finite relaxation time τ R in k → 0 limit, that is, it is not a hydrodynamic mode. Once we remove v 1,2 from the critical dynamical modes, there is no non-linear coupling that gives renormalization for the transverse conductivity along (1, 2) directions: the physical transverse conductivity remains finite in infrared. Note that the conductivity appearing in (1.2) is the transverse conductivity. The naive scaling dimension of the conductivity is negative when z < 4 (see (3.12)) and the transverse conductivity therefore becomes irrelevant in the infrared in the renormalization group. On the other hand, the physical longitudinal conductivity will be shown to diverge near the critical point due to the non-linear coupling, and its renormalized value goes to a finite fixed point in the renormalization group (see the discussion in section 4).
We will find that these variations of model H belong to dynamic universality classes that are different from the original model H, and the dynamic critical exponents such as x λ , x η and z are different: we compute these exponents in first order of the -expansion by using the Wilsonian renormalization group method [3,4]. Using the Feynman diagram method of Ref. [11] one can in principle continue to higher orders in . The reason why we get the different dynamic universality classes in these models is that the symmetries are enlarged from those of the original model H. In the absence of charge diffusion along d T , the integrated charge density along d L at any given point in d T is conserved: the symmetry group is infinite dimensional.
Despite that the real-time dynamics in these models breaks rotational symmetry, we assume that the static universality class remains to be the isotropic 3D Ising model described by a theory of a scalar field ψ with quartic interaction. This is reasonable since the only relevant or marginal perturbation near critical regime that breaks rotational symmetry is of a form K ij ∂ i ψ∂ j ψ and we can make it isotropic by diagonalizing and rescaling the coordinates.

Description of the model
We will follow the notations in Ref. [1,2], and denote a conserved order parameter by ψ and the momentum density vector by j which has components only along the space x of dimension d L = 4 − − d T in our model. The static equilibrium thermal distribution is given by e −F [ψ,j ] where F is the free energy functional of 3D Ising universality class where Λ is the physical UV-cutoff in momentum space. Writing the parameters of the model in the above way, (r, u) are dimensionless. The dynamical model is a simple twist of the model H with frustrated diffusion along x T , where P projects the vector components onto the subspace perpendicular to a momentum vector k in Fourier space, that is, it keeps only the shear components. Since j is already perpendicular to k ⊥ , it is practically a projection operator in x space; The random noises (θ, ξ) are Gaussian with the strength determined by Fluctuation-Dissipation relation, which ensures that the equilibrium distribution is e −F . The λ is the conductivity,η ,⊥ are the shear viscosities in x and x ⊥ spaces respectively, which are in general different.
We will see that their values at the IR fixed point are indeed different, but comparable to each other so we need to keep them both. The g is the non-linear coupling between two dynamical modes in model H, which drives a non-trivial IR fixed point. The exponent z is the scale dimension of the frequency ω relative to the momentum k. All parameters defined as above are then dimensionless.
Although we can perform the Wilsonian renormalization group at the level of the above equations of motion, we choose to work in a language of stochastic field theory or a path integral, where the renormalization procedure looks more organized (at least to the eyes of the author). For that purpose, we introduce the "a-type" fields (ψ a , j a ) and call the original variables in the equations of motion the "r-type" fields (ψ r , j r ), then consider a path integral of e S with an action S = dt d 4− x L, The correspondence to the usual (r, a)-fields in the Schwinger-Keldysh formalism is clear: the "r-type" fields are classical variables. The path integral over the "a-type" field localizes the path integral of "r-type" fields to the solutions of the original stochastic Langevin equations of motion with the Gaussian random noises (θ, ξ). The "wave-function" at time t is precisely then the probability distribution of the classical variables at time t generated by solutions of the stochastic equations of motion.
We can integrate out the noise variables (θ, ξ) as they are Gaussian to obtain Upon "quantization" the "a-type" fields are canonical conjugate to the "r-type" fields:

Renormalization group in -expansion
We follow the standard procedure (see Ref. [1,3,4]) of thinning the momentum shell around the cutoff Λ, and integrate over the shell with a constant b > 1 close to 1. After this we rescale the coordinates or equivalently the momenta to get back to the same cutoff Λ in the rescaled momentum space k ; Without the non-linear couplings such as g and u, the parameters of the theory would change at each step simply by their naive scaling dimensions: There are similar equations for the static parameters (r, u). The non-linear couplings by (u, g) give rise to additional contributions to the above. These are what we need to compute.
In the critical regime, the IR cutoff set by the correlation length Λ IR ∼ ξ −1 is far separated from the UV cutoff Λ. After N steps of the above procedure where b N ξ −1 = Λ, the UV and IR cutoffs in the renormalized theory become comparable and the scaling behavior is lost. This is the point where the hydrodynamic regime sets in, and the renormalization group running of the parameters of the theory stops and further contributions to these parameters are IR-finite. The number of steps of the above procedure to be performed to reach this is N = log(Λξ)/ log b. The renormalized parameter r starts close to its fixed point value which is of order r * ∼ , and we can neglect it in the propagator in the leading order perturbation in . In each step of the above procedure, the deviation of the renormalized r from r * grows, and only after performing the same N = log(Λξ)/ log b steps it becomes of order 1: this is because the scaling behavior is lost precisely when the renormalized r deviates from r * ∼ by order 1. Therefore r stays of order in most of the N -steps, and we can neglect r 1 in the propagators in all N -steps in the leading perturbation in -expansion [3,4].
When z > 1, the cutoff in frequency space can be taken to be infinite. Even if we started with a same cutoff Λ in the frequency space, a step of above procedure would change it to b z−1 Λ. In the critical regime where N 1, this cutoff quickly becomes much larger than Λ in most of the N steps.
The observed physical parameters of the theory are the ones without rescaling the coordinates and the time: they are the parameters measured in terms of the original coordinates and time. The effect of rescaling coordinates for them is simply given by (3.12), and therefore the observed physical parameters are obtained from the renormalized parameters by undoing the naive scaling transformation (3.12) [1,3,4]: where the starred parameters are the renormalized parameters after performing the large N -steps (they could be a finite fixed point value or not). In effect, the physical parameters in the above capture only the contributions from the non-linear couplings in the relevant momentum range ξ −1 < |k| < Λ as they should.
With these general discussions reviewed, we now compute the contributions from the non-linear couplings (g, u) to the renormalization group equation. We can do a perturbation theory in these couplings since their fixed point values will be of order . After integrating over the momentum shell, the action S + δS with a new cutoff Λ/b is expected to be given by a new set of parameters (3.14) The contribution to g will be shown to be absent. The contributions from the non-linear couplings (δλ , δη ⊥, , δF ) are of first order in log b in small log b limit. The effect of rescaling (3.12) to restore the original cutoff Λ and the above non-linear contribution are therefore additive to each other to first order in log b.
First, the δF should be identical to what one would have in a static renormalization group [1], because the equilibrium distribution after integrating over the momentum shell is e −(F +δF ) by Fluctuation-Dissipation relation, and this must agree with what one would have in a static renormalization after integrating over the same momentum shell.
Although this is guaranteed by the Fluctuation-Dissipation relation in the original unintegrated theory, it is assuring to see this explicitly for a few lowest order contributions, which we have checked. From the known static renormalization group [3,4], we have where η/2 ∼ 2 is the anomalous dimension of ψ, and the other terms are not of our interest for the renormalization of the transport coefficients (r + δr for example is order in the scaling regime and will be neglected in the propagator). Looking at how F appears in the charge diffusion term in the action (2.9), this wave-function renormalization contributes to δλ as The other leading contributions to (δλ , δη ⊥, ) come from the coupling g. One can compute these by looking at the contributions to the ψ a ψ r and j a j r terms in the action as they contain these transport coefficients. The real-time Feynman diagrams generating a ψ a ψ r term in δS are shown in Figure 1. There are two real-time diagrams to be summed.
In Figure 2 we show the real-time diagram generating a j a j r term in δS. The solid lines are the propagators of ψ-field and the wavy lines are for the j -fields. These propagators are easily found from the quadratic part of the action S to be where k = (ω, k) = (ω, k , k ⊥ ) is a frequency-momentum in Fourier space. Each vertex in the diagrams is from the g coupling in the action, which is We first write down the expression for the first diagram in Figure 1.
where p = (Ω, p) is the external momentum, and we denote Computing ω integration by closing the contour in the lower half-plane, we obtain where k integration to be performed in the momentum shell Λ/b < |k| < Λ. The second diagram in Figure 1 is similarly computed to be 3.22) and the sum of the two becomes where we take a limit of a small external momentum compared to the loop momentum: Here comes an important step of approximation in the model H which is self-consistent.
We will see thatη ⊥, grows large in the renormalization group running, while λ approaches a finite fixed point in the IR: after many steps of renormalization we haveη ⊥, λ . Since the loop momentum k in the shell is near the cutoff Λ, the term with λ in the denominator is then negligible compared to the terms withη ⊥, . This gives finally the expression (3.24) We now describe the k integral in the above. We have with an angular function on the sphere F (Ω). The |k| integral gives (3.27) and in the leading order in -expansion, the angular integral can be performed in = 0 limit, that is, on the S 3 sphere. We parameterize the k space in this limit as k 1 = |k| cos θ , k 2 = |k| sin θ cos φ , k 3 = |k| sin θ sin φ cos χ , with 0 < (θ, φ) < π and 0 < χ < 2π, and the measure is dΩ = sin 2 θ sin φdθdφdχ . The result should not depend on these choices as one can check. For d T = 2, we choose k ⊥ = (k 1 , k 2 ) and k = (k 3 , k 4 ), and p pointing to the k 4 direction. Then F (Ω) becomes With these the (3.24) becomes at leading order in +i where the last line is a space-time expression. Looking at the diffusion term in the action (2.9), this corresponds to a non-linear contribution we are looking for (3.33) We next compute the diagram in Figure 2 for the renormalization of shear viscosities.
We have (neglecting an external frequency) and performing the ω-integration, we obtain Recall that j (p) lies in k space, and is also perpendicular to p (and hence to p) inside the k space to be a shear component. Therefore we can replace in the numerator where k ⊥ is the subspace inside k perpendicular to p , and its dimension is Also expanding the remaining part of the integrand up to the first relevant leading order, the generated action (3.35) then becomes which is logarithmically sensitive to the cutoff. The k integral is done similarly as before.
There are two types of terms in the result: a term proportional to p 2 ⊥ and the other proportional to p 2 . They correspond to a renormalization ofη ⊥ andη respectively. The first type of term is and the second type is (3.40) What remains in these terms is an angular integration on S 3 in = 0 limit. Using the parametrization (3.28), we obtain for d T = 1 case From these, we finally obtain the non-linear contributions to the renormalization ofη ⊥ andη , The non-linear contributions to δg can potentially arise from the diagrams in Figure 3.
Explicit computations show that these diagrams generate only higher dimensional terms, and a correction to δg is absent [1].
In summary of all these, the renormalization group equations are where the angular function F (Ω,η ⊥ ,η ) is given in (3.30) and (3.31).

Fixed point and the critical exponents
Following the model H analysis [1,2], we define the two coupling constants, Using these variables, the flow equations (3.45) can be written as where the primed parameters are the renormalized ones after performing one step of integrating a momentum shell and rescaling, and the front factors of b are from naive dimensional scaling, which should be undone for the physical parameters measured in the original coordinate-time. Therefore the physical parameters receive contributions only from the non-linear effects. We see that these effects are given solely by the parameters f ⊥ and f . More explicitly, for infinitesimal log b, where λ 0 ,η 0 ⊥, are the parameters at the cutoff Λ in the original theory, and f ⊥, (log b) are the running couplings by solving the flow equations (3.45). If there exists an attractive fixed point for f ⊥, as we will see shortly, the integral is dominated by the fixed point value for a large Λξ 1 in a scaling regime. This motivates us to look at the flow equations for f ⊥, . From (3.45), we obtain It is not difficult to find a fixed point (f * ⊥ , f * ) that makes the right-hand side of the above equations vanish. To leading order in , let's ignore η ∼ 2 . First define a constant C by and inserting these back to (4.50) gives a self-consistent equation for C, that is, . (4.55) We find thatη ⊥ andη share the same critical behavior, and this justifies why we need to keep both. The critical exponents x λ and xη are defined by λ phys ∼ ξ x λ ,η ⊥, ∼ ξ xη , (4.56) and we have up to order , x λ = C , xη = (1 − C) , (4.57) which satisfies the well-known scaling relation in the model H; x λ + xη = − η [1,2]. For a comparison, the original model H has the value C = 18/19 ≈ 0.947. We find that d T = 1 case is somewhat close to the original model H, but for d T = 2 case the difference in C is about 10% which is significant. We recall that d T = 2 case is relevant for the QCD critical point in a background magnetic field. It is also easy to see thatη ⊥, grows much faster than λ , which justifies the approximation in obtaining the equation (3.24).
The relaxation frequency for the charge diffusion mode in hydro-regime k ξ −1 is given by where χ ∼ ξ 2−η is the susceptibility that is the inverse of the parameter r phys . Matching to the critical regime near k ∼ ξ −1 , we have up to order Since the charge diffusion modes relax more slowly than the shear modes, they define the critical slowing down. This gives the dynamic critical exponent [1,2] z = 4 − x λ = 4 − C . (4.61) With this, the flow equation for λ in (3.45) has a finite fixed point for λ * < ∞. For d T = 2 case, we have z ≈ 3.15 when = 1 (three dimensions). This is notably larger than the original model H value, z ≈ 3.05.
As a further direction, one could try to compute other refined quantities in these models such as scaling functions. It should also be possible to go to a next order in expansion using the method in Ref. [11] Nuclear Physics, within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration.