Numerical Measurements of Scaling Relations in Two-Dimensional Conformal Fluid Turbulence

We present measurements of relativistic scaling relations in $(2+1)$-dimensional conformal fluid turbulence from direct numerical simulations, in the weakly compressible regime. These relations were analytically derived previously in Westernacher-Schneider, Lehner, Oz (2015) for a relativistic fluid; this work is a continuation of that study, providing further analytical insights together with numerical experiments to test the scaling relations and extract other important features characterizing the turbulent behavior. We first explicitly demonstrate that the non-relativistic limit of these scaling relations reduce to known results from the statistical theory of incompressible Navier-Stokes turbulence. In simulations of the inverse-cascade range, we find the relevant relativistic scaling relation is satisfied to a high degree of accuracy. We observe that the non-relativistic versions of this scaling relation underperform the relativistic one in both an absolute and relative sense, with a progressive degradation as the rms Mach number increases from $0.14$ to $0.19$. In the direct-cascade range, the two relevant relativistic scaling relations are satisfied with a lower degree of accuracy in a simulation with rms Mach number $0.11$. We elucidate the poorer agreement with further simulations of an incompressible Navier-Stokes fluid. Finally, as has been observed in the incompressible Navier-Stokes case, we show that the energy spectrum in the inverse-cascade of the conformal fluid exhibits $k^{-2}$ scaling rather than the Kolmogorov/Kraichnan expectation of $k^{-5/3}$, and that it is not necessarily associated with compressive effects. We comment on the implications for a recent calculation of the fractal dimension of a turbulent $(3+1)$-dimensional AdS black brane.


I. INTRODUCTION
Relativistic hydrodynamics has become a subject of increased interest in recent years. Beyond its relevance in astrophysical scenarios (e.g. [2][3][4][5][6][7]), it has become relevant to the description of quark-gluon plasmas (e.g. [8,9]) and, through the fluid-gravity correspondence, it has found its way into the realm of fundamental gravity research [10,11]. Intriguingly, this correspondence has revealed that gravity can exhibit turbulent behavior, and studies of its possible consequences are gaining interesting momentum [12][13][14][15][16]. The understanding of turbulence in any regime is a difficult task given its intrinsic complexity, and despite a long history of efforts in the subject, our knowledge of this rich phenomenon is still incomplete. Important headways into this subject have been made thanks to statistical analysis complemented with numerical simulations (e.g. [17][18][19][20][21]). Of particular interest is to understand the possible onset of turbulence, and especially to derive scaling relations in fully-developed scenarios, since those cases are amenable to a statistical description from which generic statements can be drawn and then tested numerically and observationally.
To date however, only limited attention has been placed on the relativistic turbulent regime, and most of what is known restricts to the behavior of turbulent, incompressible flows in the non-relativistic regime. There has been some work on the analytical front [1,22,23] and several numerical investigations [1,13,[24][25][26][27][28][29][30][31]. Because correlation functions can indeed be measured in relevant scenarios -perhaps even in QG plasma [8,9,[32][33][34] -and interesting implications for the gravitational field follow from holography, it is of interest to further investigate relativistic turbulence.
In the current work we measure scaling relations in (2 + 1)-dimensional relativistic conformal fluids in the weaklycompressible turbulent regime and compare them to the predictions in [1] and various limits thereof. The (2 + 1)dimensional case is especially relevant to draw intuition for related phenomena in (3 + 1)-dimensional gravity with the help of the fluid-gravity correspondence (e.g. [12,16,25]). This work is largely a continuation and completion of [1], which contained a numerical study which was inconclusive at the time.
This work is organized as follows. Sec. (II) provides some background material, discussing both the inverse-and direct-cascade ranges that could ensue in fully-developed turbulence and, in particular, relevant scaling relations which we have measured. In Sec. (III) we provide details of our numerical implementations. This includes the use of a random external force to generate the turbulent flow, as well as considerations specific to simulating either a conformal fluid or an incompressible Navier-Stokes fluid. The equivalence between previously known results and the incompressible limit of the scaling relations derived in [1] is explicitly demonstrated. We give our results in Sec. (IV), where our numerical measurements of the scaling relations derived in [1] are presented. In Sec. (V), we provide additional discussion and ancillary numerical results, including the demonstration of a k −2 energy spectrum in the inverse-cascade of a turbulent conformal fluid which is not due to compressive effects.
Throughout this work, angle brackets . will refer to ensemble averages. Letters at the beginning of the alphabet (a, b, c, . . .) will represent spacetime indices (0, 1, 2), while letters in the middle of the alphabet (i, j, k, . . .) will represent spatial indices (1,2). We follow Einstein summation convention. In the context of correlation functions, which often depend on two points r 2 and r 1 , we define r = r 2 −r 1 . To avoid cumbersome notation, we denote quantities evaluated at r 2 with a prime (eg. T ij ) and quantities evaluated at r 1 without one (eg. T ij ). The metric signature is (−, +, +) for our (2 + 1)-dimensional setup.

II. BACKGROUND
We will make extensive connections with the work presented in [1], where specific scaling relations were derived analytically for (2+1)-dimensional relativistic hydrodynamic turbulence. We will compare these relations with suitable limits in order to make contact with previously known results, as well as to gauge the importance of relativistic vs compressible contributions in our simulations of the specific case of a conformal fluid.
A. Incompressible non-relativistic limit of the scaling relations In this section we explicitly demonstrate that the incompressible Navier-Stokes limit of the relativistic scaling relations presented in [1] can be written in terms of known results. We will use the particular form of a barotropic perfect fluid stress-energy tensor with equation of state P = ρ/w, where P is the pressure, ρ is the energy density, and w is the equation of state parameter. In doing so, we obtain incompressible counterparts to the relativistic scaling relations we measure in simulations, which act as a point of reference against which to gauge the relative performance of the relations derived in [1].

Inverse-cascade range
The first scaling relation, which is valid in the inverse-cascade range, reads [1] . For a P = ρ/w perfect fluid, where T ab = ((1 + w)/w)ρu a u b + (ρ/w)η ab , the scaling relation expands to where γ is the Lorentz factor and v i is the spatial velocity (u a = (−1, v i )). In the extreme incompressible nonrelativistic limit, ρ → constant and γ → 1. Thus ρ = ρ, and Eq. (2.2) becomes where we have defined NS ≡ ∂ 0 v i v i /2 as the incompressible Navier-Stokes version of . Note that the second term on the left-hand side vanishes due to statistical isotropy, yielding the final result Since Eq. (2.4) is the incompressible Navier-Stokes limit of the relativistic scaling relation in Eq. (2.1), they can be compared in the inverse-cascade range of relativistic or compressible turbulence in order to gauge their relative performance. Notice one can also arrive at Eq. (2.4) using known results in the theory of (2 + 1)-dimensional incompressible Navier-Stokes turbulence. In the derivations presented in [35], an intermediate result is displayed as valid in the inverse-cascade range. Here, δ denotes a difference, i.e. δv j = v j − v j . Expanding out the left-hand side of Eq. (2.5) and using statistical homogeneity yields By incompressibility, the second term on the right-hand side is divergence-free. Thus, assuming isotropy and regularity at r = 0, it must vanish [36] (this argument will be used repeatedly in Sec. (II A 2)). Therefore, Eq. (2.5) becomes which is the incompressible non-relativistic limit we obtained in Eq. (2.4).

Direct-cascade range
The second relativistic scaling relation that we consider, which is valid instead in the direct-cascade range, reads [1] ω Note that f k is a random external force. In the incompressible non-relativistic limit, ε becomes proportional to the enstrophy dissipation rate ω [35], namely ε → ((1 + w)/w) 2 ρ 2 ω . Expanding the left-hand side of Eq. (2.8) and setting ρ = constant and γ = 1 as before, we obtain The second term on the right-hand side is proportional to the average vorticity, which vanishes by parity invariance. Thus Eq. (2.8) becomes where the non-relativistic vorticity is mn ∂ m v n ≡ ω NR . The left-hand side needs to be manipulated further in order to compare with standard results (e.g. [35]). First, notice that the ensemble average on the left-hand side expands under the product rule to We can show that the second term on the right-hand side is zero as follows: where we used the identity ik mn = δ im δ kn − δ in δ km in the first line. Again, isotropy and regularity at the origin will imply this vanishes, provided it is divergence-free [36]. Thus, we can compute its divergence and show that it vanishes: where we used incompressibility in the second and last lines. The relativistic scaling relation Eq. (2.8) thus reduces in the incompressible Navier-Stokes limit to (2.14) As before, this relation is equivalent to an intermediate standard result from [35], namely To see this, expand the left-hand side and use statistical symmetries to obtain and then note that the second term on the right-hand side vanishes by incompressibility, isotropy, and regularity at r = 0 [36]. Thus Eq. (2.15) is the same as Eq. (2.14). Since Eq. (2.14) is the incompressible Navier-Stokes limit of Eq. (2.8), it can be compared in the direct-cascade range of relativistic or compressible turbulence in order to gauge their relative performance. Finally, we demonstrate that the relativistic correlation derived in [1], which reads also reduces to a known result in the incompressible non-relativistic limit. Note that the subscripts (L, T ) refer to the longitudinal ( r) and transverse (⊥ r) directions, respectively. Once again, setting ρ = constant and γ = 1 yields Using statistical symmetries, the left-hand side expands to 4 v i v i v j + 2 v j v i v i , and the second term vanishes due to incompressibility, isotropy, and regularity at the origin [36]. Thus, setting j = L in Eq. (2.19), we obtain We can eliminate the first term on the left-hand side using the well-known 1/8-law, also derived in [35] and valid in the direct-cascade range, (δv L ) 3 = 6 v L v L v L = (1/8) ω r 3 . This substitution finally yields Eq. (2.18).

III. IMPLEMENTATION
As stated, our goal is to explore scaling relations in conformal fluid turbulence. To ensure a clean inertial regime is established to compute the appropriate quantities, we include a driving source. Additionally, we ensure the numerical methods employed are consistent with the statistical properties of the flow we want to study. In this section we describe key aspects of our numerical implementation, beginning with general considerations in Sec. (III A). Following this, we present specific considerations for the incompressible and relativistic cases in Secs. (III B) and (III C), respectively.
A. General considerations

Stochastic Runge-Kutta
In order to implement a random white noise force in a simulation, a special integration algorithm must be used. Based on the work of Honeycutt [37], we use a second-order Stochastic Runge-Kutta algorithm (SRKII). The Gaussian random force we use, defined later in Eq. (3.2), is homogeneous, which means the average and variance of the force at every point in space is the same. Thus the prescription described in [37] is applied to each real space point, producing control over the injection rates in an aggregate sense.

Pseudorandom number generation
The random force we employ requires pseudorandom number generation at every time step. For this purpose, we implement the Intel MKL Vector Statistical Library. In particular, we use the Mersenne Twister 1 [38] and blocksplitting for parallel applications [39]. We have checked that the energy spectrum E(k) in steady-state is unaffected by the choice of random number generator by comparing the Mersenne Twister (VSL_BRNG_MT19937) and the 59-bit multiplicative congruential generator (VSL_BRNG_MCG59). We also checked that the output of our code is system-independent [39] by running it on two independent clusters.

Defining an injection length scale
In studies of turbulence, the energy/enstrophy injection and scale play a crucial role in establishing and identifying particularly relevant dynamical ranges. One can define an injection length scale associated with the external force in terms of the injection rates of energy and enstrophy as follows. Given Kraichnan-Batchelor [40] scaling of the energy spectrum in the inverse and direct cascades, respectively, one can take the injection scale to be the wavenumber at which E(k) transitions between these two scalings. Thus, set f and solve to find k f = η 0 / 0 . This definition will accurately represent the injection scale up to a numerical factor of order ∼ 1, so long as the energy spectrum transitions between these two behaviours over a short range of wavenumbers.

Formulation
In the incompressible Navier-Stokes case in 2D, the entire dynamics is determined by a single pseudo-scalar quantity, the vorticity ω = ∇ × v. Thus, it is computationally more efficient to evolve the vorticity equation directly, rather than the components of the velocity. We write the vorticity equation in "flux-conservative form", where f ω is the random force defined in the next section, and the dissipative term −ν 4 ∂ 4 ω ≡ −ν 4 ∇ 4 ω on the righthand side is often referred in the turbulence literature as "hyperviscosity of order 4". Hyperviscosity is frequently used in simulations of an incompressible Navier-Stokes fluid [20], since it limits the range of scales over which dissipation is active (yielding wider inertial ranges for a given grid resolution).

Random force and injection rates
The external force appears as f ω ≡ ∇ × f , and we wish to construct f ω directly with the appropriate statistical properties. Given a Gaussian random force with a two-point correlation in real space given by for some function g(r), the injection rate of enstrophy will be given by g(0)/2 ≡ η 0 [41], owing to the delta function (i.e. white noise) and to the choice of Gaussian randomness. Ignoring the temporal part of the correlation, we have in Fourier space where reality of the force in real space requires f ω (−k) = f * ω (k). In order to specify the enstrophy injection rate η 0 , we use a rescaling strategy as follows. First, define two random scalar fields A(k), B(k), with zero average A = B = 0 and unit variance A 2 = B 2 = 1 at all wavenumbers, and setf ω (k) = A(k) + iB(k). We first seek an isotropic rescalingf ω →g(k)f ω that gives the profile of Eq. (3.3) up to a constant factor. Under this rescaling, A, B →gA,gB, so the zero average is unchanged but the variance transforms to A 2 , B 2 →g 2 A 2 ,g 2 B 2 =g 2 . Thus, Thus choosingg ∝ ĝ/2 gives the desired spatial profile up to a constant factor. To fix the enstrophy injection rate (as η 0 = g(0)/2), we seek a second rescalingf ω → Rf ω with R = constant determined as follows. As it stands, Eq. (3.4) will produce an enstrophy injection rate given by half of its inverse Fourier transform evaluated at r = 0, Under the second rescaling, Eq. (3.4) becomes 2R 2g2 (k). Thus the appropriate rescaling is If one wishes instead to specify the energy injection rate, simply note that for a solenoidal force ∇ · f = 0, we have the spatial part of Eq. (3.2) given by So by solving the Poisson equation ∇ 2 f (0) · f (r) = −g(r) one finds the energy injection rate 0 from the relation The rescaling factor R can be chosen appropriately in this case. Extracting these a priori injection rates of energy and enstrophy allows one to define an injection length scale as per Sec. (III A 3). For our incompressible simulations of the direct-cascade we use a 'rectangular' profile, namelyĝ(k) = 1 in a narrow range of wavenumbers around k f , zero otherwise.

Dealiasing
The Navier-Stokes equation has a quadratic nonlinearity. Thus, two wavenumbers k 1 , k 2 can interact to populate a third wavenumber k 3 = k 1 + k 2 . Since we have a finite range of scales resolved in any simulation, k 3 could exceed the largest resolved wavenumber, and thus would become represented on the grid as a lower wavenumber N − k 3 (where N is the grid resolution). In this case, we say k 3 has been aliased. Prescriptions exist to avoid such aliasing errors. For a quadratically nonlinear term F × G, if we filter out all wavenumber modes with k > N /3 in F and G prior to multiplication, then filter F × G in the same manner, we will eliminate all aliasing errors. Such a prescription is known as the 2/3-dealiasing rule, since one retains 2/3 of the domain in Fourier space. Analogous dealiasing rules exist for higher-order nonlinearities, with less and less of the domain being retained as the order increases. Thus, full dealiasing becomes computationally prohibitive for higher-order nonlinearities, such as for a relativistic fluid flow.
C. Relativistic conformal fluid case

Formulation
The system of equations is given by ∇ a T ab = f b and the conformal perfect fluid stress-energy tensor T ab = (3/2)ρu a u b + (1/2)ρη ab , which uses the conformal equation of state P = ρ/2 in (2 + 1) dimensions. Defining the conservative variables as (D, S i ) = (T 00 , T 0i ), they appear in terms of the primitive variables as where v i is the spatial velocity and γ is the Lorentz factor. In terms of these variables, the equations of motion appear in flux-conservative form as We use finite differences to discretize the derivatives, with RK4 in space and SRKII (see Sec. (III A 1)) in time. The system is damped at short wavelengths using a 4th-order dissipation scheme discussed in Sec. (III C 3).

Random force and injection rates
We choose the Gaussian white-noise force f i to be divergence-free by deriving it from a stream function ψ, (f x , f y ) = (∂ y ψ, −∂ x ψ). Thus, numerically we build ψ directly in the manner described in Sec. (III B 2). For simulations of the inverse-cascade, we choose where l f is the characteristic length scale of the correlation, and = T 0i f i [1] is a constant. One can verify the equality = T 0i f i by applying the 2-dimensional Laplacian to Eq. (3.11), then noting that the spatial part of f i (r)f i (0) , written as F i i ≡ trF , is given by trF = −∇ 2 ψ(r)ψ(0) and trF = 2 T 0i f i [1]. In the weakly compressible regime, is approximately the injection rate of (1/2) T 0i T 0 i , whereas in the incompressible regime it fixes the Newtonian kinetic energy injection rate.
For simulations of the direct-cascade, we instead choose

Dealiasing
As alluded to in Sec. (III B 3), in the relativistic case a full dealiasing is computationally prohibitive. Since the computation of the velocity from the conservative hydrodynamic variables, followed by the computation of the flux, amounts to forming a product of up to 5 fields, there is a quintic nonlinearity. In the weakly-compressible regime, however, a 2/3-dealiasing rule would likely eliminate a satisfactory amount of aliasing, since the density and Lorentz factor have a small amount of power at all wavenumbers k = 0. However, in a future study we wish to explore the strongly compressible and ultrarelativistic regimes where a 2/3-rule would be inadequate. Thus we opt instead to use a 4th-order numerical dissipation scheme to suppress large wavenumber modes (since we want to explore the suitability of alternative dealiasing strategies for that future study) and employ a sufficiently high resolution (so that possibly spurious effects stay mainly confined at very high frequencies). For a variable U , this scheme amounts to including a term −ν num (∂ 4 x + ∂ 4 y )U on the right-hand side of its evolution equation, where ν num > 0 is the strength of the dissipation. It is numerically convenient to write this term as −κ(dx 3 ∂ 4 x + dy 3 ∂ 4 y )U and control the dissipation strength κ, as its magnitude will be closer to 1 and the dissipation length scale will move with the resolution [42].

IV. RESULTS
In all simulations we use periodic boundary conditions with a box size of L = 2π and resolution of N 2 = 2048 2 , with a variable step size determined by a CFL condition. This resolution has proven quite adequate for studying correlation functions in both the inverse-cascade (eg. [43]) and direct-cascade (eg. [44,45]) in incompressible fluid turbulence. We find it is also adequate for the weakly compressible regime studied here.
The time scale over which a turbulent flow is presumed to erase knowledge of its initial conditions is the large-eddy turnover time, which has various interpretations in the literature. Borue [46] estimates it as T = 2π/ω rms , where ω rms is the root-mean-squared vorticity. More generally, we have T = L/U where L is the scale of the largest eddies and U is a characteristic speed at that scale. L is estimated as 2π/k i , where k i is the infrared "cutoff" (∼ largest energy-containing scale), and we estimate U as the root-mean-square of the velocity. In our simulations, these time scales will be quoted for reference.
Averages will be computed over time, or over independent simulations, or both. The adequacy of the sample sizes is gauged via comparison of the average with the statistical error σ/ √ N , where σ is the sample standard deviation and N is the sample size. For example, a correlation function f (r) will have an ensemble of values for each r, and σ(r) is computed as the standard deviation of that collection of values.

A. Inverse-cascade simulations
We simulate the inverse-cascade of a (2 + 1)-dimensional conformal fluid with an external force described by Eq. (3.11), and an injection scale k f ≡ 2π/L f ∼ 203 defined by k f = η 0 / 0 , as in Sec. (III A 3). We consider three cases with the numerical dissipation strength given by κ = (0.05, 0.03, 0.02) (so as to compare results among them) and when quoting properties of each case we will present them in this order. Since the force is somewhat broadband, it has power in the dissipation range of scales. Thus, decreasing the dissipation strength is enough to increase the energy growth rate, and thus the rms Mach number of the flow, v rms /c s , where c s is the sound speed (1/ √ 2 of the speed of light, in our case). Statistical quantities are averaged over ensembles of independent simulations, as well as averaged over an interval of time after the energy passes k = 10 and before it reaches the box size. Table (I) contains various parameters of the flows, as well as the sample sizes for the joint average over an ensemble and over time.
In Fig. (1) we characterize the flows by presenting the probability distributions functions (pdfs) of the energy density and Mach number. The pdfs are observed to widen as the energy growth rate increases, as one would expect. For comparison, in both cases we also plot Gaussian distributions (black, dashed) with average and standard deviation matched to the data from the κ = 0.02 case. The Gaussian provides a good fit to the Mach number pdf (although with a slight hint of non-Gaussianity in the tail), whereas the energy density pdf exhibits a stronger, exponential tail towards smaller values. In Fig. (2) (Left) we display the angle-averaged Newtonian kinetic energy spectra E(k) ≡ π v 2 (k) (both the full spectrum and the potential part, obtained by projecting the velocity ontok in Fourier space). We observe a steepening of the inertial range scaling towards E(k) ∼ k −2 , which we note is steeper than the Kolmogorov/Kraichnan power law of k −5/3 . The spectra are not changed significantly (< 1%) by instead using density-weighted velocities ρ 1/3 v or ρv, the former having been suggested in the (3 + 1)-dimensional context in [47] to restore Kolmogorov/Kraichnan scaling from the observed spectral exponent of k −2 . The spectrum of the potential component of the velocity exhibits a bump beginning at k ∼ 30, with scaling of k −2.2 and k 0.78 on either side. Such a bump towards large k is commonly observed in spectra in simulated compressible flows in (3 + 1) dimensions, eg. [19,24,47,48], and is attributed in those cases to an artefact of high-order numerical dissipation known as the bottleneck effect [49]. This effect has also been observed in simulated compressible 2D flows which exhibit transfer of energy to small scales [50]. The bump we observe in Fig. (2) is likely due to the same effect, although we cannot make a conclusive statement since we have not performed the specific resolution studies necessary to do so, nor have we used dissipation of a different order. The late-time spectra obtained from our simulations of the direct-cascade (not shown) also exhibit such a bump, and in that case we note that reducing the time step by half does not change the bump perceptibly.
With regard to the full inverse-cascade spectra in Fig. (2) (Left), it is worth noting that there is no large-scale friction. In [51], it was shown that the presence of large-scale friction can affect the inertial range spectrum in the incompressible Navier-Stokes case. In the same study it was also shown that measurements of the inertial range spectrum are not reliable without a sufficiently resolved enstrophy cascade (k max /k f ∼ 16, where k max is defined as N /3). We do not have the direct-cascade range resolved to this degree in Fig. (2) (k max /k f ∼ 3.4). The approach of the full spectrum towards k −2 is generally expected for compressible turbulence in both (3 + 1) dimensions (see eg. [48]) and (2 + 1) dimensions (see eg. [52]), although usually for much larger Mach numbers than our current simulations. With that said, (2 + 1)-dimensional conformal fluids are special (eg. having a very large sound speed and no mass density), and its turbulent regime is seldom studied (see eg. [13,25]), so one may not expect the same TABLE I. Parameters of inverse-cascade simulations: κ is the dissipation parameter described in Sec. (III C 3); is the growth rate of (1/2) T0iT i 0 ; vrms/cs is the rms Mach number; Nt is the number of snapshots averaged over time; Nens is the number of independent runs (ensemble size); 2π/ωrms is the eddy turnover time defined by the rms vorticity; L/vrms is the eddy turnover time defined by vrms and L = 2π/10; δT is the time interval between snapshots of the flow that are averaged over; T1 and T2 are respectively the first the last times over which the temporal average is computed. For comparison, note that the light-crossing time is 2π. energy spectra a priori. We elaborate more on this in Sec. (V), where we demonstrate that the k −2 spectrum is not necessarily associated with compressive effects. In Fig. (3) we plot the relativistic correlation function appearing in Eq. (2.1), T 0i T i L , compensated for the expected scaling r −1 , with a linear vertical scale to help distinguish different power laws. We also plot the incompressible limit of that correlation function, obtained by setting ρ = γ = 1 (herein "the incompressible correlation"), as well as a non-relativistic but compressible version obtained by setting γ = 1 only (herein "the compressible correlation"). The former is equivalent to known results from incompressible Navier-Stokes turbulence (see Sec. (II A 1)), while the latter can be obtained from the left-hand side of Eq. (2.1) using the non-relativistic perfect fluid energy-momentum tensor, which is just the relativistic one with γ = 1. We use these comparisons to separately gauge the degree to which compressive and relativistic effects are important. In addition, we also include the predictions for each case, in matching colour, obtained from Eq. (2.1) and evaluations at γ = 1 and γ = ρ = 1 thereof. Error bars correspond to the statistical uncertainty σ/ √ N . For ease of comparison across cases, each plot has the same vertical axis range. As it is clear from the figure, we observe a progressive degradation of the scaling of the incompressible and compressible correlation functions as dissipation is weakened (and thus Mach number grows), while the relativistic one predicted in [1] outperforms. This is shown quantitatively in Fig. (4), where we display power law fits performed over the shaded interval of Fig. (3). The shaded interval is the same across all cases in order to make a fair comparison, and is chosen to capture the power law observed in the κ = 0.02 case (which has the narrowest scaling range). As dissipation is weakened, a monotonic shallowing of the best-fit power law is observed for both the incompressible and compressible correlation functions. This trend is more significant for the incompressible correlation function. The absolute performance of the relativistic κ =0.05 3. The relativistic correlation function T 0i T i L (solid blue) and its non-relativistic compressible and incompressible counterparts, T 0i T i L |γ=1 (dashed green) and T 0i T i L |γ=ρ=1 (dotted red), respectively, compensated by r −1 . From left to right: cases with dissipation strength κ = 0.02, 0.03, 0.05, respectively. Each prediction for the inverse-cascade range r/L f ∼ (10 0 , 10 1 ) is plotted as a horizontal line with matching line style. The predictions follow from Eq. (2.1) and evaluations at γ = 1 or γ = ρ = 1 thereof. Note that the centre and right plots have a nearly indistinguishable prediction for the relativistic and incompressible correlation functions (solid blue and red dotted lines, respectively). Error bars correspond to the statistical uncertainty σ/ √ N for each value of r/L f , where N is the sample size and σ is the sample standard deviation. The shaded grey area indicates the range of r/L f over which we fit a power-law, and we use the same range across all cases to ensure a fair comparison. Note the linear vertical scale, which accentuates deviations from the expected power law.
correlation function is superior to the compressible and incompressible correlation functions across all cases, and its relative performance improves as dissipation is decreased (i.e. differences in best-fit power law become larger).
The relativistic correlation function, although exhibiting power-law scaling ∼ r in all cases, nonetheless exhibits an increasing disagreement with the magnitude of the prediction in Eq. (2.1) (see Fig. (3)). In the most extreme case (κ = 0.02, v rms /c s = 0.19), the overall magnitude is less than the prediction by ∼ 4%. Our numerical ensemble of flows may be biased towards lower magnitudes, since runs with sufficiently large fluctuations from the random force can become numerically unstable and fail. As dissipation is weakened, this occurs more often. Thus, it is possible that the increasing disagreement of the magnitude of T 0i T i L and (1/2) r is an artefact of this bias. A high-resolution shock-capturing implementation could determine whether this increasing disagreement is a real effect. Power-law scalings 4. Least-squares power law fits for the relativistic correlation function T 0i T i L (circles) and its non-relativistic compressible and incompressible counterparts, T 0i T i L |γ=1 (triangles) and T 0i T i L |γ=ρ=1 (squares), respectively. All three dissipation cases are displayed. Error bars correspond to the standard deviation of the fitted power law scaling, obtained via random resampling with replacement (10 3 trials). The relativistic correlation function outperforms its compressible and incompressible counterparts across all cases, with a monotonic degradation observed for the latter two as the dissipation strength is decreased (and correspondingly, as the rms Mach number is increased).

B. Direct-cascade simulation
To simulate the direct-cascade of a (2+1)-dimensional conformal fluid, we instead use an external force with support only around k f = 7, as described by Eq. (3.12). As in our inverse-cascade simulations, we use 4th-order numerical dissipation as in Sec. (III B 3), with the choice κ = 0.01. However, in contrast to our inverse-cascade simulations, here we use a large-scale dissipation mechanism known as 4th-order hypofriction, which takes the form of a term −µ∇ −4 S i on the right-hand side of Eq. (3.10). We compute the inverse Laplacians spectrally, setting constant modes to zero. Such a term has power restricted to large scales, and terminates the brief inverse cascade from k = 7 towards k = 0. We find the value µ = 0.15 to be adequate for preventing a build-up of energy (and eventual condensation) at large scales. An energy condensate would be characterized by continued energy growth and the emergence of two dominant vorticies of opposing parity superposed on a noisy flow (see eg. [53]). Statistical quantities are averaged over the shaded interval of time indicated in Fig. (5) (Left), which consists of N = 10 snapshots separated by δT ∼ 0.68. The interval is chosen to maximize the number of snapshots available (to minimize statistical fluctuations), while remaining in a regime that roughly resembles a steady-state. Note that the measured correlations are not significantly affected if the temporal average begins slightly earlier or later. The average time step over this interval is 10 −3 t eddy , where t eddy = 2π/ω rms ∼ 0.5 is also averaged over the shaded interval. The injection rate of (1/2) T 0i T i 0 is measured initially to be = 2 × 10 −4 . For reference, the characteristic time as per v rms is L/v rms ∼ 26, where we take L = 2π/3 since the maximum of the energy spectrum occurs at k = 3. The rms Mach number over the shaded interval is ∼ 0.11, once again indicating the weakly compressible regime. In Fig. (5) (Left), we display the average Newtonian specific kinetic energy of the fluid as a function of time. As mentioned, we average various quantities over the shaded interval of time. The energy is beginning to plateau over this interval, however it continues to grow slowly. If evolved longer, the compressive component of the velocity begins to dominate over the curl-free part. To study such a regime more accurately, a Riemann solver would be desirable in order to more faithfully capture the dominant shockwave phenomena. Since we are instead using artificial high-order numerical dissipation, we choose to restrict our analysis to earlier times, when the compressive component of the velocity field is still subdominant (∼ 10% of the total energy at a given scale k -see Fig. (6) (Left)). Our high-order dissipation also results in large bottleneck effects at later times, which contaminate a rather large portion of the inertial range.
In Fig. (5) (Centre and Right), we display the probability distributions of the energy density (Centre, blue, solid) and Mach number (Right, blue, solid). For comparison, in both cases we also plot Gaussian distributions (black, dashed) with average and standard deviation matched to the data. Similar to the inverse-cascade simulations, the Gaussian provides a good fit to the Mach number pdf (although with weaker hints of a non-Gaussian tail in this case), whereas the energy density pdf exhibits a stronger, exponential tail towards smaller values.
In Fig. (6) we display the power spectra of the velocity (Left) and energy density (Right). The full energy spectrum of the flow (blue, solid), together with the energy spectrum of the compressive, curl-free, potential part of the velocity (cyan, dashed). The latter is seen to be subdominant by a factor of ∼ 10 over the range k ∈ [10, 100], which qualitatively corresponds to the direct-cascade interial range. The potential spectrum is fit by a power law k −3.13 over this range, while for the full spectrum we observe k −3 scaling with the multiplicative logarithmic correction ln(k/k f ) −1/3 . The inset shows the full spectrum compensated by k 3 with and without the logarithmic correction, with the presence of the logarithmic correction being favoured (flatter curve). In the literature, the presence of this correction seems to depend on several factors, including the length of time over which the average is taken, and the presence of large-scale dissipation [44,45,54,55].
In Fig. (7), we display the two measured correlation functions for which we have predictions in the direct cascade (Eqs. (2.8) and (2.17)). Errors again correspond to the statistical uncertainty σ/ √ N . We find reasonable agreement in the case of Eq. (2.8) (Left), and less so in the case of Eq. (2.17) (Right). The measured power laws are r 0.84 (Left) and r 2.4 , as compared to the predictions of r and r 3 , respectively. However, we note that the non-trivial factors of 1/2 (Left) and 1/24 (Right) yield marked agreement in magnitude. We suspect that by decreasing contamination from our modified large-and small-scale dissipation mechanisms, agreement with our predictions would improve, since in our inverse-cascade simulations we found that removing large-scale dissipation altogether improved agreement with our predictions significantly. We also note that the statistical uncertainty at short length scales is large enough that the sign of the correlation functions is uncertain there. We estimate the quantity ε by its incompressible limit ((1 + w)/w) 2 ρ 2 ω , with a further substitution of ρ 2 = 0.96 in place of ρ 2 (which improves agreement slightly). This estimate is justified by the fact that the correlation functions we measure do not differ significantly from their incompressible counterparts (i.e. setting ρ = γ = 1) in this regime. is fit by ∼ r 1.05 at shorter distances and ∼ r 2.40 at longer distances, with the transition occurring near r/L f ∼ 10 −1 . Error bars correspond to the statistical uncertainty σ/ √ N for each value of r/L f , where N is the sample size and σ is the sample standard deviation. Note that the statistical uncertainty at short distances is sufficiently large that the sign of the correlation functions is uncertain there. The corresponding compressible and incompressible limits of these correlation functions (obtained by setting γ = 1 or γ = ρ = 1, as in Sec. (IV A)) do not behave significantly differently. This allows us to approximate ε by its incompressible limit ((1 + w)/w) 2 ρ 2 ω as in Sec. (II A 2), and we take ω to be the initial enstrophy growth rate (before dissipation mechanisms become important). However, we substitute ρ 2 ∼ 0.96 in place of ρ 2 , which slightly improves agreement with the predictions.

V. DISCUSSION
As mentioned in Sec. (IV A), we observed that our energy spectra approach a k −2 scaling in the inverse-cascade range. We point out that it was observed in [51] that the incompressible case exhibits the same scaling provided largescale friction is absent and the direct cascade is sufficiently resolved (k max /k f ≥ 16, where we define k max = N /3). It thus becomes a prescient question whether a sufficiently resolved direct cascade in the conformal fluid case will yield the same result. To answer this, we perform an ensemble of 20 simulations with N = 2048 and k max /k f = 16, with a forcing profile given by Eq. (3.12). The resulting energy spectrum is displayed in Fig. (8), with the inset displaying the same spectrum compensated by either k 2 or k 5/3 . After filtering out the modes k ∈ [0, 5] (i.e. all modes less than the maximum of the spectrum), the rms Mach number for this flow is ∼ 0.11. We perform this filtering so as to have a more fair comparison of rms Mach number with our inverse-cascade simulations in Sec. (IV A), which we remind were ∼ 0.14, 0.17, and 0.19. As evident from Fig. (8), the spectrum clearly favours a k −2 description in the inverse-cascade range, and not the Kolmogorov/Kraichnan k −5/3 description. The best-fit power law over the range k ∈ [5,40] yields k −2.047 .
Interestingly, we note in passing that this k −2 scaling of the energy spectrum would change the result of the purported calculation of the fractal dimension of a turbulent (3 + 1)-dimensional AdS-black brane presented in [56] to D = 3 (rather than D = 3 + 1/3). An analysis of this will be reported elsewhere [57].
We also point out that, despite the narrow inverse-cascade range, we nonetheless observe a similarly narrow ∼ r 0.95 power law scaling in the same correlation functions analyzed in Sec. (IV A) (not shown). This suggests that the scaling relation Eq. (2.1) continues to hold with a more resolved direct-cascade range.
In Sec. (IV B), we observed hints of the predicted scaling of the correlation functions displayed in Fig. (7). By instead simulating an incompressible fluid (as per Sec. (III B)), for which shockwave phenomena are not present, we can measure the incompressible limits of Eqs. (2.14) and (2.18) with greater statistical significance. Conformal fluids have been shown to possess a scaling limit to an incompressible Navier-Stokes fluid [58,59]. We present the results of our simulation in Fig. (9) (Right). The enstrophy injection rate is set a priori to ω = 28. We use 4th-order Power laws k −2 and k −5/3 (black, solid) are shown for comparison. The inset displays compensated spectra E(k)k 2 (green, solid) and E(k)k 5/3 (red, dashed). The spectrum is well-represented by k −2 , rather than k −5/3 , consistent with [51] (albeit for a conformal fluid in our case). The best fit power-law slope over the range k ∈ [5,40] is −2.047.
hypofriction −λ∇ −4 ω and hyperviscosity −ν 4 ∇ 4 ω, with λ = 0.15 and ν 4 = 10 −10 . The forcing profile is given by and the injection scale is set to k f = 7. The average specific kinetic energy is plotted as a function of time in Fig. (9) (Left), with the averaging interval shaded gray. The energy spectrum compensated by k 3 and k 3 ln (k/k f ) 1/3 is displayed in Fig. (9) (Centre), with the shaded envelopes indicating 5× the statistical uncertainty σ/ √ N . The logarithmic correction is clearly favoured. In Fig. (9) (Right), we plot the correlation functions v T v L v T and − ω NR ω NR v L , together with their respective predictions in the direct cascade (1/24) ω r 3 and (1/2) ω r. The agreement is very much improved over Fig. (7), which suggests that despite the low Mach number in our direct cascade simulation of the conformal fluid, that situation is nonetheless quite different from the incompressible case (at least insofar as the numerical challenges are greater in the former case, eg. small scales being contaminated by bottleneck effects).
Finally, in Fig. (10) we display snapshots of the vorticity from several of our simulations. By doing so, we intend to provide intuition as to how the inverse and direct cascades appear in real space. In particular, the stretching and mixing of vorticity isolines characteristic of the direct-cascade range are readily identified as 'turbulence' qualitatively, where coherent features are seen over a variety of scales (see Fig. (10) (Bottom Left) and (Bottom Right)). By contrast, the inverse-cascade range has a much noisier appearance, as in Fig. (10) (Top Left). We have not observed an explicit acknowledgement of this qualitative fact in the literature, since numerical studies of the inverse-cascade range seldom include plots of the vorticity (eg. [43]). Unless the direct-cascade range is resolved, the turbulent flow qualitatively appears as random noise -but even if it is resolved, a clear hierarchy of scales is not apparent in the inverse-cascade range. In Fig. (10) (Top Right) we show a mixed case with the forcing acting at an intermediate scale k f = 42. This case is a simulation targeting the inverse-cascade range, but the direct-cascade range is just beginning to be resolved as well. Consequently, vorticity isoline mixing is beginning to be apparent, superposed on top of a more noisy structure.