Scaling Relations in Two-Dimensional Relativistic Hydrodynamic Turbulence

We derive exact scaling relations for two-dimensional relativistic hydrodynamic turbulence in the inertial range of scales. We consider both the energy cascade towards large scales and the enstrophy cascade towards small scales. We illustrate these relations by numerical simulations of turbulent weakly compressible flows. Intriguingly, the fluid-gravity correspondence implies that the gravitational field in spacetimes with anti-de Sitter asymptotics should exhibit similar scaling relations.


I. INTRODUCTION
Turbulence is a ubiquitous phenomenon of fluid flows which plays a key role in many physical scenarios. At a broad level, turbulence takes place when non-linear interactions of a large number of degrees of freedom dominate over dissipative effects (the so-called high Reynolds number regime). Due to its intrinsic strongly non-linear and far from equilibrium character, a thorough understanding of this phenomenon from first principles remains elusive. The goal of understanding goes well beyond academic interests, as a deeper grasp of this phenomenon would impact a broad range of areas including weather dynamics, astrophysical processes, aerodynamics, etc.
Making this enterprise difficult, as thoroughly discussed in e.g. [1][2][3][4][5][6], is turbulence's chaotic nature, the flow of energy towards smaller or larger wavelengths, and its non-linearity. Common approaches in the analytical study of fluid turbulence rely on dimensional and statistical arguments, often assuming as many statistical symmetries as are possible. These include rotational, parity, and translational invariance, as well as stationarity in time. These efforts have been aided and complemented by numerical and physical experiments which provide important clues as to the extent to which such analytical results are robust with respect to departures from these simplifying assumptions.
To date, the majority of our understanding of fluid turbulence is for non-relativistic, incompressible fluid flows described by the Navier-Stokes equation. The relativistic regime, which is necessarily compressible, has received less attention. Nevertheless, many applications of interest naturally require its consideration. Examples include astrophysical fluid flows (e.g. [7]), as well as applications of the fluid/gravity correspondence (see [8] and references therein). This correspondence indicates that the behaviour of large black holes in asymptotically Anti-deSitter spacetimes disturbed by long wavelength perturbations can be studied by considering the relativistic hydrodynamics of a conformal fluid. In particular, the correspondence relates the fluid stress tensor to the asymptotic behaviour of the spacetime metric, as well as to intrinsic and extrinsic geometrical data of the black hole horizon. Thus, the understanding of turbulent relativistic fluids bears relevance also to the study of gravity.
In the present work, we add to a handful of steps already taken to understand relativistic turbulence [9][10][11][12][13][14][15][16] by performing both analytical and numerical analysis of such fluid flows, placing particular emphasis on two spatial dimensions. In part, we build upon previous work by Fouxon and Oz [9], who derived some scaling relations for relativistic hydrodynamic turbulence applicable to d ≥ 3 spatial dimensions. Firstly, we present some useful remarks regarding the special case of spatial dimension d = 2, together with new scaling relations in this case in both the inverse-and direct-cascade ranges. Secondly, we describe the current state of our numerical simulations of forced turbulence on a toroidal spatial domain, with the full spacetime topology being given by T 2 × R.
This work is organized as follows. Sec. (II) describes background information on energy scaling and velocity correlations which are standard results in non-relativistic turbulent fluids. Sec. (III) provides a discussion of some analogous concepts in the relativistic case, and derivations of scaling relations for this regime, with a particular attention to the dependence on dimensionality. New relativistic scaling relations will be derived for the hydrodynamic stress-energy tensor and vorticity both in the inverse energy cascade (3.20) and in the direct enstrophy cascade (3.28), (3.34) (see also (D15)). They reduce in the non-relativistic limit to known scaling relations of incompressible fluid turbulence. Sec. (IV) describes the numerical implementation employed, the initial conditions adopted as well as the statistical properties of the resulting weakly-compressible turbulent flow. We illustrate such compressibility through Fig. (1), which shows that both absolute and relative velocities are on the order of 20% of the speed of sound. In this regime, we demonstrate that the term ρ ργ 2 v L in Eq. (4.6) is highly sensitive to compressive effects (see Fig. (5) and (6)), at least insofar as it has a much wider probability distribution than its incompressible counterpart v L . This term might give a non-negligible contribution upon an increase in sample size since it cannot be argued to vanish by statistical symmetries, unlike its incompressible counterpart. We summarize in Sec. (V), and we have also included relevant information in the appendices, which we hope will prove useful for newcomers to the subject.
In this work, letters in the beginning of the alphabet {a, b, c...} will denote spacetime indices, while those beginning from {i, j, k...} will denote purely spatial ones. We adopt a Minkowski metric with signature (−, +, +, +), and we will either denote spatial vectors with a bold symbol r or with index notation r i , where appropriate. Furthermore, square brackets [.] will be used in Sec. (II A) to refer to a quantity's units. Angle brackets . will denote ensemble averages. Finally, we use units in which the speed of light c = 1.

II. BACKGROUND: NON-RELATIVISTIC FLUID TURBULENCE
The characteristics of turbulence are most cleanly studied within inertial ranges, which are length scales far from any friction, forcing, or viscous scales. In inertial ranges, the transfer of an inviscidly conserved quantity is independent of scale. Consequently, key aspects of the analysis are often simplified. One illustration is the possibility of using simple dimensional arguments to derive the famous Kolmogorov scaling as described in Sec. (II A).
Of particular interest for our discussion is the observation that the number of distinct inertial ranges that can exist has a dependence on dimensionality. In spatial dimensions d > 2, if energy is being injected at a large scale L f and is being dissipated (e.g. by viscosity) at a small scale L ν , then there will be an inertial range at length scales L such that L ν L L f for which the rate of energy flow between scales is independent of scale. On the other hand, for d = 2, there exists an additional inviscidly conserved quantity called enstrophy which gives rise to a second inertial range [5]. In what follows, we discuss some classic results in these ranges for d = 2 which are particularly relevant for our discussion (the interested reader may consult e.g. [17] and references therein for further details).

A. Energy scaling
In the inertial range for energy in d = 2, the specific kinetic energy E = v 2 /2 transfers preferentially toward larger length scales (since this behaviour is opposite to the d > 2 case, this inertial range is referred to as the inverse-cascade range). On the other hand, the specific enstrophy, defined by Z = ω 2 /2 (where ω = ∇ × v is the vorticity, a pseudoscalar quantity in d = 2), is associated with the direct-cascade range, so-named since it transfers preferentially toward smaller length scales. In these ranges, power-law scaling of the specific energy spectrum E(k) can be obtained by dimensional analysis [17] as reviewed below. (Note that unless otherwise specified, energy and enstrophy are given per unit mass).
The energy spectrum has units of energy per wavenumber, or Let us restrict to the energy inertial range where the rate of energy transfer through scale k is given by = (k). Now, assuming that the only relevant scales are and k, the ansatz E(k) ∝ p k q for some p and q allows us to solve for the powers p, q through a comparison with Eq. (2.1): This yields p = 2/3, q = −5/3, which is the famous Kolmogorov scaling (see e.g. [17]), This scaling, theoretically obtained for all dimensions, has been reported in early experiments of 3D turbulence, such as in a jet of air under laboratory conditions [18], and in effectively-2D turbulence, such as in planetary atmospheres [19], electromagnetic-layer experiments where a thin layer of electrolyte is externally forced by magnetic fields [20], etc. Numerical experiments have also shown such behavior in, e.g. simulations of forced, steady-state, 2D incompressible Navier-Stokes tubulence [21]. It is known to be violated, however, both numerically and experimentally in d > 2, leading to an anomalous scaling exponent (see the discussion in [22] and references therein e.g. [23]).
In the special case of d = 2, a further relation can be obtained which is valid in the direct-cascade range. Here, the transfer of enstrophy η towards small scales is independent of scale. Using an analogous ansatz in this range, E(k) ∝ η p k q , gives This yields p = 2/3 and q = −3, thus giving a different scaling of the energy, for which there is a dimensionless multiplicative logarithmic correction which we will not discuss (see [17]). This scaling has been observed simultaneously with the −5/3-scaling in both 2d turbulence in a soap film [24], and more tentatively in the limit of very high spatial resolution in numerical simulations [25].

B. Velocity structure functions
Another classical result in the theory of turbulence is the scaling of velocity structure functions in the 2D inverse cascade range, which highlights important statistical correlations in a turbulent flow. A velocity structure function of order n is a Galilean invariant, defined as where eachê i is a unit vector oriented in some fixed direction with respect to the spatial separation r, and the angle brackets . denote an ensemble average. In statistically isotropic conditions, it suffices to consider only longitudinal and transverse directionsê L ,ê T , which are parallel and perpendicular to the vector r, respectively (see [26]). For brevity, let us define δv (r) = (v(r) − v(0)) ·r and δv ⊥ (r) = (v(r) − v(0)) ·r ⊥ , where isotropy now serves not only to make the notion of transverse unambiguous, but also implies these quantities depend only on the distance r = |r|. For non-relativistic, incompressible turbulent flows, one can derive scaling relations for velocity structure functions by introducing a statistically homogeneous, isotropic, random external force. The external force helps to establish the inertial ranges, and its statistical characteristics allow for a clean calculation (see [27] for the d = 2 case). The force has an energy injection rate I , in terms of which one finds in the inverse-cascade range r L f orcing in d = 2 that δv (r) If we suppose that this relation implies a scaling of the individual velocity differences as δv ∝ r 1/3 , then this immediately implies a scaling for all orders of structure functions S n (r), with any mixing between longitudinal and transverse components, given by provided only an even number of transverse components appear (see Appendix (C) for an elaboration of this point). This general scaling has been observed in various experiments [20,28,29], as well as in forced 2D Navier-Stokes turbulence [21]. Note that the scaling in Eq. (2.8) is known to be violated in all direct cascades, except for the n = 3 structure function (see [22] and references therein, e.g. [23]). We stress that this is by no means a complete list of references, and the interested reader should see [17] for a survey of previous work.

III. RELATIVISTIC HYDRODYNAMIC TURBULENCE
We now turn our attention to the case of interest, namely relativistic hydrodynamics. Let us concentrate on the equations of motion given by the conservation of the stress-energy tensor T ab , where a, b, c . . . are spacetime indices ranging from 0 . . . d, with d the spatial dimension. Our goal is to derive the scaling behaviour of correlations which are analogous to those found in Navier-Stokes turbulence.

A. Relativistic relations I: Fouxon and Oz derivation
For the sake of our presentation, we now reproduce in a more detailed manner the derivation of the scaling relations presented by Fouxon and Oz [9] for the particular case of relativistic hydrodynamics (see also [30] for compressible non-relativistic turbulence). Our notation, however, will differ: quantities evaluated at the point r 2 will have a prime, while quantities evaluated at the point r 1 will not.
As for Navier-Stokes turbulence, by including a random, homogeneous, and isotropic external force in the equation of motion, the inertial regime can be explored. We begin with and assume a steady-state condition, with no sum on i. We stress that this condition is stronger than if we were to sum over i, since in the single-point limit it enforces the stationarity of average momentum T 0i in separate directions individually, whereas summing would enforce the stationarity of the total. Acting with the derivative gives Notice that interchanging the points 1 and 2 amounts to inverting the spatial coordinate axes, but this leaves the product T 0i T 0i unchanged. This can be easily seen by considering a perfect fluid T ab = (ρ + p)u a u b + pη ab , expressing where v i is the spatial velocity, and realizing that the switch changes the sign of each v i but the product vv remains unchanged. We can therefore switch points 1 and 2 in the first term of Eq. (3.4) without consequence, which combines the two terms to give where we have used Eq. (3.2) to replace the time derivative, and a sum on j is understood. It is important to stress that the spatial derivative here is with respect to the coordinates of point 2 (we have denoted this with a prime, ∂ ). This means that it views functions of r 1 as constant, and so can be brought out of the ensemble average. Thus, Now, notice that homogeneity implies that these averaged quantities are functions of the separation r ≡ r 2 − r 1 only, so ∂ j ≡ ∂/∂(r 2 ) j = ∂/∂r j ≡ ∂ j when acting on them. This gives, Assuming that r L f ≡ the correlation length of the forcing allows the approximation i , which is now constant with respect to r. Using the fact that the left-hand side is a gradient and T 0i (t)T ij (t) is isotropic (not a function of angle), Eq. (3.6) can be integrated over a disc using the divergence theorem. This yields This completes the derivation. Assuming instead that r L f does not allow one to integrate Eq. (3.6) without more information, as the result would depend on the details of the forcing at all scales up to r. Also, note that if one wished to enforce that f i is divergence-free, then f i would only be isotropic in the sense that f i (k)f i (k) is a function of k only, so one should sum over i in Eq. (3.7) in that case.

B. Relativistic relations II: the case of d=2
In this section we concentrate on the behaviour of T 0i (t)T ij (t) for the special case of d = 2 in the inverse-cascade range, as well as an additional correlation function in the direct-cascade range which involves a quantity resembling vorticity.
A special treatment of d = 2 is required, since whether the steady-state condition (3.3) is appropriate depends upon whether the energy injected by the external force can be removed. For d > 2, since injected energy transfers to small scales, it will encounter the viscous scale and be dissipated. There is evidence that this behavior persists even for arbitrarily small viscosity, and is known as the energy dissipation anomaly (see [6]). This can be understood heuristically as a result of the direct cascade of energy; a finite viscosity, no matter how small, will produce strong energy dissipation below the viscous scale, and the direct cascade of energy guarantees that this scale will eventually be encountered. Mathematically, this can be understood in the incompressible case as a result of the unboundedness of enstrophy. One can derive the energy balance equation with no forcing or friction [6], where E ≡ v 2 /2 is the mean energy and Ω ≡ ω 2 /2 is the mean enstrophy.
If Ω can become comparatively large as the viscosity ν becomes small, then the right-hand side of Eq. (3.8) can remain non-zero. The balance equation for Ω contains a source term which is due to vortex stretching, preventing one from bounding its growth. Thus, if the anomalous energy dissipation persists in the relativistic regime for d > 2, then the energy injected by the external force can be removed and the steady-state condition (3.3) is appropriate.
On the other hand, the situation is different when d = 2. Vortex stretching is absent in this case, which means no source term appears in the enstrophy balance equation, thus being given by [17] where P = |∇ω| 2 /2 is the mean palinstrophy, and where we have again restricted to the incompressible case with no forcing or friction. Eq. (3.9) says that the mean enstrophy is dissipated in time, which means it is bounded from above. It follows that the energy dissipation vanishes in the limit of zero viscosity, since the enstrophy cannot grow comparatively. If this fact remains true in the relativistic regime, then the steady-state condition (3.3) is inappropriate for d = 2 in the inviscid limit. Even without taking that limit, when the forcing and viscous scales are sufficiently separated, one would expect energy dissipation to be small (and to remain small over time, due to the inversely-cascading energy), thus easily failing to balance the injection of energy. Thus, in the next section we present an alternative derivation which gives rise to different scaling behaviour of the same correlation functions. (It is worth mentioning large-scale energy may transfer towards the viscous scale through the formation of shocks or large gradients, where it would be dissipated. This might allow the steady-state condition (3.3) to hold in d = 2.) Lastly, in d = 2 the balance equation for the palinstrophy P reads which has a source term of indefinite sign. This means that the palinstrophy P cannot be bounded from above, so there may be an enstrophy dissipation anomaly in d = 2 [17,27] if P can become large enough that the right-hand side of Eq. (3.9) remains non-zero for arbitrarily small ν. This will allow us to derive new correlation functions in the relativistic case by considering a different steady-state condition involving quantities that resemble vorticity.

Scaling in the inverse-cascade range
We now derive a relativistic scaling relation in the inverse-cascade range by adapting a strategy used in the incompressible Navier-Stokes case [27]. Let us begin by defining a quantity by where we are summing over i this time. What follows does not require to be independent of time. Consider a new form of stationarity, weaker than Eq. (3.3), which is consistent with a lack of removal of energy, Notice that expression (3.12) reduces in the Newtonian limit to the stationarity of a second-order velocity structure function. Now, homogeneity implies T 0i T i 0 = T 0i T i 0 , and recall that we have already evaluated the third term on the right-hand side of this equation in Eq. (3.5). Upon replacement, one obtains At this point, we must relate T 0i f i to our choice of external force. We adopt a divergence-free homogeneous Gaussian random field with zero mean, characterized by its two-point correlation function (see e.g. [31]), such that F ij decays rapidly beyond the forcing scale L f . We impose isotropy in the sense that trF ≡ F i i is a function of r only and, as shown in Appendix (A), trF = 2 T 0i f i , which gives us greater control over this term. Furthermore, since f i is divergence-free, trF = ∂ k Θ k for some appropriate Θ k . This can be seen by first noting that f i , if divergence-free, can be written in terms of a stream function ψ as f i = ij ∂ j ψ. Thus, We can then relate this to trF by integrating Eq. (3.14) with respect to time to eliminate the δ-function, then defining Θ k ≡ dτ ( ik ψ)( i n ∂ n ψ ) , thereby showing trF = ∂ k Θ k by construction. Our expression Eq. (3.13) thus becomes which, under the isotropic conditions assumed, integrates to a result which holds for all r. For the inverse-cascade range, this result further simplifies since Θ j is negligible there. To see this, first note that for j = T the transverse direction, T 0i T i j vanishes by isotropy and r j vanishes by definition, so Θ T must also. For Θ L , recall that we have stipulated that trF decays rapidly beyond the forcing scale L f . Thus, integration of trF over a disc of radius r will approach a constant as r exceeds the forcing scale L f , whereas applying the divergence theorem yields Thus the longitudinal component Θ L decays at least as quickly as 1/r, becoming negligible at large distance. Consequently, in the inverse-cascade range, where we have neglected the subleading term. Notice that this result has the opposite sign with respect to the d > 2 case, which in the incompressible limit is known to reflect the inverse cascade of energy. As a word of caution, note that this scaling is usually presented as positive since the points r 2 and r 1 are switched.
In other words, Eq. (3.19) is equivalent to Thus, when comparing the overall signs in Eqs. (3.7) and (3.19) with the literature, one should be mindful of this point.
It is interesting to note that in the incompressible limit, the constant = v i f i is the lowest order term in the Taylor series of v i f i [27]. Thus, in the short-distance limit, the first non-zero term in the Taylor series of Θ j − r j is proportional to r 2 r j . This gives the cubic scaling of the third order velocity correlation familiar from the statistical theory of incompressible turbulence. However, incompressibility plays a crucial role in this result, so we cannot make a similar inference in the relativistic case without additional assumptions.
Finally, had we used the slightly weaker steady-state condition that ∂ 0 (T 0i − T 0i )(T i 0 − T i 0 ) = constant not necessarily zero, we would clearly still obtain linear scaling in the inverse-cascade range, although with a different proportionality constant. This weaker assumption might hold on a periodic 2D spatial domain, such as a torus, in the absence of any removal of energy. However, since energy would cascade towards the longest available length scale, anisotropy would grow as energy condensates into the lowest mode (see Appendix (E) for a numerical simulation of this scenario). Thus, the linear scaling obtained here would be expected to hold only in the intermediate stage when the flow is still isotropic.

Scaling in the direct-cascade range
In the incompressible, non-relativistic case, the statistics in the direct-cascade range can be cleanly studied using a steady-state condition of correlations involving the vorticity. A similar strategy can be adopted here, although subtleties arise with regard to the precise expression of vorticity adopted. In what follows we describe what we consider the most straightforward path and refer to Appendix D for a related option. First, consider the spatial component of Eq. (3.2), and apply the 2-dimensional curl to obtain The incompressible limit of this is the standard equation for vorticity. However, it is interesting to note that Eq. (3.22) does not describe what is normally regarded as the relativistic vorticity, even though it has the same incompressible limit. (We describe the behaviour of the relativistic vorticity in Appendix D, together with a mention of the subtleties related to deriving scaling relations with it.) We identify the right-hand side of Eq. (3.23) as the curl of the external force, which we will denote as F. For brevity, let us also define the first two quantities in brackets as ω = ik ∂ k T 0i andω j = ik ∂ k T ij , giving the suggestive expression We may now multiply this expression by ω and take the ensemble average, which gives, We have explicitly shown that the spatial derivative is with respect to the point r 1 . It acts on a correlation which, by the assumption of homogeneity, is a function of separation r = r 2 − r 1 only. Thus, we can change ∂/∂(r 1 ) j to −∂/∂r j , which herein we write simply as ∂ j , thus giving Assuming the existence of a dissipation anomaly for the quantity ω 2 /2 , which would balance the injection from the external force, we can impose the steady-state condition even for arbitrarily small viscosity. Thus Eq. (3.25) yields In the direct-cascade range r L f , Fω ≈ Fω ≡ ε, which allows us to integrate Eq. (3.27) using isotropy, obtaining Summarizing so far, we find that T 0i T i j scales linearly in the inverse-cascade range with the opposite sign relative to the d > 2 case, and its linear scaling in the inverse-cascade range ought to be robust with respect to the background topology, subject to the assumption of isotropy. Furthermore, we found that ω j ω scales linearly in the direct-cascade range.
Finally, it is possible to integrate Eq. (3.28) twice more to obtain a cubic scaling of T 0T T LT in the direct-cascade range, but through this procedure one obtains no information about the purely longitudinal correlation T 0L T LL . To see this, begin by writing the left-hand side of Eq. (3.28) withω j and ω appearing explicitly in terms of the stress-energy tensor, where we have again used ∂/∂(r 2 ) i = −∂/∂(r 1 ) i = ∂/∂r i ≡ ∂ i in the last line, which is true when the derivative acts on functions of the separation r = r 2 − r 1 only. For cleanliness, define A n j ≡ ∂ k mn ik T mj T 0i , so that Eq. (3.28) now reads We wish to integrate this over a disc using the divergence theorem, so let us obtain a scalar equation by projecting the j-index onto the longitudinal direction L, giving Integration over a disc, assuming isotropy so that A L L = A L L (θ), yields A further application of the divergence theorem yields Using the identity mL iL = δ mi δ LL − δ mL δ iL , one obtains the final result, valid in the direct-cascade range.

IV. IMPLEMENTATION DETAILS
In order to test the derived scaling relations, we numerically implement the relativistic hydrodynamical equations subjected to an external force with suitable statistical properties. We then extract relevant quantities from the numerical solution, as described below. In what follows we provide details of our implementation.

A. Flux-conservative formulation
For convenience we express the equations of motion in flux-conservative form. In the absence of driving-sources, this helps to ensure energy-momentum conservation at the discrete level. As discussed in [32,33], the combination of discrete operators satisfying summation by parts together with a Runge-Kutta integrator of third order guarantees an energy conserving scheme. Eq. (3.2) gives two expressions, already in the desired form, where i and j are spatial indices. These equations fully determine the system in the ultra-relativistic regime where the conservation of particle number becomes irrelevant at the classical level. We take {T 00 , T 0i } to be our set of conservative variables, and evolve them directly. Using a perfect fluid with the conformal equation of state p = ρ/2, the equations of motion become where we have used u a = γ(−1, v). We define our conservative variables as D ≡ (ρ/2)(3γ 2 − 1), S i ≡ (3/2)ργ 2 v i and our primitive variables as (ρ, v i ) for i = 1, 2. Note that the second equation provides the time-evolution of S i , which then sources the first equation for the time-evolution of D. The forcing function f i , which is completely spatial, is described in Appendix (B). The transformation from conservative variables (D, S i ) to primitive variables (ρ, v i ) is given by Solving for the Lorentz factor in terms of the conservative variables amounts to solving a quadratic equation for γ 2 .
The presence of ρ in the denominator presents a potential problem in the recovery of v i when ρ = 0. In general applications, this technical issue can be circumvented by artificially maintaining a non-zero floor or atmosphere for ρ, small enough so as not to affect the dynamics appreciably. However, in our simulations the density never reaches zero, so this mechanism is never invoked.
B. Spatio-temporal reduction of the ensemble average .
It is often impractical to calculate . as an ensemble average. In practice, one exploits statistical symmetries and the assumption of ergodicity to reduce . to a spatial or temporal average. For instance, for a statistically homogeneous and isotropic flow, one computes . as an average over pairs of points with a scalar separation r = | r|. Further details on the mathematical subtleties involved in doing these reductions rigorously are given in [6], and will not be discussed here.
Taking Eq. (2.6) as an example, the homogeneous and isotropic averaging process means that all quantities in the product are projected onto directions defined relative to the separation vector r. Thus, when computing the average spatially on a numerical grid, although r itself may vary in direction from term to term, the relative directions between r and the projection directions must remain the same.
With this in mind, we understand that the spatial indices in Eq. (3.7) stand for projections onto those directions. This implies that when the j-th direction,ĵ, is perpendicular with r, the correlation in Eq. (3.7) vanishes since r j = r ·ĵ = 0. This is a consequence of isotropy, and we elaborate on it in Appendix (C).

C. Numerical experiments
Our simulations take place on a torus and, as mentioned, unless some some form of large-scale extraction of energy is employed, energy will build up in the longest mode. To address this issue, we adopt a convenient approach by augmenting Eq. (4.2) with a linear 'friction' term −αT 0i on the right-hand side 1 , giving The friction term causes the system to evolve towards an approximately constant total energy. For sufficiently small α > 0, this final state exhibits inertial range scaling in the Newtonian spectral energy E(k), and thereby exhibits fully developed turbulent behaviour. By dimensional analysis, α can be related [17] to the energy injection rate I and the friction length-scale L α through α ∼ For a given I , we choose α such that the friction length scale L α is a few times smaller than the spatial extent of the domain. We next describe the setup of our simulation of fully-developed, steady-state turbulence, described by Eq. (4.4).
The initial conditions adopted are {ρ = 1, v i = 0}, and the uniform spatial grid has N 2 = 800 2 points (which admits the Nyquist wavenumber k max = 400, expressed in grid units; one can convert to real units via 2πk max /L) and periodic boundary conditions are imposed. For concreteness, all reported times will be given in multiples of the light-crossing time t LC = L/c, with L ≡ the size of the box, which we set numerically to 10. The random external force employed is described in Appendix (B), and with its strength controlled by the parameter Ψ(0) = 3 × 10 −5 we find a suitable value of the friction strength to be α = 1.8 × 10 −2 , producing a large-scale energy cutoff around k = 5, as shown in Fig. (3). Fig. (1) displays the velocity and density distributions at a representative time when the flow appears statistically stationary, as well as the standard deviation of the velocity differences δv L , δv T as a function of separation. This is intended to convey the degree to which this flow differs from the non-relativistic, incompressible case. In particular, notice that the peak of the velocity distribution at v ∼ 0.14c corresponds to a Lorentz gamma factor γ ∼ 1.01, while the largest velocity is v ∼ 0.52c, corresponding to γ ∼ 1.17. The bulk of the flow can therefore be considered nonrelativistic (for comparison, the sound speed for this 2 + 1 dimensional conformal fluid is c s ∼ 0.71c, so the flow is also subsonic). The density distribution shows a standard deviation of 0.055 and a peak at 0.97. The velocity differences are roughly Gaussian distributed with zero mean, and their widths σ are comparable to the sound speed. One may thus describe this flow as being in the weakly-compressible regime. Notice that the density distribution peaks at a value less than 1 and has a stronger tail at lower values, which means that a bias is formed in favour of under-density with respect to the initially uniform value of 1. Given the characteristics of the flow described, employing the Newtonian energy and enstrophy to connect with known results in the Newtonian regime is justified.
The total specific Newtonian kinetic energy of a representative simulation is shown as a function of time in Fig. (2). The energy plateaus after approximately t = 20t LC , indicating a statistically steady-state. Correlation functions are computed at t = 25t LC . In order to obtain snapshots of the flow which are statistically independent, one can choose the temporal spacing between samples to be at least one large-eddy turnover time, determined through T = U/L, where U is a typical large-scale speed and L is the large length scale. We estimate U by applying a low-pass Fourier filter to the velocity field at a representative time t = 25t LC , with all wavenumbers larger than the friction scale k α being set to zero, then choosing U as the mode of the resulting velocity distribution. The large length scale L is chosen as 2πL α . We find this procedure gives roughly T = 25t LC , which is the same amount of time required to evolve the fluid from rest to a steady-state. Thus, we opt to evolve the fluid from rest to obtain each flow realization, rather than evolve from a steady-state at time t to a later time t + T . This reduces the risk that each flow realization is not statistically independent.
The spectral energy E(k) and flux Π(k) (averaged over 200 flow realizations) are shown in Fig. (3). Π(k) is computed using the formula Π , as described in [6]. Here, the superscripts >, < denote Fourier-filtered quantities with all wavenumbers set to zero below or above the given k, respectively. The familiar Kolmogorov scaling of the spectral energy E(k) ∝ k −5/3 seems to hold, and the spectral energy flux exhibits qualitative behaviour similar to that displayed in [25]. The intersection of Π(k) with the horizontal axis at k = 100 indicates the injection of energy there, since energy is flowing away from that length scale. Negative values of Π(k) for k < 100 indicate the inverse-cascade of energy.
To provide evidence that the flow is indeed statistically isotropic, we compute 2nd-order velocity correlation functions of purely longitudinal and mixed types, the latter of which is expected to vanish under isotropic conditions. The mixed correlation δvLδvT is supposed to vanish under isotropy, and it is measured to be less than the non-vanishing purely longitudinal correlation δv 2 L by a factor of more than 10 3 . All correlations are computed over 10 4 flow realizations. Right: Numerical evidence that fif i and fiT i 0 are proportional, in particular vanishing quickly with increasing r. Note that the former has been scaled by an overall constant in order to match with the latter at r = 0. We show only the longitudinal correlations here -the transverse ones look the same. Figure. (4) (left) illustrates the results obtained with an average over 10 4 flow realizations, which shows the mixed type being negligible with respect to the purely longitudinal case. In addition, the right plot demonstrates that we have achieved f i f i ∝ f i T i 0 , which is a crucial part of the derivation in Sec. (III B 1) of the linear scaling of T 0i T ij in the inverse-cascade range. However, the proportionality factor is on the order of 10 3 rather than 2 as the argument in Appendix (A) would suggest. This discrepancy is not surprising, as our force is not δ-correlated in time as the argument in Appendix (A) requires. A proper numerical implementation of δ-correlated statistics requires a modified integration algorithm, as described in [34].
As a further display of the properties of this flow, we also compute velocity structure functions of orders 1 through 4, but with the absolute value of the velocity differences taken in the case of odd orders (this has been argued [35] to preserve the scaling properties, though it obscures the overall magnitude of the correlation). By taking the absolute value, all contributions to the correlation add constructively, which improves the convergence drastically (this is what was done in [13] for a relativistic fluid in d = 3). The scaling behaviour reflects the Kolmogorov expectation |δv L | n ∝ r n/3 increasingly poorly as n increases, but the same phenomenon has also been reported in [25] for positive-definite velocity structure functions.
Lastly, it is interesting to closely examine the non-vanishing correlations, T 0L T LL and T 0T T LT , together with their incompressible limits, (9/4) [17] for details about this expectation). The odd orders have an absolute value operation performed on δvL, which causes them to converge more rapidly with sample size N . Right: The two relativistic correlations T 0L TLL and T 0T TLT which are not expected to vanish under isotropic conditions, along with their incompressible limits (9/4) v L vLvL and (9/4) v T vLvT . The factors of 9/4 are left over after the limit is taken. All correlations are computed over 10 4 flow realizations with a grid size of 800 2 .

D. Discussion
The correlations shown in Fig. (5) (right) are still unresolved, even with a sample size of 10 4 flow realizations. We observe this by computing, at each r/L, the standard deviation of the sample divided by the square root of the number of samples, σ/ √ N . We find it to be comparable with the value of the correlation itself, and thus conclude that the fluctuations have not averaged down sufficiently. Poor signal-to-noise constitutes the main difficulty in numerically measuring odd-order correlations in compressible turbulence. For our current simulation, a much larger sample size is required, as we will discuss later.
Nevertheless, some relevant conclusions can still be made. Firstly, notice in Fig. (5) (right) that there is a great disparity in how well T 0T T LT and T 0L T LL match with their incompressible limits. The former is indistinguishable from its limit in the plot, though zooming in reveals that there are small differences. The latter, on the other hand, bears little resemblance to its incompressible limit. To gain insight about this, consider the two correlations written in terms of the primitive variables: In the incompressible limit, γ, ρ → 1. The 2nd term on the right-hand side of Eq. (4.6) will therefore become ∝ v L which is zero by statistical symmetries. Indeed, we find numerically that the overall magnitude of 3 4 ρ ργ 2 v L is roughly 10 5 times larger than v L . This indicates that the underlying probability distribution for this term is highly sensitive to compressive effects, which cause it to become considerably wider, translating into much larger fluctuations. We can also implicate this term in the disagreement between Eq. (4.6) and its incompressible limit by subtracting it from Eq. (4.6). We display the result in Fig. (6), where it is seen that the agreement improves considerably. It remains to be seen whether this spoiler term will average down to become negligible in the weakly compressible regime we are exploring here.  6. The purely longitudinal relativistic correlation plotted without its spoiler term 3 4 ρ ργ 2 v L , as compared with its incompressible limit. The agreement is considerably improved, demonstrating that the spoiler term is sensitive to compressive effects.
Secondly, note that the third-order velocity correlation v L v L v L has been successfully resolved in simulations of an exactly-incompressible fluid in [21] when averaged over only tens of flow realizations, albeit with a larger grid size of 2048 2 . As a proof of principle, switching to the less-costly case of an exactly-incompressible fluid 2 we obtain similar results for our current grid size of 800 2 , shown in Fig. (7) after averaging over 7 × 10 4 flow realizations. An investigation into the dependence of the signal-to-noise of the correlations on compressive effects and the nature of the random external force is left for future studies. For such work, it is important to estimate the sample size required to resolve a given correlation. Such an estimate can be obtained in terms of the standard deviation of the underlying distribution and the scaling prediction. For instance, Eq. (3.20) provides the prediction for the strength of the signal. Supposing the underlying distribution for T 0i T i L | r=r I has a standard deviation σ(r I ), where r I is a separation within the inverse-cascade range. Then the signal-to-noise ratio SNR would be given by where N is the sample size. Solving for N yields The flat interval on the right plot also corresponds to a linear trend, which is a similar result to that of [21]. Negative values have been indicated in red on the left. Bottom: the sample average and uncertainty σ/N 1/2 plotted versus sample size N for the purely longitudinal velocity structure function (δvL) 3 |r=r I , for a separation rI in the inertial range (indicated with an arrow in the top row). This separation corresponds to a wavenumber of k = 35 in grid units. The uncertainty reduces to roughly 1/2 of the average at N = 7 × 10 4 , translating into a signal-to-noise ratio of ∼ 2.

V. SUMMARY
In this work we derived scaling relations in fully-developed relativistic turbulence in two spatial dimensions. We considered both the inverse-and direct-cascade ranges, and the relativistic results reduce in the non-relativistic limit to the corresponding scalings in the incompressible case. This derivation bridges known results in the field of incompressible fluid turbulence with ongoing work in the relativistic case.
We have also begun a numerical experiment in an effort to measure the derived scaling relations through direct numerical simulations. We showed through Fig. (1) that the flow displays Mach numbers around 0.2 in both absolute and relative velocities, and is thus weakly-compressible. In this regime, the probability distribution underlying the term ρ ργ 2 v L in Eq. (4.6) acquires a large standard deviation as compared with its incompressible counterpart, v L , being 10 5 times larger with the same sample size. While the latter can be argued to vanish by isotropy, the former cannot. This opens the possibility that this 'spoiler term' provides the dominant deviation of Eq. (4.6) from its incompressible limit to leading order in compressive effects, although this must be verified by increasing the sample size considerably. It would be interesting to observe deviations of the relativistic correlations derived here from their incompressible limits in the highly-compressible or relativistic regimes.
The signal-to-noise of odd-order correlations is the overarching difficulty in measuring them accurately. Indeed, previous studies of compressible Navier-Stokes turbulence (eg. [36][37][38]) and relativistic turbulence (eg. [13]) sidestepped this issue for the case of velocity structure functions by taking the absolute value of the velocity differences, |δv| n . The result of this procedure is that every term in the average adds constructively, and the scaling behavior has been argued in [35] to be preserved. In our case, an analogous work around is not available for the relativistic correlations in Eq. (3.20), since the coupling of factors of the velocity in T ij prevent writing the correlations in terms of velocity differences. We are continuing our effort to resolve the correlations, and those results will be left for a future communication.
As a final comment, our work examining the behaviour of relativistic, conformal fluids undergoing turbulence has a natural connection with both conformal field theories and gravity through holography and the fluid-gravity correspondence (e.g. [39][40][41]). The correspondence relates the fluid stress tensor in d dimensions to the asymptotic behaviour of a d + 1 dimensional black hole spacetime metric (up to counter-terms to obtain a finite expession) as, where γ ab and K ab are the intrinsic and extrinsic curvatures respectively of a timelike surface at r → ∞. It implies that in the turbulent gravitational regime, correlations involving the metric tensor itself should obey scaling behaviour of the form discussed here 3 . The implications of such an intriguing observation are still unexplored.
future time. For instance, the integrand on the right-hand side evaluated at a time τ cannot depend on the force at a time t > τ , and so the lower limit of integration cannot extend below t . Upon evaluating Eq. (A3) at equal time, the integral appearing on the right-hand side vanishes so long as the integrand remains finite, and one obtains Setting R[f ] = T 0i (r , t) and then using Eq. (A4) and Eq. (A1), one obtains the desired result, Appendix B: The forcing function We wish to construct a Gaussian random forcing function which is divergence-free and statistically homogeneous and isotropic, whose two-point correlation decays quickly with increasing distance, and which is delta-correlated in time (sometimes called white noise in time). In symbols, ∂ i f i = 0 and where F i i ≡ trF = trF (r). This latter condition is the sense in which this is an isotropic vector field. For f i divergence-free, the stronger form of isotropy where F xx = F xx (r) and F yy = F yy (r) forces F xx = F yy = 0. One can see this by moving to Fourier space, where the xx and yy two-point correlation functions become f x (k)f * x (k) ≡ g(k) and f y (k)f * y (k) ≡ h(k), for some functions g, h of the magnitude of the wavevector only. The divergence-free condition reads k xfx + k yfy = 0, which allows us to convert between these correlation functions. Thus, for k x = 0 which contradicts f y (k)f * y (k) ≡ h(k) unless h(k) = 0. This is why we chose to sum over i in Sec. (III B 1), since the divergence-free nature of the force played a role in the argument.
In our simulations, to generate a divergence-free force we derive it from a stream function ψ such that f = (∂ y ψ, −∂ x ψ). We thus specify ψ itself as a Gaussian random, homogeneous, isotropic scalar field which is deltacorrelated in time. In symbols, with Ψ(r) a thin Gaussian function, which ensures a short correlation length. In practice, ψ is built in Fourier space, where the reality conditionψ * (k) = ψ(−k) is imposed, and where each mode receives a complex amplitude drawn from zero-mean Gaussian distributions whose widths are given byΨ 1/2 (k). As constructed, ψ satisfies Eq. (B2). The force itself then has trF whose Fourier transform is a wide Gaussian weighted by k 2 . In real space, trF behaves as is plotted in Fig. (4) (right). At each step in the Runge-Kutta integration this procedure is repeated anew, thus giving different individual realizations of the random force. This is the sense in which we have approximated a delta-correlation in time. A proper implementation requires a modified algorithm, as described in [34].
Appendix C: The consequences of isotropy In their seminal paper, Karman and Howarth [26] argue as follows (we reproduce their argument for d = 2). Let the two points under consideration lie on the x-axis. We say that the x-direction is the longitudinal direction, pointing directly between the two points, while we say that all other perpendicular directions are transverse directions. The triple velocity correlation functions can be listed as Now, both directionsŷ and −ŷ are transverse, and by isotropy (or parity-invariance in d = 2) every correlation will be invariant under a switch between them. However, any correlation with an odd number of y-components will undergo a change of sign when the y-axis is inverted. Those correlations therefore vanish. The remaining correlations are The incompressibility condition i ∂ i v i = 0 further implies that only one of these remaining three correlations is independent (see [26]). In the relativistic case, there does not necessarily exist an analogous condition, so all three correlations might be independent.
The same arguments about sign-flipping apply in the case of homogeneous, isotropic turbulence in a special relativistic perfect fluid. For a correlation such as T 0i T ij , one can see for example with a perfect fluid energy-momentum tensor where u = γ(1, v), that for i = j we will have T ij = (ρ + p)γ 2 v i v j undergo a change in sign when one of the i-or j-axes is inverted. On the other hand, if i = j, then T ii = (ρ + p)γ 2 v i v i + pδ ii does not change sign when the i-axis is inverted. Furthermore, T 0i = (ρ + p)γ 2 v i changes sign when the i-axis is inverted. Thus all the facts are in place to run the same arguments presented in [26]. This allows us to conclude that the only non-vanishing correlations of this type are where L and T are the longitudinal and transverse directions, respectively.

Forced case
Suppose now there is a force acting in the problem, so that the conservation equation reads ∇ a T ab = f b . Projecting this equation along u a and orthogonal to it gives, Now, using the above relations, we have from Eq. (D1), We are interested now in exploring the condition ∂ a W a = 0 as in the previous section under the influence of a force. To proceed, let us observe that We thus arrive at the equation, This equation however does not relate the vector W a in the way we sought, i.e. an equation of the form ∂ a W a = F . Nevertheless, it does motivate what the right vorticity-related vector should be, namely W a ≡ adc P b d Ω bc which, from the derivation above, satisfies where we have used ρ = T 3 in the second line. Notice that with O(f ) denoting terms depending linearly or quadratically on f b , thus in the absence of forcing W a W a = Ω ab Ω ab and we recover the conservation of vorticity, Eq. (D5), as expected from this quantity.

Scaling Argument
A scaling argument in the direct-cascade range involving the relativistic vorticity W a can also be made, following Sec (III B 2) closely. First, define the right-hand side of Eq. (D11) as a forcing termF = acd ∂ a 1 3T 2 P b c u d f c , where we use a tilde to distinguish this force from the one in Sec (III B 2). This gives the equation of motion succinctly as ∂ a W a =F.
Second, consider the steady-state condition ∂ 0 W 0 W 0 = 0, and apply the time derivative and use the equation of motion to obtain, where ∂ i stands for the derivative with respect to the separation r = r 2 − r 1 . Lastly, notice that far below the forcing scale, r L f , the right-hand side is constant, F W 0 ≈ F W 0 ≡˜ , so upon integration (using isotropy) which is valid in the direct-cascade range. Notice that it is more difficult to integrate this expression twice than it is for the expression Eq. (3.28) due to the presence of the projector in the definition of W a , which prevents taking the derivative operator outside without picking up additional terms. This means that obtaining an r 3 scaling relation from this linear one is not as straightforward as in the case of Eq. (3.28).

Appendix E: Energy condensate
In the absence of large-scale removal of energy in 2D, energy will build up in the gravest mode. Such a state is called an energy condensate. For completeness, we present the energy condensate and a method for removing it from the analysis. As a concrete example, we adopt a periodic doman with grid size N 2 = 400 2 and a homogeneous, isotropic, random external force acting at k f = 50, normalized to a real-space amplitude β = 0.6. After a sufficiently long time, the inverse cascade leads to an energy condensate, as shown in Fig. (8). The figure displays the progression of the vorticity, with all times quoted in multiples of the light-crossing time t LC . The colour scale has been omitted. Notice the late-time appearance of two dominant vortices of opposing sign.
To analyze the resulting energy condensate one can make use of wavelets [45]. We perform a similar analysis here as far as decomposing the velocity field into coherent and incoherent parts and computing the spectral energy of each. Fig. (10) illustrates the obtained results. The decomposition is performed using Coiflet-12 wavelets, which are a complete set of functions which are localized in both real and Fourier space. Their first two moments vanish (as well as their third and fourth moments), thus they couple weakly to Gaussian features. In other words, a relatively large number of basis elements with relatively low weights are required to represent Gaussian features of the data, whereas non-Gaussian features are represented by fewer basis elements with higher weights. Assuming the incoherent part of the velocity field is closer to Gaussian than the coherent part, one can therefore extract the incoherent part by imposing a threshold on the field in wavelet space, setting to zero all wavelet weights above a certain value, and then transforming back to real space. The remainder is the coherent part.
To get a sense of what this procedure does, Fig. (9) displays such a decomposition of the vorticity at t = 960t LC using a threshold value of 3. Notice the increased blurriness of the coherent part of the vorticity (a common feature of compressed images, being represented by a small number of basis elements), and the dominant overall amplitude of the coherent part with respect to the incoherent part. The energy scalings displayed in Fig. (10) are approximately consistent with [45]. Note that there is a significant amount of arbitrariness in the choice of threshold value and wavelet type which we do not attempt to address here. Right: E(k) for the coherent part of the velocity field. The scaling behaviour of approximately k −3 , k −1 , and k −3 , respectively, is consistent with [45].