Anton Svirsky

Corentin Herbert

Anna Frishman

Statistics of inhomogeneous turbulence in large scale quasi-geostrophic dynamics

A remarkable feature of two-dimensional turbulence is the transfer of energy from small to large scales. This process can result in the self-organization of the flow into large, coherent structures due to energy condensation at the largest scales. We investigate the formation of this condensate in a quasi-geostropic flow in the limit of small Rossby deformation radius, namely the large scale quasi-geostrophic model. In this model potential energy is transferred up-scale while kinetic energy is transferred down-scale in a direct cascade. We focus on a jet mean flow and carry out a thorough investigation of the second order statistics for this flow, combining a quasi-linear analytical approach with direct numerical simulations. We show that the quasi-linear approach applies in regions where jets are strong and is able to capture all second order correlators in that region, including those related to the kinetic energy. This is a consequence of the blocking of the direct cascade by the mean flow in jet regions, suppressing fluctuation-fluctuation interactions. The suppression of the direct cascade is demonstrated using a local coarse-graining approach allowing to measure space dependent inter-scale kinetic energy fluxes, which we show are concentrated in between jets in our simulations. We comment on the possibility of a similar direct cascade arrest in other two-dimensional flows, arguing that it is a special feature of flows in which the fluid element interactions are local in space.

I. INTRODUCTION

Multi-scale, non-linear, interactions are one of the defining properties of turbulent flows, posing a considerable challenge both for theoretical understanding and numerical modeling. In particular, they imply that the flow at large scales is coupled to smaller scale fluctuations, and that, for instance, the structure of the large scale flow depends on the transport of momentum and the dissipation of energy by such small scales. For statistically homogeneous and isotropic flows such interactions are, to leading order, well captured within phenomenological theories, though many theoretical questions remain open and there are very few results which can be obtained from first principles [1, 2, e.g.]. However, most real flows break such symmetries at large enough scales, either because of external fields such as gravity or a magnetic field, due to the effect of rotation, or the existence of boundaries. In such flows, often the key question is to characterize the large scale mean-flow. In turn, this requires the prediction of energy and momentum transfers due to turbulent fluctuations, across scales and spatially. Generally, this is a challenging task, requiring ad-hoc assumptions. However, there is growing evidence that the non-linear interactions in the presence of a strong mean flow may in fact be easier to treat than those in a homogeneous and isotropic flow (see e.g. [START_REF] Marston | Recent developments in theories of inhomogeneous and anisotropic turbulence[END_REF] for a review). This has been particularly evident in two-dimensional (2D) and quasi-2D flows, where dimensionality imposes strong constraints upon the nature of multi-scale interactions.

Two-dimensional flows exhibit a remarkable tendency to spontaneously self-organize into a coherent mean flow when excited at small scales. The mechanism behind * frishman@technion.ac.il this self-organization is an inverse transfer of a quadratic invariant (e.g energy) from small to large scales, in a process called the inverse cascade [START_REF] Kraichnan | Inertial Ranges in Two-Dimensional Turbulence[END_REF][START_REF] Leith | Diffusion Approximation to Inertial Energy Transfer in Isotropic Turbulence[END_REF][START_REF] Batchelor | Computation of the Energy Spectrum in Homogeneous Two-Dimensional Turbulence[END_REF]. The inverse cascade arises due to the existence of a second inviscid invariant of the dynamics, which is simultaneously transferred to small scales in a so-called direct cascade. In a finite system this inverse transfer results in the accumulation of energy at the largest available length scale, forming a system-size coherent mean flow termed a condensate [START_REF] Kraichnan | Inertial Ranges in Two-Dimensional Turbulence[END_REF][START_REF] Smith | Bose condensation and small-scale structure generation in a random force driven 2d turbulence[END_REF][START_REF] Chertkov | Dynamics of energy condensation in twodimensional turbulence[END_REF][START_REF] Chertkov | Universal velocity profile for coherent vortices in two-dimensional turbulence[END_REF]. In this condensate regime, the direct interactions between the mean-flow and turbulence can dominate over local-in-scale interactions. Indeed, theoretical ideas and numerical methods which do not explicitly resolve the fluctuation-fluctuation interactions have been shown to be applicable in this system (and its variants) [START_REF] Farrell | Structural Stability of Turbulent Jets[END_REF][START_REF] Farrell | Structure and Spacing of Jets in Barotropic Turbulence[END_REF][START_REF] Tobias | Direct Statistical Simulation of Out-of-Equilibrium Jets[END_REF][START_REF] Marston | Generalized Quasilinear Approximation: Application to Zonal Jets[END_REF][START_REF] Srinivasan | Zonostrophic Instability[END_REF][START_REF] Bouchet | Kinetic Theory of Jet Dynamics in the Stochastic Barotropic and 2D Navier-Stokes Equations[END_REF]. Moreover, analytical results describing both the spatial structure of the mean flow [START_REF] Laurie | Universal Profile of the Vortex Condensate in Two-Dimensional Turbulence[END_REF][START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF][START_REF] Frishman | The culmination of an inverse cascade: Mean flow and fluctuations[END_REF] and the turbulent kinetic energy density [START_REF] Frishman | Turbulence Statistics in a Two-Dimensional Vortex Condensate[END_REF], were successfully derived from first-principles in this regime.

The interest in two-dimensional flows is not limited to the theoretical understanding of turbulent interactions, as many flows in nature become effectively twodimensional and thus exhibit a similar phenomenology. This occurs when the fluid motion is constrained in one of the directions either because the fluid is contained within a thin layer, is stratified, or is rapidly rotating [START_REF] Xia | Spectrally condensed turbulence in thin layers[END_REF][START_REF] Boffetta | Two-dimensional turbulence[END_REF][START_REF] Rivera | The direct enstrophy cascade of two-dimensional soap film flows[END_REF][START_REF] Bouillaut | Experimental observation of the geostrophic turbulence regime of rapidly rotating convection[END_REF]. Astrophysical and large-scale geophysical flows often have one or more of these properties, with rotation playing a particularly important role in constraining the motion, called the geostrophic regime [START_REF] Young | Forward and inverse kinetic energy cascades in Jupiter's turbulent weather layer[END_REF][START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF]. A minimal model for the flow in this regime, capturing the main features of the large-scale dynamics and serving as an important theoretical tool, is the shallow water quasi-geostrophic equation (SWQG). In SWQG, there is a typical scale which determines the range of interactions between fluid elements, called the Rossby deformation ra-dius. The special cases within which the condensate had been studied in detail so far are two-dimensional incompressible Navier-Stokes (2DNS) with the beta effect (differential rotation) [START_REF] Farrell | Structure and Spacing of Jets in Barotropic Turbulence[END_REF][START_REF] Tobias | Direct Statistical Simulation of Out-of-Equilibrium Jets[END_REF][START_REF] Srinivasan | Zonostrophic Instability[END_REF][START_REF] Woillez | Theoretical prediction of Reynolds stresses and velocity profiles for barotropic turbulent jets[END_REF][START_REF] Woillez | Barotropic theory for the velocity profile of Jupiter's turbulent jets: An example for an exact turbulent closure[END_REF], or 2DNS without differential rotation [START_REF] Laurie | Universal Profile of the Vortex Condensate in Two-Dimensional Turbulence[END_REF][START_REF] Frishman | The culmination of an inverse cascade: Mean flow and fluctuations[END_REF][START_REF] Frishman | Turbulence Statistics in a Two-Dimensional Vortex Condensate[END_REF]. Both these cases capture dynamics at scales much smaller than the deformation radius, so that interactions span the entire domain.

Here we consider the opposite limit, where interactions are local-the so-called large-scale quasi-geostrophic (LQG) equation [START_REF] Larichev | Weakly decaying turbulence in an equivalent-barotropic fluid[END_REF]. This model captures the long-time dynamics at scales much larger than the deformation radius. It contains two inviscidly conserved quantities: the potential energy, transferred to large scales, and the kinetic energy, which cascades to small scales in a direct cascade. The main question we investigate is what type of condensate does this system support, and how does the locality of interactions affect the properties of turbulent fluctuations and their interaction with the mean-flow.

We begin by reviewing the derivation of the LQG equation from SWQG and discuss the conditions necessary for the emergence of a condensate in this system. Building upon our results for the mean flow of an LQG jettype condensate [START_REF] Svirsky | Statistics of inhomogeneous turbulence in large scale quasi-geostrophic dynamics[END_REF], we fully characterize the two-point second-order statistics, combining analytical derivations and results from direct numerical simulations. Turning to the direct cascade and the kinetic energy balance, we show that in the presence of a mean flow it includes a spatial flux of fluctuating kinetic energy. Such a flux is absent in 2D Navier-Stokes, and seems to be related to the locality of interactions in LQG. Using our analytical results for the second-order statistics, we show that this flux carries most of the kinetic energy away from regions of strong mean flow, effectively arresting the direct cascade there. We confirm the in-homogeneity of the direct cascade induced by the mean flow by examining the flux of energy between scales within a smooth filtering approach [START_REF] Eyink | Localness of energy cascade in hydrodynamic turbulence. I. Smooth coarse graining[END_REF]. In particular, we consider the local inter-scale flux of kinetic energy and potential energy for the LQG system. Measuring this flux in simulations, we demonstrate that the flux of kinetic energy to small scales is indeed locally suppressed in regions where the jets are strong, an effect so strong it is evident under short-time averaging.

This work also serves as the companion to the paper [START_REF] Svirsky | Statistics of inhomogeneous turbulence in large scale quasi-geostrophic dynamics[END_REF]. Here we provide detailed derivations and discussions of some of the results stated in [START_REF] Svirsky | Statistics of inhomogeneous turbulence in large scale quasi-geostrophic dynamics[END_REF] alongside new results and analysis not contained in [START_REF] Svirsky | Statistics of inhomogeneous turbulence in large scale quasi-geostrophic dynamics[END_REF].

II. FRAMEWORK

The LQG equation can be derived as a limit of the shallow water quasi-geostrophic equation which is given by [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF] ∂ t q+v•∇q = ∂ t q+J(ψ, q) = 0; q = ∇ 2 -L -2 d ψ, [START_REF] Frisch | Turbulence: The Legacy of A.N. Kolmogorov[END_REF] where q is the potential vorticity, ψ is the stream-fucntion which is related to the fluid velocity via v = ẑ × ∇ψ,

ω = ∇ 2 ψ = (∇ × v)
• ẑ is the vorticity, J(ψ, q) is the Jacobian operator defined as J(ψ, q) = ∂ x ψ∂ y q -∂ y ψ∂ x q = ϵ ij ∂ i ψ∂ j q with ϵ ij the 2D Levi-Civita symbol. This equation describes the dynamics of a rapidly-rotating homogeneous fluid layer, wherein the pressure gradient force due to fluctuations of the free surface is balanced by the Coriolis force (the so-called geostrophic balance). The stream function both determines the velocity and is proportional to the surface height perturbations of the fluid layer. The scale L d is called the Rossby deformation radius and sets the range of influence of a surface perturbation on its surroundings. When L d /L → ∞, where L is a characteristic scale for the domain size, surface perturbations have a long-range influence on the fluid, and their equilibration is fast compared to the rotation period, giving an incompressible 2D fluid. The opposite limit L d /L → 0, corresponds to a very rapidly rotating fluid, where the effect of surface perturbations is strictly local. In this limit, the long time dynamics are given by the LQG equation [START_REF] Larichev | Weakly decaying turbulence in an equivalent-barotropic fluid[END_REF]

∂ τ ψ + v ω • ∇ψ = ∂ τ ψ + J(ω, ψ) = f + α∇ 2 ψ -ν(-∇ 2 ) p ψ, (2)
with v ω = ẑ × ∇ω and we have included forcing f and dissipation, with α the friction coefficient (corresponding to linear drag on velocity v) and ν the (hyper) viscosity.

This advection equation is similar to 2DNS but with the roles of the vorticity and the stream-function reversed. Here the vorticity acts as the "effective stream function", and the stream-function is advected by an effective velocity v ω = ẑ × ∇ω. The integral invariants of (2) without forcing and dissipation are the kinetic energy

Z = 1 2 (∇ψ) 2 d 2 x = 1 2 |v| 2 d 2
x and all moments of ψ, in particular the potential energy E = 1 2 ψ 2 d 2 x. The existence of the two quadratic invariants results in the inverse cascade of E and a direct cascade of Z [START_REF] Smith | Turbulent diffusion in the geostrophic inverse cascade[END_REF].

Let us briefly comment on the conditions for the LQG limiting dynamics to be consistent, a more complete discussion can be found in the Appendix A. The SWQG equation is derived from the rotating shallow water equations in the limit of a small Rossby number, Ro = U/(ΩL) with U a typical velocity scale, and Ω the fluid rotation rate. The derivation also requires that Ro(L/L d ) 2 ∼ o(1) so that height perturbations are small compared to the mean fluid thickness. The limit L d /L → 0 which we take to obtain LQG is consistent with this assumption provided that L d /L ∼ Ro β with 0 < β < 1/2, under which condition LQG can be derived as a limit of SWQG . We remark that traditionally SWQG is derived assuming L d /L ∼ O(1) rather than L d /L ∼ Ro β . However, we show in the Appendix A that LQG can be directly derived from the rotating shallow water equations in the latter limit, that requires rescaling time by τ = t(L d /L) 2 ∝ tRo 2β and expanding the height field in powers of Ro n+2β (instead of Ro n like the velocity). Thus, the forced LQG equation should be able to capture the large scale dynamics of SWQG with a small but finite L d with a forcing scale which is larger than L d [START_REF] Smith | Turbulent diffusion in the geostrophic inverse cascade[END_REF]. Note that there is evidence that the inviscid equation eventually develops motions on smaller scales, so the LQG equation may become inadequate [START_REF] Burgess | Potential vorticity fronts and the late-time evolution of large-scale quasigeostrophic flows[END_REF].

We wish to explore the LQG system (2) in the condensation regime, where the potential energy condenses at the largest available scale. This requires that the rate of energy removal at the box scale is much slower than the transfer rate by the inverse cascade. The inverse cascade rate (eddy-turnover time) can be found using dimensional analysis, requiring that this rate depends only on the scale and the potential energy injection rate ϵ = ⟨ψf ⟩ [START_REF] Smith | Turbulent diffusion in the geostrophic inverse cascade[END_REF]. We have [ϵ] ∼ [ψ 2 ]/t ∼ l 8 /t 3 , so that the rate of the inverse cascade of E at scale l is τ E (l) ∼ ϵ -1/3 l 8/3 . Similarly, the rate of the direct cascade of Z at scale l is τ Z (l) ∼ η -1/3 l 2 , where η = ⟨∇ψ∇f ⟩ is the kinetic energy injection rate. Assuming the forcing acts in a narrow band of scales around l f , the injection rates can be simply related by η ≈ ϵ/l 2 f . Note that the eddy turnover time decreases with the scale much faster than in 2DNS, with a factor of l 2 between the two. This limits the available resolution for simulations and thus also the separation of scales between the forcing and box scale.

The dissipation rates due to the drag and viscous terms are τ α (l) ∼ α -1 l 2 and τ ν (l) ∼ ν -1 l 2p respectively. Assuming p ≥ 2 (integer), the former will serve as the large scale dissipation mechanism arresting the inverse cascade of E while the latter as the small-scale dissipation arresting the direct cascade of Z. For the potential energy to condense at the box scale, L, requires

δ ≡ τ E (L)/τ α (L) = α(L 2 /ϵ) 1/3 ≪ 1.
Note that δ grows with the scale, i.e. that the ratio between the nonlinear time-scale and the dissipative time scale grows with the length scale (like in 2DNS with linear friction, factors of l 2 appearing in both time-scales cancelling out to give the same ratio). Additionally, in order for a significant fraction of the (potential) energy to be transferred to large scales, we require that at the forcing scale the dissipation rate τ ν (l f ) is low compared to the non-linear transfer rate, resulting in the requirement Re

≡ τ ν (l f )/τ E (l f ) = l 2p-8/3 f ϵ 1/3 /ν ≫ 1.
The kinetic energy cascade is arrested at the Kolmogorov scale l ν where the inverse cascade rate and the viscous dissipation rate are comparable, resulting in l ν ∼ (ν 3 /ϵ) 1/(6p-8) . We can then estimate the potential energy dissipation at small scales: assuming a constant kinetic energy flux down to the Kolmogorov scale, the energy dissipation rate is given by ϵ ν = l 2 ν η where η = ϵ/l 2 f is the injected kinetic energy. Thus we get ϵ ν = (l ν /l f ) 2 ϵ. This implies that the ratio between the dissipated energy and the energy transferred to large scales is

ϵ ν /ϵ α = ϵ ν /(ϵ-ϵ ν ) = l 2 ν /(l 2 f -l 2 ν ) ≈ l 2 ν /l 2 f
which is indeed small in the limit of a large Re number.

The sharpness of the small-scale cutoff is determined by p, higher values will increase the kinetic energy removal rate at scales l < l ν (and decrease the removal rate for l > l ν ) resulting in a sharper cutoff of the spectrum at l ν . We perform direct numerical simulations (DNS) of the LQG equation (2) using the Dedalus framework [START_REF] Burns | Dedalus: A flexible framework for numerical simulations with spectral methods[END_REF]. The pseudo-spectral method is implemented using the 3/2 dealiasing rule and time stepping using a third-order, four-stage DIRK/ERK method. We focus on a jet-type LQG condensate, which simplifies the analysis in the following. In a doubly periodic domain such a condensate emerges if the symmetry between the x and y directions is broken [START_REF] Bouchet | Random Changes of Flow Topology in Two-Dimensional and Geophysical Turbulence[END_REF][START_REF] Frishman | Jets or vortices-What flows are generated by an inverse turbulent cascade?[END_REF]. Note that such jets differ from those which emerge due to differential rotation (beta-plane turbulence [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF]), which is absent here. We therefore use a doubly periodic box of dimensions L ≡ L y = 2L x = 2π. The spatial resolution is taken to be 64 × 128, which is relatively low, restricted by the rapid decrease of the eddy-turnover time with decreasing scale in LQG. We use a white-in-time forcing which is localized in Fourier space at a wavenumber k f = 2π/l f = (10, 13, 15) (forcing in an annulus of width 2dk = 2 with a constant amplitude A = 10 -3 and a random phase). We use hyperviscosity with p = 7 and ν = (0.5, 7.3, 10) × 10 -19 and take α = (0.5, 1, 2) × 10 -3 . Simulation parameters are chosen such that a significant fraction of the (potential) energy is transferred to large scales Re ≫ 1 and such that potential energy condenses at large scales δ ≪ 1. Each simulation is run until the system reaches a statistically steady state, and statistics are then gathered over many large scale turnover-times τ E (L). The full list of simulations performed is presented in Table I and the choice of the temporal and spatial resolutions are discussed in Appendix C.

The resulting condensate, with two alternating jets along the short side (x direction) of the domain, is shown in Fig. 1. Between the jets, there are two small vortices, similarly to what was found in 2DNS [START_REF] Frishman | Jets or vortices-What flows are generated by an inverse turbulent cascade?[END_REF], possibly due to instabilities of the mean flow. In the jet region, the flow is statistically homogeneous in x. Small magnitude oscillations along the y direction of the jet amplitude can also be seen in Fig. 1. In steady state, no significant drift of the profile is observed over time, so the mean profile is simply determined from the average of the snapshots without shift. After averaging, we set the axis such that the mean velocity is zero on the y = 0 line with U > 0 above it and U < 0 below it.

To obtain a statistical description of the steady state LQG jet condensate we decompose the flow into the mean Ψ = ⟨ψ⟩ and fluctuations ψ ′ = ψ -Ψ focusing on the jet region, where the flow is statistically homogeneous in x, and the mean flow depends on y only. The mean flow ∂ y Ψ ≡ -U (y) and the mass flux ⟨v ω ′ y ψ ′ ⟩ can be obtained from the mass flux balance (average of equation ( 2)) and the potential energy balance, neglecting kinetic energy dissipation for the fluctuations and cubic-in-fluctuations terms. The latter assumes that non-linear interactions are dominated by mean-flow-turbulence interactions at the relevant scales, and is also known as the quasi-linear approximation, see [START_REF] Marston | Recent developments in theories of inhomogeneous and anisotropic turbulence[END_REF] for a review. The derivation and the comparison to DNS are presented in [START_REF] Svirsky | Statistics of inhomogeneous turbulence in large scale quasi-geostrophic dynamics[END_REF], and here we only cite the resulting leading order solution

∂ y Ψ = -U = ± ϵ α , (3)
⟨ψ ′ v ω′ y ⟩ = ± √ ϵα. (4)
In agreement with (3), the simulated mean velocity U (y) is indeed constant in the region where the jets are strong, rapidly switching sign in a thin transition region between the jets, as can be seen in Fig. 1(b).

III. SECOND ORDER STATISTICS: TWO-POINT CORRELATION FUNCTIONS A. Analytical results

Given an expression for the mean flow, we are now in a position to go further in the perturbation theory and consider the full second order (single-time) statistics. It is sufficient to consider the two-point correlation function ⟨ψ ′ 1 ψ ′ 2 ⟩ ≡ ⟨ψ ′ (r 1 )ψ ′ (r 2 )⟩ where r i = (x i , y i ), from which other single and two-point second-order correlation functions can be subsequently derived. To obtain an expression for ⟨ψ ′ 1 ψ ′ 2 ⟩, we will use that in a statistically steady state

0 = ∂ τ ⟨ψ ′ 1 ψ ′ 2 ⟩ = i̸ =j ψ ′ j ∂ τ ψ ′ i . ( 5
)
The evolution equation for the fluctuations ∂ τ ψ ′ i = ∂ τ (ψ -Ψ) is obtained by subtracting the average of equation ( 2) from (2), giving

∂ τ ψ ′ = -v ω • ∇ψ + ∂ y v ω′ y ψ ′ + f + α∇ 2 ψ ′ -ν(-∇ 2 ) p ψ ′ , (6)
where we have used that V ω y ≡ ∂ x ⟨ω⟩ = 0 (due to homogeneity in x). Evaluating the derivative at point r i , multiplying by ψ ′ j and averaging gives

ψ ′ j ∂ τ ψ ′ i = -ψ ′ j v ω i • ∇ψ i + f i ψ ′ j + α ψ ′ j ∇ 2 ψ ′ i -ν ψ ′ j (-∇ 2 ) p ψ ′ i , (7)
where v ω i ≡ v ω (r i ). Note that no summation over i is implied here. The cubic term reads

ψ ′ j v ω i • ∇ψ i = ∂ y Ψ i ψ ′ j v ω′ y (r i ) + V ω x (r i )⟨ψ ′ j ∂ x ψ ′ i ⟩ + ψ ′ j v ω′ i • ∇ψ ′ i , (8)
again using that V ω y = 0. As the derivatives act on r i ̸ = r j we can take them out of the average, resulting in

ψ ′ j ∂ τ ψ ′ i = -∂ y Ψ i ∇ 2 i ∂ xi + V ω x (r i )∂ xi -α∇ 2 i + ν(-∇ 2 i ) p ψ ′ j ψ ′ i + f i ψ ′ j -∇ i • v ω′ i ψ ′ j ψ ′ i . (9)
where ∇ i and ∇ 2 i denotes the gradient and Laplacian with respect to r i . Finally using (5) and (9) we get i=1,2

∂ y Ψ i ∇ 2 i ∂ xi + V ω x (r i )∂ xi -α∇ 2 i + ν(-∇ 2 i ) p ⟨ψ ′ 1 ψ ′ 2 ⟩ =2χ 12 -∇ 1 • ⟨v ω′ 1 ψ ′ 1 ψ ′ 2 ⟩ -∇ 2 • ⟨v ω′ 2 ψ ′ 1 ψ ′ 2 ⟩ , (10)
where we have used that the force two-point correlation function is given by ⟨f (x 1 , t)f (x 2 , t)⟩ = 2χ 12 δ(t -t ′ ). Having derived the general equation for the two-point correlation function, equation ( 10), we will proceed using a perturbative approach. For the mean flow, we will use the leading order solution cited above. In the same manner as was done for the single point correlation function, we shall neglect the viscous dissipation term (there is no dissipative anomaly) as well as that by linear friction, since in the condensate regime we expect the dissipation of the fluctuations of ψ to be a sub-leading effect, e.g.

α∇ 2 i ⟨ψ ′ 1 ψ ′ 2 ⟩ ∂ y Ψ i ∇ 2 i ∂ xi ⟨ψ ′ 1 ψ ′ 2 ⟩ ∼ α 3/2 √ ϵ l 2 = δ 3/2 l 2 L ≤ δ 3/2 ≪ 1, (11)
where l 2 ≤ L is the length scale of the two-point function and assuming the condensate regime with δ ≪ 1. We will further use the quasi-linear approximation, expecting that at leading order the cubic fluctuation terms are negligible compared to the mean-flow-fluctuations term

∂ y Ψ i ∇ 2 i ∂ xi ⟨ψ ′ 1 ψ ′ 2 ⟩. This gives ∂ y1 Ψ 1 ∇ 2 1 ∂ x1 + ∂ y2 Ψ 2 ∇ 2 2 ∂ x2 ⟨ψ ′ 1 ψ ′ 2 ⟩ = 2χ 12 . (12)
Using homogeneity in the x direction (also for χ 12 which only depends on x 1 -x 2 ), i.e. that ∂ x1 = -∂ x2 when acting on the two-point function, and ∂ y Ψ = ϵ/α at leading order, simplifies the advection operator to

∇ 2 1 ∂ x1 + ∇ 2 2 ∂ x2 = (∂ 2 x1 + ∂ 2 y1 )∂ x1 -(∂ 2 x1 + ∂ 2 y2 )∂ x1 = (∂ 2 y1 -∂ 2 y2 )∂ x1 .
The equation then reads

(∂ y1 + ∂ y2 ) (∂ y1 -∂ y2 ) ∂ x1 ⟨ψ ′ 1 ψ ′ 2 ⟩ = 2 α ϵ χ 12 . ( 13
)
Changing variables to y + = (y 1 + y 2 )/2 and y -= (y 1y 2 )/2 = ∆y/2 so that

∂ y+ = ∂ y1 + ∂ y2 and ∂ y-= ∂ y1 - ∂ y2
, finally gives the equation for the two-point function in compact form

∂ y+ ∂ y-∂ x1 ⟨ψ ′ 1 ψ ′ 2 ⟩ = 2 α ϵ χ 12 . (14)
We now briefly outline the solution of equation ( 14), leaving the detailed derivation to Appendix B. The solution will be a sum of the particular and the homogeneous solutions of equation ( 14). We begin with the former, first noting that the forcing correlation function χ 12 in equation ( 14) should be replaced by χ12

χ12 = χ 12 - Ly 2 - Ly 2 ds L y χ 12 (∆x, s) - Lx 2 -Lx 2 ds L x χ 12 (s, ∆y). ( 15
)
where the ∆x and ∆y independent parts (the respective k x = 0 and k y = 0 Fourier modes) are subtracted. This is necessary as these modes do not satisfy the Fredholm alternative, so the particular solution for them must be determined at next order, see Appendix B. The modified equation can now be straightforwardly integrated to obtain the particular solution. Note that for small separations ∆x, ∆y ≪ l f χ12 ≈ χ 12 , while for ∆x, ∆y ≫ l f χ12 ≪ 1 so the influence of the forcing is limited to scales ∆x, ∆y < l f , see Appendix B.

While the forcing provides the leading order contribution to the odd in ∆x part of the correlation function, corresponding to parity+time reversal symmetry breaking, the even contribution at leading order must come from the homogeneous solutions to [START_REF] Srinivasan | Zonostrophic Instability[END_REF]. Those are the zero modes of the advection operator

L 1 + L 2 = ∇ 2 1 ∂ x1 + ∇ 2 2 ∂ x2 = ∂ y+ ∂ y-∂ x1 : ⟨ψ ′ 1 ψ ′ 2 ⟩ hom = C(∆y, ∆x) + C 1 (y + , ∆x) + C 2 (y + , ∆y). ( 16
)
The relevant form of the solution in our case is only C(∆y, ∆x) as we detail in Appendix B. Thus, the full solution reads

⟨ψ ′ 1 ψ ′ 2 ⟩ = C(∆y, ∆x) + 2y + α ϵ ∆x 0 dz ∆y/2 0 dz ′ χ12 (z, z ′ ) , (17)
Note that for the homogeneous part C(∆x, ∆y) = C(-∆x, ∆y) = C(-∆x, -∆y), where the first equality is a consequence of the invariance with respect to parity (x → -x) + time reversal (t → -t) (PT) which we expect the zero modes to have, and the second of the exchange symmetry r 1 → r 2 of the two-point correlation function. In addition, we get the prediction that for ∆x = 0, ∆y ̸ = 0 the correlation function is independent of y + , as confirmed in DNS Fig. 3. Our approach lacks information about the boundary conditions to be applied and treats the differential operator perturbatively. It is thus unclear how to determine C(∆x, ∆y), and it may require going to the next order in perturbation theory, which is beyond the scope of the present work. As a consistency check, we can compute the mass flux ⟨v ω′ y ψ ′ ⟩ directly from our result for the two-point function [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF]. In particular, we directly confirm that, being an odd correlator, it is determined by the inhomogeneous solution to the two-point function equation. We shall compute v ω′ y (r 1 )ψ ′ 2 and will subsequently merge the two points.

ϵ α v ω′ y (r 1 )ψ ′ 2 = = ∇ 2 1 ∂ x1 2y + ∆x 0 dz ∆y/2 0 dz ′ χ12 (z, z ′ ) + {odd} = ∂ 2 y1 2y + ∆y/2 0 dz ′ χ12 (∆x, z ′ ) + {odd} = χ12 (∆x, ∆y) + {odd} ( 18
)
where {odd} denotes terms odd in ∆y and ∆x which will vanish when we take the single point limit ∆x, ∆y → 0. Note that the zero mode indeed produces only odd contributions in ∆x (since C(∆x, ∆y) is even under x 1 → -x 1 while in [START_REF] Frishman | The culmination of an inverse cascade: Mean flow and fluctuations[END_REF] there is an odd derivative with respect to this variable), which do not contribute. Taking the limit ∆x, ∆y → 0, χ12 → ϵ, up to the contribution to the energy injection rate from modes with k x = 0 and k y = 0, assumed to be O(l f /L). We thus get the expected result of ⟨v ω′ y ψ ′ ⟩ = √ αϵ.

B. Simulation results for the two-point correlation function

We now present results from DNS for the two-point function ⟨ψ ′ 1 ψ ′ 2 ⟩(∆y, ∆x, y + ). Note that a-priori the correlation function also depends on x + (or x 1 ), but taking into account statistical homogeneity in the x direction, we also average over x + (in addition to time). Using the fact that both jet regions are statistically identical, we compute the two-point function for each of them and average the two to obtain better statistics. In Fig. 2(a) we present the variance ⟨ψ ′2 ⟩ as a function of y normalized by its value at y = 0 at the center of the jet. The jet region, where the leading order solution for the mean flow (3) applies, is defined by |∂ y U |/ ϵ/αL 2 < 1 and is delimited by dashed lines, and we expect that ⟨ψ ′2 ⟩ = C(0, 0) in this region. The large peaks in the variance ⟨ψ ′2 ⟩ outside the jet region are related to the vortices in between the jets. While they are coherent structures with a large amplitude, which we would normally associate with a mean flow, since they freely move across the domain they contribute to the fluctuations in our averaging procedure. In Fig. 2(b) we present the normalized correlation function ⟨ψ ′ 1 ψ ′ 2 ⟩(∆y, ∆x = 0, y + = 0)/⟨ψ ′2 ⟩(y = 0), showing how ψ correlations decay with ∆y when the separation between the points is taken symmetrically around a jets center. In Fig. 2(c) we present the normalized correlation function ⟨ψ ′ 1 ψ ′ 2 ⟩(∆y = 0, ∆x, y + = 0)/⟨ψ ′2 ⟩(y = 0), showing the correlations in the x direction for points at the center of the jet y 1 = y 2 = 0. In Fig. 3 we show that the shape of these correlations is to leading order independent of y + in the jet region, presenting ⟨ψ ′ 1 ψ ′ 2 ⟩(∆y, ∆x = 0, y + ) and ⟨ψ ′ 1 ψ ′ 2 ⟩(0, ∆x, y). Note that according to [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF] setting y + = 0 or ∆x = 0 or ∆y = 0, as we do in Fig. 2 In the upper panel of Fig. 4 we present the structure of the correlation function as a function of ∆x and ∆y at a few fixed y + . According to [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF] the zero modes C(∆x, ∆y) can be observed by setting y + = 0, as presented in the center panel in Fig. 4. Qualitatively, from Fig. 4 it appears that the zero mode is symmetric with respect to reflection of ∆x, as expected. To quantify this we decompose the two-point function into its even and odd parts with respect to ∆x (and separately ∆y) the decomposition given by:

G ∆z even = G(∆z, ...) + G(-∆z, ...) 2 , (19)
G ∆z odd = G(∆z, ...) -G(-∆z, ...) 2 . ( 20
)
To quantify the symmetry of the zero mode we compute the relative power:

R ∆z [G] = d 2 xG 2 ∆z odd d 2 xG 2 ∆z even , (21)
at y + = 0. For all simulations considered we get that

R ∆x [⟨ψ ′ 1 ψ ′ 2 ⟩] ≪ 1 (as well as R ∆y [⟨ψ ′ 1 ψ ′ 2 ⟩] ≪ 1
) at y + = 0 (and in fact for any y + inside the jet region). Specifically, in the case of the simulation considered in Figs. [START_REF] Pope | Turbulent Flows[END_REF][START_REF] Marston | Recent developments in theories of inhomogeneous and anisotropic turbulence[END_REF][START_REF] Kraichnan | Inertial Ranges in Two-Dimensional Turbulence[END_REF]

we get R ∆x [⟨ψ ′ 1 ψ ′ 2 ⟩] = 0.092 ( L/4 0 L/4 x L/2 L/4 0 L/4 L/2
R ∆y [⟨ψ ′ 1 ψ ′ 2 ⟩] = 0
.099) at y + = 0. Thus, our results are consistent with the presence of zero modes of the form predicted in [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF], and support that these zero modes are even with respect to PT as expected from theoretical considerations.

For y + ̸ = 0 one expects contributions both from the even and the odd correlators, however, in practice we observe that the two-point function is independent of y + (as we expect for the zero mode) and looks identical to that at y + = 0. Thus it seems that the two-point function is dominated by the zero mode. In the lower panel of Fig. 4 we present the part even with respect to ∆x of the correlation function. The lower and upper panels indeed appear identical, meaning that the full correlation function is dominated by the even part. We note that from DNS we also get that R ∆y,∆x [⟨ψ ′ 1 ψ ′ 2 ⟩] = 0.102, while it should vanish from the exchange symmetry, suggesting that the odd contribution to ⟨ψ ′ 1 ψ ′ 2 ⟩ is comparable to the numerical noise.

An important question is how does the level of fluctuations scale with the parameters of the problem. In particular, there are at least two small parameters which can be important, l f /L and δ. In order for the perturbation theory to be consistent, fluctuations should be suppressed compared to the mean flow, the ratio expected to scale as a power of a small parameter. To analyze the scaling of the fluctuations we will focus on single point quantities, but instead of considering ⟨ψ ′2 ⟩ we will consider ⟨u ′2 ⟩ = ⟨(∂ y ψ ′ ) 2 ⟩ and ⟨v ′2 ⟩ = ⟨(∂ x ψ ′ ) 2 ⟩. First, since U = -∂ y Ψ = Const this will make the comparison to the mean flow much cleaner. In addition, it will allow to compare between the scaling of even under PT correlators and that of the odd correlator ⟨u ′ v ′ ⟩ = -⟨∂ x ψ ′ ∂ y ψ ′ ⟩ (while ⟨ψ ′2 ⟩ is an even correlator so there is no odd part to compare to).

For the latter we can derive an analytic expression:

⟨v ω′ y ψ ′ ⟩ = ⟨(∂ x ∇ 2 ψ ′ )ψ ′ ⟩, = -⟨∂ 2 x ψ ′ ∂ x ψ ′ ⟩ -⟨∂ 2 y ψ ′ ∂ x ψ ′ ⟩, = -∂ x ⟨(∂ x ψ ′ ) 2 ⟩ 2 + ⟨∂ y ψ ′ ∂ x ∂ y ψ ′ ⟩ + ∂ y ⟨uv⟩, = ∂ x ⟨(∂ y ψ ′ ) 2 ⟩ 2 + ∂ y ⟨uv⟩ = ∂ y ⟨uv⟩, (22)
where we have used the homogeneity in x repeatedly.

From the leading order solution for the mass flux we thus know that ∂ y ⟨vu⟩ = √ αϵ, meaning that ⟨vu⟩ = √ αϵy which can be written as ⟨u ′ v ′ ⟩ = (ϵL) 2/3 (y/L)δ 1/2 . This is confirmed by our DNS, shown in Fig. 7(b).

On the other hand, we expect the even correlators to be determined by the zero modes, and using equation [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF] we have that ⟨v ′2 ⟩ = -∂ 2 ∆x C| (0,0) and ⟨u ′2 ⟩ = -∂ 2 ∆y C| (0,0) so that both variances are expected to be constant in the jet region. Fig. 5 confirms this expectation. Note in passing that it is not a-priori clear that we can use [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF] to compute single point correlators of the derivatives of ψ ′ . Indeed, the presence of a direct cascade of the kinetic energy (∇ψ ′ ) 2 implies that non-linear interactions become important for correlators of derivatives of ψ ′ at small enough distances ∆x, ∆y ≪ l f , which would invalidate the approximations leading to equation ( 14) and its solution [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF] for such correlators. However, as we will see in the next section, in the region where the mean flow is strong the direct cascade is arrested for LQG, which may explain why the solution [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF] can still be used.

While lacking a prediction from analytic considerations, we can use the DNS results to determine the scaling of the fluctuations with the parameters of the model. We find that ⟨v ′2 ⟩ and ⟨u ′2 ⟩ scale differently, probably because of the asymmetry introduced by the mean flow, and we therefore examine them separately. We find that the following scalings lead to a collapse of data with different run parameters:

⟨u ′2 ⟩ ∼ (ϵL) 2/3 δ -1/2 , ⟨v ′2 ⟩ ∼ (ϵL) 2/3 δ 1/4 . (23)
Fig. 5 demonstrates the collapse of the variance profile when normalized by this scaling for three runs where only δ is varied. In Fig. 6 we demonstrate the collapse for runs with varying forcing and viscous scale. Since the velocity variance inside the jet is uniform (Fig. 5), here we take the mean value in the middle of the jet as representative for the run.

All in all, we find using that the mean jet velocity U 2 ∼ (ϵL) 2/3 δ -1 . Thus, the perturbation theory is indeed consistent, with the fluctuations suppressed compared to the mean flow with powers of δ. It is also worth noting that there isn't one characteristic scaling for the fluctuations, but rather a hierarchy with ⟨u ′2 ⟩ ≫ ⟨v ′2 ⟩ ≫ ⟨u ′ v ′ ⟩, in particular the odd correlator is suppressed compared to the even ones. That was also the case for the condensate state in 2DNS [START_REF] Laurie | Universal Profile of the Vortex Condensate in Two-Dimensional Turbulence[END_REF][START_REF] Frishman | Turbulence Statistics in a Two-Dimensional Vortex Condensate[END_REF]. This emphasizes that one cannot straightforwardly use a kinetic theory approach and justify quasilinear dynamics based on the naive scaling for the fluc-tuations (coming from the odd correlator).

⟨u ′2 ⟩/U 2 ∼ δ 1/2 , ⟨v ′2 ⟩/U 2 ∼ δ 5/4 , ( 24
)
⟨u ′ v ′ ⟩/U 2 ∼ (y/L)δ 3/2 , (25)

IV. THE DIRECT CASCADE

So far we have discussed the potential energy balance and derived the equation for the fluctuations two-point function based on the assumption that potential energy is transferred to large scales, and the only effective way by which it is dissipated is by the formation of the condensate. At the same time, we expect there to be a direct cascade of kinetic energy to small scales and we do not expect the condensate to have a significant influence on this process. Below we will show that it is not the case for LQG turbulence.

A. Spatial kinetic energy balance

We first derive the spatial kinetic energy balance. To obtain the total kinetic energy balance we act with (∂ i ψ)∂ i on (2). Note that from here, summation is implied over repeated indices. For the non-linear term we have:

∂ i ψ∂ i ∂ j (v ω j ψ) = ∂ i (∂ i ψ∂ j (v ω j ψ)) -ω∂ j (v ω j ψ) = ∂ i ∂ i ψ∂ j (v ω j ψ) -ωv ω i ψ = ∂ i J i , (26)
where we have used that since v ω i = ϵ ij ∂ j ω, from symmetry

ψv ω i ∂ i ω = ψϵ ij ∂ j ω∂ i ω = 0. ( 27
)
As expected from the inviscid conservation of kinetic energy, this contribution takes the form of a divergence of a flux (of kinetic energy) which we denote by J . We then decompose the stream-function into its mean and fluctuations ψ = Ψ + ψ ′ , average, and assume homogeneity in x (implying there is only a flux in the y direction):

⟨J y ⟩ = ⟨∂ y ψ∂ j (v ω j ψ)⟩ -⟨ωv ω y ψ⟩ = ∂ y Ψ∂ y ⟨v ω′ y ψ ′ ⟩ + ∂ y Ψ⟨v ω′ y ∂ y ψ ′ ⟩ + ⟨v ω′ j ∂ y ψ ′ ∂ j ψ ′ ⟩ -∂ 2 y Ψ⟨v ω ′ y ψ ′ ⟩ -Ψ⟨ω ′ v ω′ y ⟩ -⟨ψ ′ ω ′ v ω′ y ⟩. ( 28
)
The term proportional to Ψ in fact vanishes as ⟨ω ′ v ω′ y ⟩ = ⟨ω ′ ∂ x ω ′ ⟩ = ∂ x ⟨ω 2′ ⟩/2 = 0. Thus, the spatial kinetic energy flux is given by:

⟨J y ⟩ =∂ y Ψ∂ y ⟨v ω′ y ψ ′ ⟩ + ∂ y Ψ⟨v ω′ y ∂ y ψ ′ ⟩ -∂ 2 y Ψ⟨v ω ′ y ψ ′ ⟩ + ⟨v ω′ j ∂ y ψ ′ ∂ j ψ ′ ⟩ -⟨ψ ′ ω ′ v ω′ y ⟩.
It is straightforward to compute the remaining linear terms and the resulting steady state balance of kinetic energy can finally be written as

∂ y [⟨J y ⟩ + I D ] = η -D, (29)
where η = ⟨∂ i ψ ′ ∂ i f ′ ⟩ is the kinetic energy injection rate, D is the kinetic energy dissipation rate (expected to be mainly due to hyper-viscous dissipation of the fluctuations) and I D is the flux due to diffusion (e.g. for the drag

I Dα = α∂ y ⟨(∂ i ψ) 2 ⟩/2).
It is also useful to write the kinetic energy balance for the fluctuations. For the non-linear term the contribution can be computed by subtracting

∂ i Ψ∂ i ∂ j ⟨v ω′ j ψ ′ ⟩ = ∂ y Ψ∂ 2 y ⟨v ω′ y ψ ′ ⟩ from ∂ y ⟨J y ⟩.
In particular, for the terms involving the mean flow ∂ y Ψ we have

∂ y ∂ y Ψ∂ y ⟨v ω′ y ψ ′ ⟩ + ∂ y Ψ⟨v ω′ y ∂ y ψ ′ ⟩ -∂ y Ψ∂ 2 y ⟨v ω′ y ψ ′ ⟩ = ∂ y ∂ y Ψ⟨v ω′ y ∂ y ψ ′ ⟩ + ∂ 2 y Ψ∂ y ⟨v ω′ y ψ ′ ⟩. (30)
The first term on the bottom line is (part of) a flux of fluctuating kinetic energy, while the second term is (minus) the transfer term of kinetic energy between the fluctuations and the mean flow which we denote by T . We therefore get

∂ y J ′ y + I ′ D = η -D ′ + T, (31)
where

J ′ y ≡ ⟨J y ⟩ -∂ y Ψ∂ y ⟨v ω′ y ψ ′ ⟩ ( 32
)
is the flux of kinetic energy of the fluctuations, I ′ D is the fluctuating flux due to diffusion and D ′ is the dissipation rate of kinetic energy fluctuations. We expect T ≡ -∂ 2 y Ψ∂ y ⟨v ω′ y ψ ′ ⟩ = ∂ y U ∂ y ⟨v ω′ y ψ ′ ⟩ ≥ 0 so that the kinetic energy is transferred from the mean flow to the fluctuations. Note that by an order of magnitude estimate we expect the transfer term and the difference between the total flux and the fluctuating flux ⟨J y ⟩ -J ′ y to be of order ϵ/L 2 ≪ η.

We can now use the leading order solutions (3), ( 4), which imply that ⟨ω⟩ = 0 and V ω y = 0 as well as that ⟨J y ⟩ = J ′ y and the transfer term vanishes, since both ∂ y Ψ and ⟨v ω′ y ψ ′ ⟩ are independent of y to leading order. Therefore

⟨J y ⟩ = J ′ y = ∂ y Ψ⟨v ω′ y ∂ y ψ ′ ⟩ + ⟨v ω′ j ∂ y ψ ′ ∂ j ψ ′ ⟩ -⟨ψ ′ ω ′ v ω′ y ⟩, (33)
and the balance for the fluctuations reads

∂ y ∂ y Ψ⟨v ω′ y ∂ y ψ ′ ⟩ + ⟨v ω′ j ∂ y ψ ′ ∂ j ψ ′ ⟩ + ⟨ψ ′ ω ′ v ω′ y ⟩ + I ′ D = η -D ′ . ( 34
)
We can now directly evaluate the spatial flux of kinetic energy mediated by the mean flow, computing ⟨v ω′ y ∂ y ψ ′ ⟩ based on our previous results for the two-point function. Indeed, it is given by the limit r 2 → r 1 of the two-point function

v ω′ y (r 1 )∂ y ψ ′ 2 = ∂ 2 y1 + ∂ 2 x1 ∂ x1 ∂ y2 ⟨ψ ′ 1 ψ ′ 2 ⟩ . ( 35
)
Note that we only need to consider the inhomogeneous part of the solution [START_REF] Kolokolov | Structure of coherent vortices generated by the inverse cascade of twodimensional turbulence in a finite box[END_REF] as the zero mode is even under the reflection symmetry x → -x, while there is an odd number of derivative with respect to x 1 appearing above. Thus, the contribution from the zero mode will be odd and vanish in the limit r 2 → r 1 . Carrying out the calculation, we get

v ω′ y (r 1 )∂ y ψ ′ 2 = α ϵ ∇ 2 1 ∂ x1 ∂ y2 (y 1 + y 2 ) ∆x 0 dz ∆y/2 0 dz ′ χ12 (z, z ′ ) + {odd} = α ϵ ∇ 2 1 ∂ y2 (y 1 + y 2 ) ∆y/2 0 dz ′ χ12 (∆x, z ′ ) + {odd} = α ϵ ∇ 2 1 - y 1 + y 2 2 χ12 (∆x, ∆y) + ∆y/2 0 dz ′ χ12 (∆x, z ′ ) + {odd} = -y + α ϵ ∇ 2 1 χ12 + {odd}, (36)
where ∆x = x 1 -x 2 , ∆y = y 1 -y 2 and we have used that the correlation function is even with respect to ∆x → -∆x, and ∆y → -∆y. In writing [START_REF] Rhines | Waves and turbulence on a beta-plane[END_REF] we have stated explicitly only the terms which will contribute in the limit r 2 → r 1 , suppressing the odd contributions. To take the limit we must evaluate ∇ 2 1 Φ12 in this limit, which can be obtained from the definition of the kinetic energy

injection rate η η = ⟨∂ i f ∂ i ψ⟩ = -⟨f ∇ 2 ψ⟩, = -lim 2→1 ∇ 2 1 ⟨ψ ′ 1 f 2 ⟩ = -lim 2→1 ∇ 2 1 χ 12 . (37)
We obtain ⟨v ω′ y ∂ y ψ ′ ⟩ by taking ∆x, ∆y → 0 of (36) using [START_REF] Eyink | Local energy flux and the refined similarity hypothesis[END_REF] and get:

⟨v ω′ y ∂ y ψ ′ ⟩ = η α ϵ y = √ ϵα l 2 f y, (38)
which is in good agreement with DNS, as presented in Fig. 7(a). Note that here we have assumed that there is a negligible amount of kinetic energy injected into the k x = 0 and k y = 0 modes (but have not assumed that the forcing is isotropic) so that η = ∇ 2 1 Φ12 ≈ ∇ 2 1 Φ 12 = η. As a whole we thus get that at leading order

∂ y J ′ y = ∂ y Ψ∂ y ⟨v ω′ y ∂ y ψ ′ ⟩ ≈ η. (39)
We could have inferred that this flux term will give a contribution of the order of η based on an order-ofmagnitude estimate. To see this, we first recall that ⟨v ω′ y ψ ′ ⟩ = ∂ y ⟨vu⟩, where u = -∂ y ψ ′ , v = ∂ x ψ ′ . On the other hand, we would like to evaluate ∂ y ⟨v ω′ y ∂ y ψ ′ ⟩ = -∂ y ⟨v ω′ y u⟩. As the fluctuations are determined at small scales, we thus expect that derivatives acting on the fluctuating fields inside the average will get a contribution from scales ∼ l f , which leads to the estimate ∂ y ⟨v ω′ y u⟩ ∼ ∂ y ⟨vu⟩/l 2 f ∼ √ αϵ/l 2 f in agreement with equation [START_REF] Chen | Physical mechanism of the two-dimensional enstrophy cascade[END_REF]. Note that the sign of the flux, implying that it carries kinetic energy away from the jet region, seems to be a non-trivial result of the calculation. The direction of the flux is evidently linked to the direction of transfer of potential energy between scales: when potential energy is transferred from the fluctuations to the mean flow, that suppresses the direct cascade in that region, and kinetic energy is carried away from this region.

The result for J ′ y , [START_REF] Chen | Physical mechanism of the two-dimensional inverse energy cascade[END_REF], suggests that all the kinetic energy which is injected locally is carried away by a spatial flux due to the presence of the mean flow. In particular, if there is no spatial flux due to non-linear fluctuationsfluctuations interactions which brings kinetic energy to this region from other regions, this implies that the dissipation of kinetic energy in the jet region is negligible, D ′ ≪ η. This is indeed in agreement with our results from DNS, as can be seen from Fig. 8 where the profiles of the terms in the kinetic energy balance (31) are shown. In the jet region we indeed see that the balance is between the kinetic energy injection and the divergence of the flux ∂ y (U ⟨⟨v ω ′ y u ′ ⟩) due to the mean flow, in accordance with Eq.( 39). This implies that the kinetic energy is carried away from the jet region, where the mean flow is strong, before it has time to cascade to small scales and dissipate there -so that the mean flow effectively arrests the direct cascade. The kinetic energy is then deposited in the region in between the jets where the divergence of this flux becomes negative in most of the region as seen in Fig. 8. This is also the region where dissipation of kinetic energy occurs. Note however, that other terms in the flux J ′ y also become important in that region, redistributing kinetic energy in the opposite direction to that of U ⟨⟨v ω ′ y u ′ ⟩, as seen in the red curve in Fig. 8. In particular cubic terms in fluctuations (not shown separately here) have an important contribution to the flux, which is probably related to the presence of a coherent vortex in that region.

B. Local scale-to-scale flux: the filtering approach

The presence of a condensate makes our problem inhomogeneous due to the effects of the large-scale mean flow. In the previous section, we have shown indications that this inhomogeneity affects the transfer of kinetic energy to small scales so that we expect the direct cascade to proceed inhomogeneously in space. In this section, we would like to confirm this scenario by directly examining the kinetic energy flux between scales for different regions in the flow. The flux in Fourier space only gives the mean flux for the entire flow, so cannot differentiate between different spatial regions. Instead, we employ a real space filtering technique [START_REF] Eyink | Localness of energy cascade in hydrodynamic turbulence. I. Smooth coarse graining[END_REF], combining local spatial information and information about transfer between scales. It relies on a coarse-graining of the fields in real space using a convolution kernel with a characteristic length scale. The convolution with the kernel effectively filters features on scales smaller than its length scale. This approach is somewhat similar to a 2D wavelet transform, that keeps the spatial dependence. It then allows for the decomposition of the fluid kinetic energy (or other quadratic integrals) into band-pass contributions from a series of length scales in real space, writing the corresponding budget equation gives the transfers of turbulent energy both in space and in scale. The main feature of this approach which will be useful here is the scale-to-scale flux term which is space dependent in this approach. It will allow us to determine the spatial distribution of the flux across scales of potential and kinetic energy. The approach was previously applied to incompressible Navier-Stokes [START_REF] Eyink | Localness of energy cascade in hydrodynamic turbulence. I. Smooth coarse graining[END_REF][START_REF] Eyink | Local energy flux and the refined similarity hypothesis[END_REF][START_REF] Chen | Physical mechanism of the two-dimensional enstrophy cascade[END_REF][START_REF] Chen | Physical mechanism of the two-dimensional inverse energy cascade[END_REF], as well as other kinds of flows, including compressible flows [START_REF] Aluie | Scale decomposition in compressible turbulence[END_REF]. Here we will perform the scale decomposition and derive the analogous balance equations for the LQG equations.

Following [START_REF] Eyink | Localness of energy cascade in hydrodynamic turbulence. I. Smooth coarse graining[END_REF], we define a smooth low-pass filter as

ψ l (r) ≡ dx ′ G l (x ′ )ψ(x ′ + r), (40)
where the convolution kernel G l (x) is taken to be smooth, non-negative, normalized dr G l (r) = 1 and spatially localized. The filter scales with l as G l (r) = l -2 G(r/l). Specifically, we will choose the Gaussian kernel G l (r) = e -r 2 /2l /(2πl 2 ) when applying filtering to DNS. We may use the filtering operator to write equations for the filtered quantities with a given length scale. In doing so, non-linear interactions cause the emergence of terms representing energy transfer between scales. Filtering is a type of averaging, so the balance equations derived in this section for the large scale kinetic and potential energy are identical in structure to those one derives for the mean flow (where the averaging is over time). In particular, acting with the filter on (2) results in the filtered equations of motion

∂ τ ψ l + ∇ v ω l ψ l + ξ l = f l + α∇ 2 ψ l -ν(-∇ 2 ) p ψ l , ( 41
)
where ξ is the space dependant flux of the stream function to small scales, defined as

ξ l ≡ (v ω ψ) l -v ω l ψ l . ( 42
)
Note that aside from the additional, small scale, spatial flux term ξ l , Eq. ( 41) is the same as the regular LQG equation ( 2) (with ψ replaced by ψl ).

We are interested in the potential e l ≡ 1 2 ψ l 2 and kinetic h l ≡ 1 2 ∂ i ψ l 2 energy balance. We start with the potential energy flux, obtained by multiplying ( 41) by ψ l and writing the non-linear terms as a divergence of a flux and a transfer term between scales:

∂e l ∂τ + ∇ • J e l = P e l -Π l -D e l , (43)
where J e l is a spatial flux term of large-scale energy, P e l is the production of large-scale energy, D e l is the dissipation of energy at large scales and Π l is the scale-to-scale energy flux, positive if the transfer is out of the large scales to small scales. The terms are given by:

Π l = -∇ψ l • ξ l , ( 44
) J e l = v ω l e l + ψ l ξ l -α∇e l + νI e,p l , (45)
D e l = α ∂ i ψ l 2 + ν ∂ i1 • • • ∂ ip ψ l 2 , ( 46
)
P e l = ψ l f l , (47)
where I e,p l is the spatial transport due to hyper-viscosity

∇ • I e,p l ≡ ψ l (-∇ 2 ) p ψ l -(∂ i1 • • • ∂ ip ψ l ) 2 .
Similarly, to derive the balance for the kinetic energy one takes the derivative of (41) ∂ i , and multiplies it by ∂ i ψ l , which gives the equation for low-pass kinetic energy density balance

∂ ∂τ h l + ∇ • J h l = P h l -Z l -D h l , (48)
where J h l is the spatial flux of large-scale kinetic energy (compare the nonlinear contribution to (26)), P h l is the production of kinetic energy at large scales, D h l is the dissipation of kinetic energy at large scales and Z l is the scale-to-scale kinetic energy flux. The different terms are given by:

Z l = -(∂ j ∂ i ψ l )∂ i ξ j , (49)
J h l = (∇ψ l )v ω j ∂ j ψ l -v ω j ψ l ∇ 2 ψ l + ∂ i ψ l ∂ i ξ l -α∇h l + νI h,p l , (50)
D h l = α ∂ i1 ∂ i2 ψ l 2 + ν ∂ i1 • • • ∂ ip+1 ψ l 2 , ( 51
)
P h l = ∂ i ψ l ∂ i f l , (52)
with ∇I h,p l

≡ ∂ i ψ l (-∇ 2 ) p ∂ i ψ l -(∂ i1 • • • ∂ ip+1 ψ l ) 2 .
The transfer terms Π l and Z l can either be positive, transfering energy from scale l to smaller scales (acting as a sink) or negative, transferring energy from small scales to l (acting as a source). The direct cascade of kinetic energy corresponds to a positive flux (large to small scales), i.e. one expects Z l > 0 on average, while an inverse transfer of potential energy implies Π l < 0 on average. In addition, we expect that Z l ≈ 0 for large enough scales l, as kinetic energy is transferred from the forcing scale to smaller scales, and similarly that Π l ≈ 0 for small enough l. The indications above that the direct cascade does not occur in regions where the jets are strong leads to the expectation that Z l ≈ 0 in those regions and that the transfer of kinetic energy to small scales is concentrated in between the jets.

We show a short time average over 15T L of Π l and Z l in Fig. 9, which demonstrates their spatial distribution. In previous studies of the inter-scale flux in turbulent flows, the flux was observed to be statistically isotropic (as expected) and had regions of both positive and negative contributions on the level of a single snapshot [START_REF] Rivera | The direct enstrophy cascade of two-dimensional soap film flows[END_REF]. A definite sign thus emerged only upon averaging. This is roughly what we observe for the potential energy interscale flux Π l averaged over short times, Fig. 9(a)(b). In the panel below, we also show the flux when averaged in the x (homogeneous) direction. Then, though the result is still highly fluctuating, a negative flux, on average independent of y in the jet region, emerges for the larger coarse-graining scale l. At the small scales, as expected the flux is homogeneous and fluctuating around zero. Note that at the large scales there is an imprint of the jet structure on the flux, with stronger fluctuations in the inter-jet region.

For the kinetic energy flux Z l the distinction between jet regions and inter-jet regions is evident already at the level of a short-time average, Fig. 9(c)(d). In particular, the flux is visibly suppressed in the jet regions, and for small l is mostly positive between them (rather than having spatially distributed patches of positive and negative contributions of almost comparable magnitude). At large scales l larger fluctuations can be seen in between the jets, but a definite sign is harder to distinguish. These observations are further quantified in the panel below, where upon averaging in the x direction the difference between the two regions is even more clearly seen for the smaller scale l. For the larger l, the flux fluctuates around zero in the jet region while in between jets a very small negative flux emerges (related to the kinetic energy of the mean flow in that region, which has large gradients there). Thus, we see that most of the direct cascade is indeed concentrated in regions between jets, where the bias between a positive and a negative transfer is so much amplified.

Finally, to more systematically quantify the effects observed in Fig. 9 we consider the spatially averaged interscale fluxes with varying coarse-graining scale l. We average both in space and in time, starting once the simulations reach the steady state and up to 100T L . To examine the difference in Z l between spatial regions inside and outside the jets, we split the spatial average of Z l (x, y) into the jet region A jet and the region outside the jet A inter-jet . We choose A jet as the region where the leading order solution for the mean flow applies, as also used in previous figures, and A non-jet = A -A jet with A being the entire domain. We can also define the total average flux (averaged over the whole domain) given by: ⟨Z l ⟩ = (A jet /A) ⟨Z l ⟩ | Ajet + (A inter-jet /A) ⟨Z l ⟩ | Ainter-jet . For truly homogeneous turbulent flow, the partition would not affect the measurement (as long as both parts are large enough so that the statistics are comparable or the averaging time is long enough). The results are presented in Fig. 10. As per our expectation, the potential energy inter-scale flux ⟨Π l ⟩ is negative everywhere, corresponding to an inverse transfer, while ⟨Z l ⟩ is everywhere pos-itive (up to a very slight negative flux for l/L > 0.3 in the inter-jet region), as expected for a direct cascade. Moreover, the inter-scale flux from small scales l ≲ 0.1 is significantly suppressed in the jet region, implying that so is the direct cascade. This means that the spatial flux J h l dominates over the inter-scale flux Z l in the jet region at scales smaller than the forcing scale. Furthermore, at the smallest scales we observe that the total inter-scale flux is completely dominated by the inter-jet region (which occupies a smaller area fraction), in agreement with our observation that the overwhelming majority of the dissipation occurs there, Fig. 8. The presence of the mean flow also affects the potential energy inter-scale flux ⟨Π l ⟩ at large scales l/L > 0.3, though less dramatically. We observe that the inter-scale flux is reduced in between jets at large enough scales. This is probably a consequence of the inverse transfer being mostly mediated by the mean flow, which takes the form of a vortex in that region. The size of the vortex being of the order of 0.1L may thus explain the observed decrease. We have seen that the condensate has a dramatic effect on the direct cascade in LQG, an effect that appears to be absent in 2DNS [START_REF] Frishman | Jets or vortices-What flows are generated by an inverse turbulent cascade?[END_REF]. That raises the question of what determines for which types of 2D flows (with an inverse cascade) the latter effect could occur. In particular, recall that LQG and 2DNS are part of a wider class of active scalar equations where a scalar q is advected by a velocity with stream-function ϕ, with the relation between the two given by q k = |k| m ϕ k [START_REF] Pierrehumbert | Spectra of local and nonlocal two-dimensional turbulence[END_REF]. Here m controls the range of the dynamics, for 2DNS, (m, q, ϕ) = (2, ω, ψ) and for LQG (m, q, ϕ) = (-2, ψ, ω) so the velocity is given by derivatives of the scalar, making the dynamics local. All these flows have two positive definite conserved quantities, we shall call E = 1 2 qϕd 2 x and Z = 1 2 q 2 d 2 x, where Z cascades to small scales while E is transferred to large scales.

We have seen that for LQG the arrest of the direct cascade occurs in the regions where the mean flow is strong due to a spatial flux of Z away from those regions. A natural question is whether this mechanism could occur for other active scalar flows. To answer this question we consider the balance of Z for the turbulent fluctuations. First note that there are two types of terms involving the mean flow which enter this balance and can interfere with a homogeneous direct cascade: a transfer term between the mean flow and the fluctuations, and a spatial flux term. We expect Z to be transferred to small scales, and therefore an exchange term would tend to remove Z from the mean flow and transfer it to the fluctuationsenhancing the Z injection into the fluctuations in the regions of a strong mean flow. Thus, it is only a spatial flux term which could arrest the direct cascade as in LQG.

We now show that a spatial flux of Z fluctuations due to the mean flow is absent in models where interactions are non-local. In particular, we demonstrate this both for SWQG with L d > l f (for L d ≪ l f the deformation radius influences the direct cascade and we expect a transition to the LQG regime for scales L d < l < l f ), and for an active scalar with m > 0, assuming the flow is statistically homogeneous in the direction of the mean flow (i.e. that there is no trivial spatial flux due to the inhomogeneity of the turbulence). This is a consequence of the Z = q 2 balance for the fluctuations in the steady state:

∂ i u ′ i q ′2 2 + ∂ i Q⟨u i q ′ ⟩ + U i ∂ i q ′2 2 = η -D, (53)
where the third term (which is a spatial flux of Z ′ = q ′2 , due to advection by the mean flow) vanishes for a flow statistically homogeneous in the direction of U . Thus, the feedback between the condensate and the direct cascade as we have demonstrated in LQG does not exist for an active scalar with m > 0 which has long-range interactions, but might ,exist in models with m < 0 where small scale interactions are amplified. Finally, let us also discuss if the transfer of Z from the mean flow to the fluctuations, ∂ i Q⟨u i q ′ ⟩ , could significantly enhance the direct cascade in regions of strong mean flow (or large Q gradients). That requires for this term to be of order η (the Z injection rate), which we now show is not the case in 2DNS and SWQG. To estimate it let us assume a jet geometry for simplicity, giving ∂ y Q⟨u y q ′ ⟩ ≡ ∂ y Q⟨vq ′ ⟩ where we denote u = -∂ y ϕ ′ , v = ∂ x ϕ ′ as we had above. For SWQG (and 2DNS) we have ⟨vq ′ ⟩ = -∂ y ⟨uv⟩:

⟨vq ′ ⟩ = ⟨∂ x ψ ′ (∇ 2 -L -2 d )ψ ′ ⟩ = -⟨ψ ′ ∂ x ∇ 2 ψ ′ ⟩ = -∂ y ⟨uv⟩, (54)
where we have already demonstrated the last equality (the Taylor identity) in equation [START_REF] Rivera | The direct enstrophy cascade of two-dimensional soap film flows[END_REF] above. Thus, an order of magnitude estimate provided that L d > l f gives ∂ y Q⟨vq ′ ⟩ ∼ U ′′′ ∂ y ⟨uv⟩ ∼ ϵ/L 2 ≪ ϵ/l 2 f = η where ϵ is the injection rate of E, meaning that the transfer term is small. This is consistent with the observations in 2DNS [START_REF] Frishman | Jets or vortices-What flows are generated by an inverse turbulent cascade?[END_REF], where the cubic-in-fluctuations spatial flux of Z was more significant (though still small) compared to the transfer term.

V. DISCUSSION

In this work we have characterized the second order statistics of a jet condensate forming in the large-scalequasi-geostrophic equation, where potential energy experiences an inverse transfer while kinetic energy cascades to small scales. We have demonstrated that in the regions where the jets are strong the quasi-linear approximation is sufficient to obtain the second-order, two-point correlation functions of all the fluctuating fields (ψ and its derivatives). This is the case since the direct cascade is effectively arrested in those regions, so that non-linear fluctuation-fluctuation interactions are unimportant even for the kinetic energy (and thus can be neglected when determining e.g. correlators of ∇ψ). Using a local coarsegraining approach we have shown that the direct cascade is indeed mostly limited to the inter-jet regions. In the regions where the jets are strong, there is instead a spatial flux of kinetic energy, mediated by the mean flow, which prevents the direct cascade from developing, and which carries the kinetic energy to the inter-jet regions. At the same time, in between the jets we find regions where the quasi-linear approximation for the potential energy necessarily cannot work, since the mean-flow-fluctuations interactions (proportional to U = -∂ y Ψ) in those regions are small, and there is no other quasi-linear terms which can facilitate a transfer between mean flow and fluctuations. This is a consequence of interactions being local in LQG, so that there are no non-local (e.g. pressure) terms related to the mean flow which can redistribute the energy. Thus we find that in LQG the domain can be decomposed into two distinct regions: one where the dynamics is quasi-linear both for potential energy and for kinetic energy and another where fluctuationfluctuation interactions overwhelm mean-flow-turbulence interactions for both, which is also where the direct cascade is concentrated. We argue that both phenomena are related to the locality of interactions in LQG, and does not occur for flows with long-range interactions, i.e. an active scalars with m > 0, as well as models with short-range interactions reaching beyond the forcing scale, namely SWQG with L d > l f . It remains to be seen if active scalars with m < 0 or SWQG with L d < l f , both having dominant interactions below the forcing scale, can exhibit an arrest of the direct cascade as we have found for the limiting LQG case (m = -2, L d = 0). More generally, understanding the similarities and differences in the condensate state between these two classes of flows away from the LQG limit is an interesting direction for future work.

For the regions where the quasi-linear approximation applies, we have found that fluctuations are suppressed compared to the mean flow with powers of δ, the parameter which quantifies the strength of the condensate. Furthermore, we find that different correlation functions scale differently with δ and that correlators which are odd with respect to parity+time reversal symmetry are significantly suppressed compared to even correlators. Such a hierarchy was previously observed in 2DNS [START_REF] Frishman | Turbulence Statistics in a Two-Dimensional Vortex Condensate[END_REF], and points to the fact that constructing a closed perturbative quasi-linear theory for the condensate may be a subtle issue, as it cannot simply rely on a uniform scaling for the fluctuations. Related to this issue, in this work we have determined that even correlators arise from zero modes of an advection operator. We found that these zero modes are homogeneous in the jet region, depending only on ∆x and ∆y. How exactly those modes are to be determined, including their scaling with δ, however, remains unclear and is left for future work. We first briefly remind the physical origin of the shallow water quasi-geostrophic equation, from which the large-scale quasi-geostrophic (LQG) system is derived. It describes a rotating shallow fluid layer, where the horizontal scale of the fluid motion, L, is assumed much larger than the layers' mean depth H, and which is under the influence of gravity g. Assuming a constant rotation rate Ω ẑ, and a characteristic velocity U , the ratio between inertia and the Coriolis force is given by the Rossby number Ro = U/ΩL. A perturbative expansion in Ro ≪ 1, while assuming Ro(L/L d ) 2 ∼ o(1) then gives the SWQG equation [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF]

∂ t q + v • ∇q = ∂ t q + J(ψ, q) = 0; q = ∇ 2 -L -2 d ψ, (A1)
where q is the potential vorticity, ψ is the stream fucntion with v = ẑ × ∇ψ, ω = ∇ 2 ψ = (∇ × v) ẑ is the vorticity, J(ψ, q) is the Jacobian operator defined as J(ψ, q) = ∂ x ψ∂ y q -∂ y ψ∂ x q = ϵ ij ∂ i ψ∂ j q with ϵ ij the 2D Levi-Civita symbol. The length scale L d = √ gH/2Ω is the Rossby deformation radius. Here, hydrostatic balance relates the variation in the layers' depth δh to the pressure gρ∇δh = ∇p. While geostrophic balance relates the stream-function to the pressure so that in total ψ = (g/Ω)δh, and there is a single equation for ψ. There are two quadratic invariants in the SWQG system: energy (potential + kinetic) d 2 x qψ and square potential vorticity d 2 x q 2 . As a consequence, it permits both a direct and an inverse cascade.

Including forcing f , friction α (linear drag on velocity) and (hyper) viscosity ν, and using the box scale L to nondimensionalize lengths we can write the SWQG equation as

∂ t ∇ 2 - L L d 2 1 L 2 ψ + ϵ ij ∂ i ψ∂ j ∇ 2 ψ = f -α∇ 2 ψ + ν(-∇ 2 ) p ∇ 2 ψ. (A2)
Taking the limit (L d /L) 2 → 0 formally gives ∂ t ψ = 0, which is a purely decaying system. Following [START_REF] Larichev | Weakly decaying turbulence in an equivalent-barotropic fluid[END_REF], to capture the emerging slow dynamics we will work in rescaled time τ = t(L d /L) 2 with the limit L d /L → 0, giving

-∂ τ ψ + L 2 ϵ ij ∂ i ψ∂ j ∇ 2 ψ = f L 2 -αL 2 ∇ 2 ψ + L 2 ν(-∇ 2 ) p+1 ψ. (A3)
Next, we define a new stream-function variable ψ = L 2 ψ and a corresponding forcing f = -f L 2 , drag α = αL 2 and viscosity ν = νL 2 . With the chosen scaling the relation between the stream-function and the height perturbation in the shallow water system becomes ψ = (gL 2 /Ω)δh. We thus arrive at the LQG equation 2.

Let us also demonstrate that the LQG equation can be consistently derived directly from the rotating shallow water equations in the geostrophic limit Ro → 0.

For a single-layer fluid, and including the Coriolis term, the inviscid shallow water equations (SW) are and with the definitions ψ ′ 0 ≡ δh ′ 0 /f ′ 0 and J(ω, ψ) = ∂ x ω∂ y ψ -∂ y ω∂ x ψ we obtain the dimensionless invicid LQG equation

∂ τ ′ ψ ′ + J(ω ′ , ψ ′ ) = 0. (A16)
Returning the dimensions using ψ ∼ U L we have that ∂ τ ψ + L 2 J(ω, ψ) = 0. We may absorb the additional factor L 2 by redefining the stream function as ψ = L 2 ψ and thus ω = ∇ 2 ψ = L 2 ω. Dropping the tilde notation we arrive at the LQG equation for the re-scaled streamfunction

∂ τ ψ + J(ω, ψ) = 0. ( A17
)
Appendix B: Solution for the two-point function

Here we describe solution to equation ( 14) for the leading order of the two-point function. The solution to [START_REF] Srinivasan | Zonostrophic Instability[END_REF] is the sum of the solution to the homogeneous equation and the particular solution to the inhomogeneous equation. We begin by describing the latter, assuming the forcing is homogeneous in y, so that χ 12 depends only on y -≡ ∆y/2 and ∆x. In principle, in the variables y + , y -, x 1 the equation can be straightforwardly integrated to obtain the inhomogeneous solution. However, there is a subtle point that has to do with the consistency of the perturbation theory for some of the fluctuations modes.

Consider the Fourier transform of equation ( 14) with respect to ∆x (equivalently x 1 ), denoting the corresponding wavenumber by k x , and with respect to ∆y, denoting the wavenumber by k y . In ∆x, this is possible since the mean flow solution is homogeneous in x and is applicable throughout the x direction, which is periodic. In ∆y, we assume that correlations decay with ∆y within the region of applicability of the mean flow solution (note that this is not necessarily the case in y + , and that y + is not a periodic coordinate since the leading order mean flow has a finite range of applicability in y). We see that for k x = 0 and k y = 0 while the left hand side of the equation turns to zero, the right hand side does not. Thus, modes of the forcing with k x = 0 or k y = 0 need to be treated separately, and equation ( 14) is not the leading order equation. This is easily understood for k x = 0: perturbations with k x = 0 are not advected by the mean flow, so cubic terms or dissipative terms must be important in balancing the injection of the forcing into such modes. That k x = 0 and k y = 0 modes cannot be treated in a quasi-linear approximation is a completely general statement for mean-flow-turbulence interactions, e.g. also for 2DNS [START_REF] Frishman | Turbulence Statistics in a Two-Dimensional Vortex Condensate[END_REF]. Therefore, in equation ( 14) we should subtract these modes from the forcing correlation function:

χ12 = χ 12 - 1 L y Ly 2 - Ly 2 χ 12 (∆x, s)ds- 1 L x Lx 2 -Lx 2 χ 12 (s, ∆y)ds. (B1)
Then, the inhomogeneous solution reads

⟨ψ ′ 1 ψ ′ 2 ⟩ inh = (y 1 + y 2 ) α ϵ ∆x 0 dz ∆y/2 0 dz ′ χ12 (z, z ′ ) .
(B2) In (B2) we choose the initial point of the integration to be at coincident points ∆y = 0, and ∆x = 0 which makes the inhomogeneous part symmetric with respect to the replacement r 1 → r 2 (i.e. even under reflection ∆x → -∆x, ∆y → -∆y). On the other hand, it is odd with respect to ∆x → -∆x (and ∆y → -∆y) separately. This is what we expect from the fact that the forcing combined with the mean flow break the parity+time reversal symmetry x → -x, t → -t of the system, see the discussion in [START_REF] Svirsky | Statistics of inhomogeneous turbulence in large scale quasi-geostrophic dynamics[END_REF].

Let us consider the influence of the subtraction of the modes with k x = 0 and k y = 0 from the forcing. Assume a typical forcing length scale l f , it will be convenient to denote χ 12 (∆x, ∆y) = ϵΦ( ∆x l f , ∆y l f ), such that Φ(0, 0) = 1. This gives

⟨ψ ′ 1 ψ ′ 2 ⟩ inh = 2y + l 2 f √ αϵ ∆x l f 0 dz ∆y 2l f 0 dz ′ Φ (z, z ′ ) , (B3) where Φ(z, z ′ ) = = Φ(z, z ′ ) - l f L y Ly 2l f - Ly 2l f Φ(z, s)ds - l f L x Lx 2l f -Lx 2l f Φ(s, z ′ )ds = Φ(z, z ′ ) -Φ(z, k y = 0) -Φ(k x = 0, z ′ ) (B4)
with Φ being the Fourier transform of Φ with respect to ∆x or ∆y. Here, Φ(z, z ′ ) has no modes with k x = 0 or k y = 0: Φ (z, k y = 0) = Φ (k x = 0, z ′ ) = 0. For the forcing we have been using in DNS, a direct calculation gives Φ (z, k y = 0) = 1 π cos 2πz and Φ (k x = 0, z ′ ) = 1 π cos 2πz ′ . Now, for ∆x, ∆y ≤ l f the replacement of Φ by Φ does not change the result at leading order: the difference between (B3) and the expression when Φ(z, z ′ ) is replaced by Φ(z, z ′ ) is of order O(l f /L) (after integration). Similarly, we expect that the contribution to two-point correlation functions from k x = 0 and k y = 0 modes of the forcing, finding which requires a fully non-linear treatment not carried out here, will be small, of order O(l f /L), compared with the leading order.

On the other hand, for e.g. ∆y ≈ L/2 (similarly for ∆x ≈ L/4) we notice that

Ly 4l f 0 dz ′ Φ (z, z ′ ) ≈ Ly 2l f 0 dz ′ Φ (z, z ′ ) = 1 2 Ly 2l f - Ly 2l f dz ′ Φ (z, z ′ ) , = L y 2l f Φ (z, k y = 0) , (B5)
where we have used that Φ is a decaying function of z ′ assuming that

Ly 2l f Ly 4l f dz ′ Φ (z, z ′ ) → 0 as L y /l f → ∞, that
we are working in the regime L/l f ≫ 1, and that Φ (z, z ′ ) is even in z ′ : Φ(z, z ′ ) = Φ(z, -z ′ ) (corresponding to the assumed statistical reflection symmetry y → -y of the forcing). We then have

∆x l f 0 dz Ly 2l f 0 dz ′ Φ (z, z ′ ) ≈ L y 2l f ∆x l f 0 dz Φ (z, k y = 0) = 0. (B6)
Thus, we see that the forcing influences two-point correlation functions only for ∆x, ∆y ≤ l f , where we can use Φ(z, z ′ ) ≈ Φ(z, z ′ ), while for ∆x, ∆y ≈ L the contribution from the inhomogeneous solution is negligible.

While the forcing provides the leading order contribution to the odd in ∆x part of the correlation function, corresponding to parity+time reversal symmetry breaking, the even contribution at leading order must come from the homogeneous solutions to (B3). Those are the zero modes of the advection operator

L 1 + L 2 = ∇ 2 1 ∂ x1 + ∇ 2 2 ∂ x2 = ∂ y+ ∂ y-∂ x1 :
⟨ψ ′ 1 ψ ′ 2 ⟩ hom = C(∆y, ∆x) + C 1 (y + , ∆x) + C 2 (y + , ∆y). (B7) To determine which zero modes contribute to the correlation function we need to take into account the boundary conditions. First, we assume that the fluctuations decorrelate as ∆y → L, as confirmed in DNS Fig. 2(b), implying that C 1 (y + , ∆x) = 0. Indeed, the odd and even parts of the correlation function should decay to zero separately in this limit. Also, since C(∆x, ∆y) is independent of y + while C 2 is independent of ∆x, C 1 (y + , ∆x) should separately decay to zero as ∆y → L, implying it must be identically zero. This gives

(B9)

For all DNS, a constant time step dt, different for each DNS was used (Table II), with the forcing amplitude normalized by √ dt so that energy injection is independent of it. The grid spacing dx = dy ≈ 0.05 is the same for all simulations considered at the default resolution of 64x128. To verify the adequacy of the choice of dt and dx, they are compared with the smallest physical time scale τ E (l ν ) and length scale l ν respectively as presented in Table II. While the temporal resolution is smaller by at least 4 orders of magnitude than τ E (l ν ), the grid spacing is relatively close to the Kolmogorov scale l ν . The large difference between the spatial and temporal resolutions required is due to the hyper-viscosity used in the evolution equation [START_REF] Pope | Turbulent Flows[END_REF]. It allows us to use a relatively large grid spacing (or low resolution), as for p = 7 the energy cutoff is extremely sharp leaving only a small fraction of the total energy at length scales l ν < l ≤ dx. To verify that the simulations are spatially fully developed and that this resolution is not too coarse, two highresolution simulations (Sim-B(*2) and Sim-B(*4)) were performed with the same parameters as Simulation-B but with dx = dy ≈ 0.025 and dx = dy ≈ 0.0125 corresponding to 128x256 and 256x512 resolutions respectively. The resulting condensate is exactly the same as for the lower resolution, as can be appreciated from the snapshot comparison in Fig. 11 and from the averaged terms in Fig. 12. Finally, to demonstrate that the choice of hyperviscosity in eq. ( 2) does not affect the mean flow condensate, we performed an additional low hyper-viscosity run with the same parameters as Sim-B(*2) but with p = 5 and ν = 6.3×10 -13 -Sim-B(p5). The value of ν was chosen so that the Kolmogorov scale for both Sim-B(*2) and Sim-B(p5) is l ν ≈ 0.2. The higher resolution (compared to Sim-B) is used in the comparison as the energy cutoff is not as sharp in the p = 5 case, requiring a larger separation of scales to ensure convergence. The comparison demonstrates that the choice of hyper-viscosity does not affect the mean flow, as presented in Fig. 13. Note that due to the longer time required to integrate the equations at the high resolution, only limited statistics were obtained amounting to ∼ 30T L . To make the comparison quantitative, the averaged terms presented in Figs. 12,13 are over 30T L for all simulations.

4 FIG. 1 .

 41 FIG. 1. LQG jet condensate, showing the velocity v = ẑ×∇ψ snapshot (a) and mean over time (b). The color corresponds to the velocity magnitude. (Simulation-B)

FIG. 2 .

 2 FIG. 2. The two-point correlation function ⟨ψ ′ 1 ψ ′ 2 ⟩ as measured in DNS (Simulation-B) with (a) ∆x = ∆y = 0 (b) ∆x = y+ = 0 and (c) ∆y = y+ = 0. The region where the leading order solution for the mean flow applies is delimited by dashed lines. (Simulation-B)

FIG. 3 .

 3 FIG. 3. The variation of the two-point function with y+ with: (a) ∆x = 0 and (b) ∆y = 0. (Simulation-B)

FIG. 4 .

 4 FIG. 4. The two-point function (averaged over x+) as a function of (∆x, ∆y) measured at different y+ in the jet region. (Simulation-B)

FIG. 5 .Re × 10 13 FIG. 6 .

 5136 FIG. 5. Rescaled variance of the velocity fluctuations in the direction (a) parallel to the jet (u ′ = -∂yψ ′ ) and (b) perpendicular to the jet (v ′ = ∂yψ ′ ), for different values of the parameter δ. The two regions where the leading order solution for the mean flow applies are delimited by vertical dashed lines.

FIG. 7 .

 7 FIG.7. The terms (a) ⟨v ω′ y ∂yψ ′ ⟩ and (b) ⟨vu⟩ as measured from the DNS (solid lines), rescaled and compared with their theoretical predictions (dashed line), for different values of the parameter δ. The two regions where the leading order solution for the mean flow applies are delimited by vertical dashed lines.

FIG. 8 .

 8 FIG. 8. Kinetic energy balance of the fluctuations (31). The dominant flux term in the jet region U ⟨v ω′ y y ′ ⟩ = ∂yΨ⟨v ω′ y ∂yψ ′ ⟩ is plotted separately from the other flux terms. The two regions where the leading order solution for the mean flow applies are delimited by vertical dashed lines. (Simulation-B)

8 FIG. 9 .

 89 FIG. 9. Potential energy flux between scales Π l normalized by the total injection rate ϵ for scales l/L = 0.11 (a) and l/L = 0.89 (b). Kinetic energy flux between scales Z l normalized by the total injection rate η for scales l/L = 0.11 (c) and l/L = 0.89 (d). Below each plot the average over x is shown. The fields shown are obtained upon taking a short-time average over 15TL in steady-state. The two regions where the leading order solution for the mean flow applies are delimited by vertical dashed lines. (Simulation-K)

FIG. 10 .

 10 FIG. 10. Average energy fluxes across length scale l inside (blue) and outside (orange) the jet region. The total flux is denoted by a dashed line (green). (a) Potential energy flux ⟨Π l ⟩. (b) Kinetic energy flux ⟨Z l ⟩.(Simulation-K)

  Appendix A: LQG from SWQG and consistency of limits

⟨ψ ′ 1 ψ ′ 2 ⟩0 2 0

 22 = C(∆y, ∆x) + C 2 (y + , ∆y) dz ′ χ12 (z, z ′ ) .(B8)Next, C 2 (y + , ∆y) is in fact a zero mode of the individual advection operators L i , reflecting the fact that k x = 0 modes of ψ ′ are not advected by a mean flow pointing in the x direction, irrespective of the shape of the mean flow (as discussed above). So, such contributions to the correlation function are not constrained by the quasi-linear approximation. While there is no a-priori reason to set them to zero, we may thus expect that modes with k x = 0 do not contribute significantly. Indeed, we see empirically in our DNS that when setting ∆y = 0 the even part of the two-point correlation function (with varying ∆x) is independent of y + , see Fig.2(a) and Fig.3(b). Thus, it does not contribute, at least to leading order. Finally, the full solution reads ⟨ψ ′ 1 ψ ′ 2 ⟩ = C(∆y, ∆x) dz ′ χ12 (z, z ′ ) .

4 FIG. 11 .

 411 FIG. 11. Comparison of the LQG jet condensate at two spatial resolutions, showing the velocity v = ẑ × ∇ψ snapshot (a) Simulation-B with dx = dy = 0.0491 and (b) Simulation-B(*2) with dx = dy = 0.0245.

FIG. 12 .

 12 FIG.12. The terms (a) ∂yΨ and (b) ⟨v ω′ y ∂yψ ′ ⟩ as measured from the DNS (solid lines), rescaled and compared with their theoretical predictions (dashed line), for simulations B, B(*2) and B(*4) with spatial resolution dx = dy = 0.0491, dx = dy = 0.0245 and dx = dy = 0.0123 respectively, with all other parameters the same. The profiles are obtained from short averaging over ∼ 30TL. The two regions where the leading order solution for the mean flow applies are delimited by vertical dashed lines.

FIG. 13 .

 13 FIG.[START_REF] Marston | Generalized Quasilinear Approximation: Application to Zonal Jets[END_REF]. The terms (a) ∂yΨ and (b) ⟨v ω′ y ∂yψ ′ ⟩ as measured from the DNS (solid lines), rescaled and compared with their theoretical predictions (dashed line), for simulations B(*2) and B(p5) with hyper-viscocity (ν = 7.3 × 10 -19 , p = 7) and (ν = 6.91 × 10 -13 , p = 5) respectively chosen such that kν ≈ 31.2 The two regions where the leading order solution for the mean flow applies are delimited by vertical dashed lines.

TABLE I .

 I Parameters of the DNS runs. All runs are performed with hyper-viscosity p = 7 on a 64×128 grid. The forcing wave number is k f , the drag coefficient α, viscosity ν, potential energy injection rate ϵ, the ratio of large eddy turnover time and dissipation time scales is δ = αϵ -1/3 L 2/3 , the ratio of viscous and forcing time scales is Re = l

						2p-8/3 f	ϵ 1/3 /ν and
	TL is the simulation time in units of large eddy turnover time
	τL = ϵ -1/3 L 8/3 .
	k f α×10 3 ν×10 19 ϵ×10 4	δ	Re×10 -13	TL
	A 13	2.0	7.30	2.44 .109	2.26	872.9
	B 13	1.0	7.30	2.42 .055	2.25	1031.7
	C 13	0.5	7.30	2.42 .027	2.25	1122.8
	D 13	1.0	7.30	1.04 .072	1.70	185.0
	E 13	1.0	0.50	2.59 .053	33.64	176.0
	F 10	2.0	10.00	1.35 .133	26.47	191.6
	G 10	1.0	10.00	1.37 .066	26.60	155.1
	H 10	0.5	10.00	1.36 .033	26.54	252.7
	I 15	1.0	7.30	2.48 .054	0.45	327.9
	J 15	2.0	0.50	2.71 .105	6.75	300.2
	K 15	1.0	0.50	2.74 .052	6.77	206.6

TABLE II .

 II Spatial and temporal resolutions used for integration of the simulations detailed in TableIwith and addition of simulations B( * 2), B( * 4) and B(p5) which are performed at higher spatial resolutions. The spatial resolutions dx are compared with the Kolmogorov scale lν while the temporal resolution dt is compared with the fastest timescale τE(lν ).

		dx	dt	lν /dx	τE(lν )/dt
	A	0.0491	4.0 • 10 -5	4.085	5.2 • 10 3
	B	0.0491	4.0 • 10 -5	4.086	5.2 • 10 3
	C	0.0491	2.0 • 10 -5	4.087	1.1 • 10 4
	D	0.0491	4.0 • 10 -5	4.185	7.4 • 10 3
	E	0.0491	1.6 • 10 -5	3.224	7.0 • 10 3
	F	0.0491	1.0 • 10 -4	4.279	2.9 • 10 3
	G	0.0491	4.0 • 10 -5	4.278	7.3 • 10 3
	H	0.0491	4.0 • 10 -5	4.279	7.3 • 10 3
	I	0.0491	2.0 • 10 -5	4.076	1.0 • 10 4
	J	0.0491	1.6 • 10 -5	3.217	6.7 • 10 3
	K	0.0491	1.6 • 10 -5	3.216	6.7 • 10 3
	B(*2)	0.0245	4.0 • 10 -5	8.169	5.2 • 10 3
	B(*4)	0.0123	4.0 • 10 -5	16.34	5.2 • 10 3
	B(p5)	0.0245	4.0 • 10 -5	8.149	5.2 • 10 3

where u = (u, v) is the horizontal velocity, h is the height of the upper free surface (where the bottom surface is assumed flat), f c = Ω ẑ and g is gravity. We apply the geostrophic scaling [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation[END_REF] -assuming that u = (u, v) ∼ U , (x, y) ∼ L and an advective time scale T ∼ L/U . We decompose the free layer height as h = h + δh with the mean height h = H = const. and the variation δh ∼ RoH(L/L d ) 2 . The assumptions so far are the same as those made to obtain SWQG. Here, we also rescale the time by τ = t(L d /L) 2 . We will consider two limits, Ro → 0 and (L d /L) → 0, and we must specify the relation between them. For what follows we assume that when both limits are taken Ro tends to zero faster than (L d /L). Thus we may assume that (L d /L) ∼ Ro b with 0 < b < 1/2. Note that with this scaling the height perturbations are still small compared to the mean height as δh/H ∼ Ro 1-2b ≪ 1. With this scaling, we obtain the non-dimensional SW momentum equation:

) and the non-dimensional SW height variation equation:

Having expressed both small parameters as a functions of Ro we expand the velocity u ′ , v ′ in ϵ i = ϵ i (Ro) such that 1 = ϵ 0 ≫ ϵ 1 ≫ ... and similarly we expand the height variation δh ′ in µ

We leave the asymptotic series arbitrary for now. Substituting the series (A8) into the re-scaled time momentum equation (A6) we get:

where f ′ c ≡ f ′ 0 ẑ = 1 ẑ. The dominant balance (for any ϵ 1 , µ 1 ≪ 1) is between the pressure and Coriolis force thus

The re-scaled mass conservation (A7) gives at leading order the same result. This allows for the definition of the stream function ψ ′ 0 ≡ δh ′ 0 /f ′ 0 . Moving on to the next order in perturbation theory to get the dynamics, we consider the next order of the momentum equation (A6).

Taking its curl gives the vorticity ω = ∇ × u equation

where the time derivative term has been neglected as it is of higher-order in Ro than the advection term. Note that the only non-trivial option, in this case, is for ϵ 1 = Ro and in general we may assume that ϵ n = Ro n . To proceed consider the next order of (A7). First, we note that (A10) gives (u ′ 0 ∇)δh ′ 0 = 0 and that Ro 1-2b ≪ 1, thus we get

Using the (∇ • u ′ 1 ) term to relate (A12) and (A13) we obtain

We wish to obtain a solution for which the leading order velocity does not vanish, thus ∂ τ δh ′ 0 must be determined at this order and µ 1 ≤ Ro 2b . If we assume µ 1 = Ro 2b ≫ Ro ≫ Ro 1+2b we obtain from (A11)∇δh ′ 1 = 0 and thus (u ′ 0 ∇)δh ′ 1 = 0. Therefore in any case the RHS term can be neglected and equation (A14) reduces to