# Correlations in a weakly interacting two-dimensional random flow

I.V. Kolokolov1,2, V.V. Lebedev1,2, V.M. Parfenyev1,2 1Landau Institute for
Theoretical Physics, RAS,
142432, Chernogolovka, Moscow region, Russia;
2National Research University Higher School of Economics,
101000, Myasnitskaya ul. 20, Moscow, Russia.

###### Abstract

We analytically examine fluctuations of vorticity excited by an external
random force in two-dimensional fluid. We develop the perturbation theory
enabling one to calculate nonlinear corrections to correlation functions of
the flow fluctuations found in the linear approximation. We calculate the
correction to the pair correlation function and the triple correlation
function. It enables us to establish the criterion of validity of the
perturbation theory for different ratios of viscosity and bottom friction. We
find that the corrections to the second moment are anomalously weak in the
cases of small bottom friction and of small viscosity and rely the weakness to
the energy balance and to the enstrophy balance. We demonstrate that at small
bottom friction the triple correlation function is characterized by a
universal scaling behavior in some region of lengths. The developed
perturbation method was verified and confirmed by direct numerical
simulations.

##  I Introduction

For a three-dimensional flow in a fluid, the only dissipation mechanism of the
kinetic energy is viscosity. Thus, there is the only dimensionless parameter,
Reynolds number, characterizing the degree of nonlinearity. In contrast, when
considering effectively two-dimensional thin fluid films one has to deal with
two dissipative mechanisms: viscosity and bottom friction. An interplay of the
dissipative mechanisms leads to a more complicated characterization of the
flow nonlinearity.

The traditional way to produce flow in $2d$ hydrodynamics is to apply an
external periodic static force (Kolmogorov force) to the fluid. Such setup is
used both in laboratory experiments with thin liquid films Sommeria86 ;
Tabeling02 ; Shats09 ; Orlov17  and in numerical simulations MCH04 ; Mishra15
; DFKL22 . Then the transition to turbulence appears to be complicated, it can
be soft transition or jump (depending on the ratio between the pumping length
and the box size) and can go through some intermediate stages. Aiming to
simplify the situation and to reveal its universal features, we analytically
consider the case when the external force exciting the flow is a random
function of space and time with homogeneous statistical properties. Our
preliminary results show that the transition to turbulence in this case occurs
gradually (without jumps). Therefore, there is a wide range of parameters
where the interaction of the random flow fluctuations is relatively weak.

The regime of relatively weak random flow can be easily realized numerically.
And we use direct numerical simulations to verify our analytical results.
Besides, the regime can be realized experimentally in thin fluid films by
applying a random fast varying in time external force to the fluid.
Unfortunately, we do not know works devoted to such setup. An interest to this
regime in thin fluid films may be related to the organization of a
controllable mixing process of impurities in the film. To characterize the
process one has to find the coefficient of turbulent diffusion.

In this paper, we investigate effects related to correlations in a weakly
interacting $2d$ random flow, excited by an external force, random in space
and time. The character of the flow interaction essentially depends on the
ratio of the bottom friction coefficient $\alpha$ and $\nu k_{f}^{2}$, where
$\nu$ is the kinematic viscosity coefficient and $k_{f}$ is the characteristic
wave vector of the external force. We examine an arbitrary ratio between
$\alpha$ and $\nu k_{f}^{2}$ in our work. Of course, arbitrary values of the
coefficients $\nu$ and $\alpha$ can be easily implemented in numerical
simulations of the two-dimensional fluid. As for experiments, the situation is
as follows.

For thin liquid films it is easy to achieve the inequality $\alpha\gg\nu
k_{f}^{2}$, since $\alpha$ can be estimated as $\nu h^{-2}$, where $h$ is the
thickness of the film. The characteristic length of the pumping force should
be larger than $h$ (otherwise the flow cannot be treated as two-dimensional),
thus $k_{f}h\gtrsim 1$. That leads to the conclusion $\alpha\gtrsim\nu
k_{f}^{2}$. However, to observe interesting effects related, e.g., to large-
scale vortices in the turbulent regime, one should make $\alpha$ as small as
possible. To achieve the goal some experimental tricks are used related to
multilayered films Tabeling02 ; Shats09 ; Orlov17 . The opposite inequality
$\alpha\ll\nu k_{f}^{2}$ can be achieved for the effectively two-dimensional
subsystem of the fast rotating three-dimensional fluid where the effective
bottom friction is related to formation of Ekman boundary layers Vergeles21 .
Other cases, where small $\alpha$ can be realized, include suspended soap or
smectic films, which have no contact with a solid boundary and therefore do
not exhibit a velocity gradient across the film thickness kellay2002two ;
Yablonskii2017 .

Our aim is to examine consistently first corrections to correlation functions
of the flow fluctuations. More precisely, we analyze corrections to the pair
correlation function of vorticity and its triple correlation function. It
enables one to establish the criterion of validity of the perturbation theory,
and, consequently, to predict conditions of exciting true turbulent state.
Note also, that the flow fluctuations are believed to be weakly interacting
inside the coherent vortices appearing in $2d$ turbulence KL15 ;
kolokolov2016structure ; kolokolov2016velocity ; frishman17 . This work will
help to clarify validity of the approach.

The rest of the paper is organized as follows. In Section II, we give the
general theoretical framework for the two-dimensional flow excited by an
external force. Section III is devoted to the development of perturbation
series enabling one to calculate effects related to the nonlinear interaction.
In Section IV, we calculate first corrections to the correlation functions. In
particular, in Subsection IV.1 we analyze the correction to the pair
correlation function of vorticity, thereby establishing parameters controlling
the perturbation series, and in Subsection IV.2 we examine the correction to
the second moment of vorticity, which has some features. Section V is devoted
to the calculation of the triple correlation function of vorticity. In Section
VI, we briefly expose results of direct numerical simulations confirming our
analytical analysis. In Conclusion, we outline the results of our work and
discuss some perspectives.

##  II General relations

In this study, we consider an unbounded two-dimensional fluid. We are
interested in fluctuations of the flow excited by an external, relatively
weak, random force, and examine effects related to their nonlinear
interaction. We take into account two dissipation mechanisms: bottom friction
and viscosity. The external force is assumed to be a random function of time
and coordinates possessing statistical properties homogeneous in time and
space.

The flow excited by the external force is described by a (two-dimensional)
velocity field $\bm{v}$, which is a random function of time and coordinates,
as well as of the external force. The flow is assumed to be incompressible,
$\nabla\bm{v}=0$. We examine the statistically stationary state of the fluid.
It is characterized by correlation functions of the velocity. For example, one
can be interested in the pair correlation function
$\langle\bm{v}(\bm{r})\bm{v}(\bm{x})\rangle$ for the velocities taken at the
points $\bm{r}$ and $\bm{x}$. The angular brackets here and below indicate
averaging over the system statistics.

Having in mind dynamics of thin films, one should add to the two-dimensional
Navier-Stokes equation the term related to bottom friction:

$\partial_{t}\bm{v}+(\bm{v}\nabla)\bm{v}+\nabla
p=-\alpha\bm{v}+\nu\nabla^{2}\bm{v}+\bm{f}.$ (1)

Here $p$ is pressure, $\alpha$ is the bottom friction coefficient and $\nu$ is
the kinematic viscosity coefficient. Taking the derivative of the equation
(1), we arrive at the following relation for pressure $p$:

$(\partial_{\alpha}v_{\beta})(\partial_{\beta}v_{\alpha})+\nabla^{2}p=\nabla\bm{f},$

where Greek indices run over $1,2$.

In two dimensions, it is convenient to describe the flow in terms of vorticity
$\varpi$:

$\varpi=\mathrm{curl}\,\bm{v}\equiv\partial_{1}v_{2}-\partial_{2}v_{1}.$ (2)

Obviously, $\varpi$ is a scalar (or, more precisely, a pseudoscalar) field.
The equation controlling the flow is derived from Eq. (1) by taking curl of
it:

$\partial_{t}\varpi+\bm{v}\nabla\varpi=-\alpha\varpi+\nu\nabla^{2}\varpi+\mathrm{curl}\,\bm{f}.$
(3)

Note that pressure drops from the equation.

To close the equation (3) one should restore the velocity field $\bm{v}$ from
the vorticity field $\varpi$. Due to the incompressibility condition
$\partial_{1}v_{1}+\partial_{2}v_{2}=0$ it is possible to introduce the stream
function $\Psi$, related to the velocity components and to the vorticity as

$v_{1}=\frac{\partial\Psi}{\partial x_{2}},\quad
v_{2}=-\frac{\partial\Psi}{\partial x_{1}},\quad\varpi=-\nabla^{2}\Psi.$ (4)

To determine the stream function, it is necessary to solve the Laplace
equation $\nabla^{2}\Psi=-\varpi$. Its solution can be written, say, as the
integral

$\Psi(\bm{r})=-\frac{1}{2\pi}\int d^{2}x\,\varpi(\bm{x})\ln|\bm{r}-\bm{x}|.$
(5)

After calculating the integral one finds the velocity components in accordance
with Eq. (4), thus expressing the velocity via the vorticity.

Multiplying the equation (1) by $\bm{v}$ and averaging it, one obtains the
energy balance

$\langle\bm{f}\bm{v}\rangle=\alpha\langle
v^{2}\rangle+\nu\langle\varpi^{2}\rangle.$ (6)

We introduce the special notation for the average in the left hand side of Eq.
(6), $\epsilon=\langle\bm{f}\bm{v}\rangle$. It is the average power of the
external force per unit mass. At deriving Eq. (6) we omitted all full
derivatives over time and coordinates, having in mind homogeneity in space and
time. The relation (6) has simple physical meaning: the energy pumped into the
fluid is dissipated by bottom friction and viscosity.

Analogously, multiplying the equation (3) by $\varpi$ and averaging, one
obtains

$\langle\phi\varpi\rangle=\alpha\langle\varpi^{2}\rangle+\nu\langle(\nabla\varpi)^{2}\rangle,$
(7)

where $\phi=\mathrm{curl}\,\bm{f}$. The physical meaning of the relation (7)
is the enstrophy balance. The enstrophy production rate
$\langle\phi\varpi\rangle$ can be estimated as $\epsilon k_{f}^{2}$ where
$k_{f}$ is the characteristic wave vector of the pumping. The relations (6)
and (7) are exact and are correct for any level of the pumping.

In the paper, we exploit the model where the pumping force $\bm{f}$ is shortly
correlated in time and has zero mean, $\langle\bm{f}\rangle=0$. Then
statistical properties of the external pumping are determined solely by the
pair correlation function of the force, that can be written as

$\langle
f_{\alpha}(t,\bm{x})f_{\beta}(\tau,\bm{y})\rangle=2\epsilon\delta_{\alpha\beta}\delta(t-\tau)\varXi(\bm{x}-\bm{y}).$
(8)

The expression implies homogeneity in space and time. Here $\varXi(\bm{r})$ is
a function of coordinates determined by details of pumping. To ensure the
property $\epsilon=\langle\bm{f}\bm{v}\rangle$ one should take
$\varXi(\bm{0})=1$. Passing to the quantity $\phi=\mathrm{curl}\,\bm{f}$,
figuring in the equation (3) for vorticity, one finds from Eq. (8)

$\langle\phi(t,\bm{x})\phi(\tau,\bm{y})\rangle=-2\epsilon\delta(t-\tau)\nabla^{2}\varXi(\bm{x}-\bm{y}),$
(9)

where the gradient $\nabla$ can be taken either over $\bm{x}$ or over
$\bm{y}$. Note that the enstrophy production rate per unit mass
$\langle\phi\varpi\rangle$ in our model is unambiguously determined by the
function $\varXi$,

$\langle\phi\varpi\rangle=-\epsilon\nabla^{2}\varXi(\bm{0}),$ (10)

as it follows from Eqs. (8) and (9).

We assume, that the pumping is statistically isotropic, i.e. $\varXi(\bm{r})$
depends solely on the absolute value of $\bm{r}$. Next, we assume that the
function $\varXi(\bm{r})$ has a characteristic length $k_{f}^{-1}$
(correlation length of the pumping), where $k_{f}$ is the characteristic wave
vector of the force $\bm{f}$. As an example, one can think about the
expression

$\varXi(\bm{r})=\exp\left(-k_{f}^{2}r^{2}/2\right),$ (11)

satisfying the condition $\varXi(\bm{0})=1$ and decaying fast at large $r$.
Its Fourier transform

$\tilde{\varXi}(\bm{k})=\frac{2\pi}{k_{f}^{2}}\exp\left(-\frac{k^{2}}{2k_{f}^{2}}\right),$
(12)

decays fast as the wave vector $k$ grows.

In the linear approximation, equation (3) is reduced to

$(\partial_{t}+\alpha-\nu\nabla^{2})\varpi=\phi.$ (13)

To find correlation functions of $\varpi$, one should solve the equation (13)
for arbitrary $\phi$ and then average the corresponding product over the
statistics of $\phi$, determined by Eq. (9). To find corrections to the linear
approximation, one has to take into account the nonlinear term
$\bm{v}\nabla\varpi$ in the equation (3). Corrections to the “bare”
correlation functions, found in the linear approximation, can be obtained as a
series over the nonlinear term.

Our approach implies that the size of the box $L$ is much larger than all
characteristic lengths of the problem. First of all, $L$ should be much larger
than the pumping length $k_{f}^{-1}$, $k_{f}L\gg 1$. However, some corrections
can be determined by the scales, much larger than the pumping length
$k_{f}^{-1}$, depending on the ratio $\nu k_{f}^{2}/\alpha$. The scales should
be much less than $L$ as well, for validity of our results.

One expects that the corrections are weak provided one of two dimensionless
parameters is small

$\beta_{\nu}=\frac{\epsilon}{\nu^{3}k_{f}^{4}},\qquad\beta_{\alpha}=\frac{\epsilon
k_{f}^{2}}{\alpha^{3}}.$ (14)

The first parameter $\beta_{\nu}$ in Eq. (14) is a power of Reynolds number
taken at the pumping scale. The second parameter $\beta_{\alpha}$ in Eq. (14)
is related to bottom friction and characterizes its strength at the same
pumping scale. Further we check the expectation, that the parameters (14)
control the perturbation series, and examine an interplay of viscosity and
bottom friction.

##  III Perturbation theory

The nonlinear interaction of the flow fluctuations can be consistently
examined in the framework of Wyld diagrammatic technique Wyld61 . The
diagrammatic technique can be derived from the representation of correlation
functions as functional integrals over the observed variables and auxiliary
fields MSR73 . The integration is performed like in the quantum field theory
qufield . A detailed description of the technique can be found in the review
HRS .

In our case, the functional integral can be written in terms of the vorticity
$\varpi$ and the corresponding auxiliary scalar field $\mu$. The functional
integral is taken with the weight $\exp(-{\mathcal{I}})$, where
${\mathcal{I}}$ is the effective action. The expression for the effective
action ${\mathcal{I}}$ is derived from the equation (3) and the relation (9).
It is the sum of two terms

$\displaystyle{\mathcal{I}}={\mathcal{I}}_{2}+{\mathcal{I}}_{int},$ (15)

where

$\displaystyle{\mathcal{I}}_{2}=\int
dt\,d^{2}x\,\mu(\partial_{t}+\alpha-\nu\nabla^{2})\varpi$
$\displaystyle+\epsilon\int
dt\,d^{2}x\,d^{2}r\,\nabla^{2}\varXi(\bm{x}-\bm{r})\mu(t,\bm{x})\mu(t,\bm{r}),$
(16) $\displaystyle{\mathcal{I}}_{int}=\int
dt\,d^{2}x\,\mu\bm{v}\nabla\varpi.$ (17)

The velocity $\bm{v}$ in Eq. (17) is implied to be expressed via the vorticity
$\varpi$, see Section II.

The pair correlation function of the vorticity can be written as the following
functional integral

$\langle\varpi(t,\bm{x})\varpi(0,\bm{y})\rangle=\int
D\varpi\,D\mu\,e^{-{\mathcal{I}}}\varpi(t,\bm{x})\varpi(0,\bm{y}).$ (18)

We also introduce the Green function

$\langle\varpi(t,\bm{x})\mu(0,\bm{y})\rangle=\int
D\varpi\,D\mu\,e^{-{\mathcal{I}}}\varpi(t,\bm{x})\mu(0,\bm{y}),$ (19)

which determines the response of the system to an additional external force
added to Eq. (1). Analogously to Eqs. (18) and (19), any correlation function
of the fields $\varpi,\mu$ can be written. Note that the average
$\langle\mu\mu\rangle$ is zero.

Vorticity (2) changes its sign under reflection. Reflection relative to the
second coordinate axis implies the following transformation

$\displaystyle x_{1}\to-x_{1},\quad x_{2}\to x_{2},\quad v_{1}\to-v_{1},$
$\displaystyle v_{2}\to v_{2},\quad\varpi\to-\varpi,\quad\mu\to-\mu.$ (20)

The effective action (15) is invariant under the transformation (20).
Therefore all odd moments of vorticity, including the third moment, are zero.
However, the triple correlation function of vorticity is nonzero, since it
changes its sign together with changing sign of the first components of the
coordinates.

Neglecting the third order term ${\mathcal{I}}_{int}$ (17) in the functional
integrals (18) and (19), we arrive at the Gaussian functional integrals
determining the “bare” correlation functions:

$\displaystyle\langle\varpi(t,\bm{x})\mu(0,\bm{y})\rangle_{0}=\int{\mathcal{D}}\varpi{\mathcal{D}}\mu\,e^{-{\mathcal{I}}_{2}}\varpi(t,\bm{x})\mu(0,\bm{y}),$
(21)
$\displaystyle\langle\varpi(t,\bm{x})\varpi(0,\bm{y})\rangle_{0}=\int{\mathcal{D}}\varpi{\mathcal{D}}\mu\,e^{-{\mathcal{I}}_{2}}\varpi(t,\bm{x})\varpi(0,\bm{y}),$
(22)

corresponding to the linear approximation in the hydrodynamic equation. The
averages (21) and (22) can be found explicitly.

The expressions (16) and (21) lead to the following equation for the “bare”
Green function

$(\partial_{t}+\alpha-\nu\nabla^{2})\langle\varpi(t,\bm{x})\mu(0,\bm{y})\rangle_{0}=\delta(t)\delta(\bm{x}-\bm{y}).$
(23)

Remind that any Green function is zero at negative times due to causality.
Therefore $\langle\varpi(t,\bm{x})\mu(0,\bm{y})\rangle$ is zero at $t<0$. The
expression for the “bare” pair correlation function can be derived using Eq.
(22):

$\displaystyle\langle\varpi(t,\bm{x})\varpi(0,\bm{y})\rangle_{0}=-2\epsilon\int
d\tau\int d^{2}rd^{2}z\nabla^{2}\varXi(\bm{r}-\bm{z})$
$\displaystyle\langle\varpi(t,\bm{x})\mu(\tau,\bm{r})\rangle_{0}\langle\varpi(0,\bm{y})\mu(\tau,\bm{z})\rangle_{0}\quad$
(24)

in real space.

Let us pass to Fourier transforms of the pair correlation function (24) and of
the Green function (23):

$\displaystyle\langle\varpi(t,\bm{x})\mu(0,\bm{y})\rangle_{0}=\int\frac{d^{2}k}{(2\pi)^{2}}e^{i\bm{k}(\bm{x}-\bm{y})}{G}(t,\bm{k}),$
(25)
$\displaystyle\langle\varpi(t,\bm{x})\varpi(0,\bm{y})\rangle_{0}=\int\frac{d^{2}k}{(2\pi)^{2}}e^{i\bm{k}(\bm{x}-\bm{y})}{F}(t,\bm{k}).$
(26)

For the Fourier transform ${G}(t,\bm{k})$ (25) we derive from Eq. (23) the
following differential equation

$\displaystyle\left(\frac{\partial}{\partial
t}+\alpha+\nu\bm{k}^{2}\right){G}(t,\bm{q})=\delta(t).$ (27)

It can be easily solved to obtain

$\displaystyle G(t,\bm{k})=\theta(t)\exp(-\alpha t-\nu\bm{k}^{2}t),$ (28)

where $\theta(t)$ is the Heaviside step function, reflecting causality. Next,
we derive from Eqs. (24), (25) and (26)

$\displaystyle{F}(t,\bm{k})=2\epsilon\int\limits_{-\infty}^{\min(t,0)}d\tau\,k^{2}\tilde{\varXi}(\bm{k}){G}(t-\tau,\bm{k}){G}(-\tau,\bm{k}),$
(29)

where

$\tilde{\varXi}(\bm{k})=\int d^{2}x\,\exp(-i\bm{k}\bm{x})\varXi(\bm{x})$ (30)

is Fourier transform of $\varXi(\bm{r})$. Therefore, the simultaneous pair
correlation function is

$\displaystyle F(\bm{k})=\frac{\epsilon k^{2}}{\nu
k^{2}+\alpha}\tilde{\varXi}(\bm{k})$ (31)

and

$F(t,\bm{k})=G(|t|,\bm{k})F(\bm{k}).$ (32)

The last expression is a consequence of the short correlation in the time of
our pumping.

The “bare” second moment of vorticity is expressed as

$\displaystyle\langle\varpi^{2}\rangle_{0}=\int\frac{d^{2}k}{(2\pi)^{2}}F(\bm{k})=\int\frac{d^{2}k}{(2\pi)^{2}}\frac{\epsilon
k^{2}}{\alpha+\nu k^{2}}\tilde{\varXi}(\bm{k}).$ (33)

If $\alpha\sim\nu k_{f}^{2}$ then one derives from Eq. (33)

$\langle\varpi^{2}\rangle_{0}\sim\epsilon/\nu\sim\epsilon k_{f}^{2}/\alpha.$
(34)

If $\alpha\ll\nu k_{f}^{2}$ then $\langle\varpi^{2}\rangle_{0}=\epsilon/\nu$
and in the opposite case $\alpha\gg\nu k_{f}^{2}$ one finds the estimate
$\langle\varpi^{2}\rangle_{0}\sim\epsilon k_{f}^{2}/\alpha$. All the values
are proportional to the power (energy production rate) per unit mass
$\epsilon$, with different factors.

The bare value of the second moment of the velocity can be found as

$\langle
v^{2}\rangle_{0}=\int\frac{d^{2}k}{(2\pi)^{2}}\frac{F(k)}{k^{2}}=\epsilon\int\frac{d^{2}k}{(2\pi)^{2}}\frac{\tilde{\varXi}(k)}{\alpha+\nu
k^{2}}.$ (35)

In the limit $\alpha\gg\nu k_{f}^{2}$ the expression (35) gives $\langle
v^{2}\rangle_{0}=\epsilon/\alpha$. In the opposite limit $\alpha\ll\nu
k_{f}^{2}$ the integral (35) is gained from the interval
$\sqrt{\alpha/\nu}<k<k_{f}$. Then one obtains

$\langle
v^{2}\rangle_{0}=\frac{\epsilon\tilde{\varXi}(0)}{4\pi\nu}\ln\frac{\nu
k_{f}^{2}}{\alpha}.$ (36)

The expression implies, that $\tilde{\Xi}(k)$ can be approximated by
$\tilde{\Xi}(0)$ inside the interval of the integration. Note that $\langle
v^{2}\rangle_{0}\sim k_{f}^{-2}\langle\varpi^{2}\rangle_{0}$, up to the
logarithmic factor at small $\alpha$.

One can also calculate the turbulent diffusion coefficient

$\displaystyle D=\int dt\,\langle\bm{v}(t,\bm{r})\bm{v}(0,\bm{r})\rangle_{0}$
$\displaystyle=\frac{\epsilon}{2}\int\frac{d^{2}k}{(2\pi)^{2}}\frac{\tilde{\varXi}(k)}{(\alpha+\nu
k^{2})^{2}}.$ (37)

It determines an evolution of an impurity field if molecular diffusion is
smaller than $D$. In the limit $\alpha\gg\nu k_{f}^{2}$, expression (37) gives

$D=\epsilon/(2\alpha^{2}).$

In the opposite limit $\alpha\ll\nu k_{f}^{2}$, we find

$D=\frac{\epsilon\tilde{\varXi}(0)}{16\alpha\nu}\sim\frac{\epsilon}{\alpha\nu
k_{f}^{2}}.$

Note that $\alpha$ does not drop from the expression since it is determined by
$k\sim\sqrt{\alpha/\nu}$.

To conclude this section, we would like to briefly examine the scenario when
the pumping has a narrow spectrum $\tilde{\varXi}(\bm{k})$ in the Fourier
space. One can consider the infinitely thin function

$\tilde{\varXi}(\bm{k})=2\pi\epsilon\delta(k-k_{f})/k_{f},$ (38)

satisfying the condition $\varXi(\bm{0})=1$, in the real space. Then the
integrals over $\bm{k}$ can be easily calculated, and we obtain

$\displaystyle\langle\varpi^{2}\rangle_{0}=\frac{\epsilon
k_{f}^{2}}{\alpha+\nu k_{f}^{2}},\;\langle
v^{2}\rangle_{0}=\frac{\epsilon}{\alpha+\nu
k_{f}^{2}},\;D=\frac{\epsilon}{2(\alpha+\nu k_{f}^{2})^{2}}.$

Note that the expressions for $\langle v^{2}\rangle_{0}$ and $D$ in the limit
$\alpha\ll\nu k_{f}^{2}$ are different from ones discussed above since
$\tilde{\varXi}(\bm{0})=0$ for the spectrum (38).

##  IV First corrections

Corrections to the bare values of the correlation functions are given by the
perturbation series. They are obtained by the expansion of the factor
$\exp({\mathcal{I}}_{2}+{\mathcal{I}}_{int})$ in the series over
${\mathcal{I}}_{int}$ in the corresponding functional integral. Each term of
the expansion can be found analytically, using Wick theorem Wick50 : the
average of a product of the fields $\varpi,\mu$ is equal to the sum of
products of the bare pair averages (21) and (22) organized by all possible
pairings. As a result, we come to multiple integrals over times and wave
vectors of some expressions that are products of (28) and (32). One should
also introduce factors, corresponding to gradients in $\nabla\varpi$ and to
converting $\bm{v}\to\varpi$.

$t$$0$$\tau$$s$$\bm{k}$$\bm{q}+\bm{k}$$\bm{q}$ Figure 1: One-loop diagram
representing the first correction to the Green function $G(t,\bm{k})$.

The terms of the expansion can be represented by Feynman diagrams of different
orders, see Ref. qufield . Let us give some examples. The diagram determining
the first order correction to the Green function is depicted in Fig. 1. One of
the diagrams determining the next order correction to the Green function is
depicted in Fig. 2. Higher order corrections correspond to more complicated
diagrams. Note that the first order corrections correspond to one-loop
diagrams, the second order corrections correspond to two-loop diagrams and so
further.

All lines in the diagrams consist of two segments: solid segments correspond
to the field $\varpi$, whereas dashed segments correspond to the field $\mu$.
Thus, a combined solid-dashed line designates the Green function
${\mathcal{G}}(t,\bm{k})$ and a solid line designate the pair correlation
function ${\mathcal{F}}(t,\bm{k})$. Two solid segments and one dashed segment
are attached to each vertex in agreement with the structure of
${\mathcal{I}}_{int}$, see expression (17). The factor corresponding to each
vertex is

$\frac{1}{2}\left(\frac{1}{q^{2}}-\frac{1}{k^{2}}\right)(q_{2}k_{1}-q_{1}k_{2}),$
(39)

where $\bm{k},\bm{q}$ are wave vectors of the solid segments attached to the
vertex.

To construct an analytical expression, corresponding to a given diagram, one
should fix the corresponding combinatorial factor, take the product of Green
functions, pair correlation functions, vertex factors and integrate the result
over “internal” wave vectors and times. The integration should be performed at
the condition of the wave vector conservation at each vertex: the sum of the
wave vectors of three segments attached to the vertex has to be zero.

Figure 2: Two-loop diagram representing a correction to the Green function.

###  IV.1 Correction to the pair correlation function

Here we examine the first correction $\delta F(\bm{k})$ to the simultaneous
pair correlation function $F(\bm{k})$. The correction contains two
contributions, determined by the diagrams depicted in Fig. 3, labelled by the
letters a and b. The diagrams correspond to the following analytical
expressions

$\displaystyle
F_{\mathrm{a}}(\bm{k})=\frac{1}{2}\int^{\infty}_{0}dt\int_{-t}^{\infty}ds\int\frac{d^{2}m}{(2\pi)^{2}}$
$\displaystyle
G(t,\bm{k})G(|s|,\bm{p})G(|s|,\bm{m})G(s+t,\bm{k})F(\bm{m})F(\bm{p})$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)^{2}(m_{2}k_{1}-m_{1}k_{2})^{2},$
(40)

and

$\displaystyle
F_{\mathrm{b}}(\bm{k})=-2\int_{-\infty}^{0}dt\int_{0}^{\infty}ds\int\frac{d^{2}m}{(2\pi)^{2}}$
$\displaystyle
G(t,\bm{k})G(s,\bm{p})G(s,\bm{m})G(s+t,\bm{k})F(\bm{m})F(\bm{k})$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)\left(\frac{1}{\bm{m}^{2}}-\frac{1}{{\bm{k}}^{2}}\right)(m_{2}k_{1}-m_{1}k_{2})^{2},$
(41)

where $\bm{p}=-\bm{k}-\bm{m}$ and we substituted expression (32).

a$0$$0$$-t$$-t-s$$-\bm{k}$$\bm{k}$$\bm{p}$$\bm{m}$.

b$0$$0$$-t$$-t-s$$-\bm{k}$$\bm{k}$$\bm{p}$$\bm{m}$

Figure 3: First-order contributions to the simultaneous pair correlation
function.

Using expressions (28) and (31), one can take the integrals over the times
$t,s$ in Eqs. (40) and (41) to obtain

$\displaystyle F_{\mathrm{a}}(\bm{k})=\frac{\epsilon^{2}}{2(\alpha+\nu
k^{2})}\int\frac{d^{2}m}{(2\pi)^{2}}\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{p})$
$\displaystyle\frac{m^{2}p^{2}}{(\alpha+\nu m^{2})(\alpha+\nu
p^{2})(3\alpha+\nu k^{2}+\nu p^{2}+\nu m^{2})}$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)^{2}(m_{2}k_{1}-m_{1}k_{2})^{2},$
(42)

and

$\displaystyle F_{\mathrm{b}}(\bm{k})=-\frac{\epsilon^{2}}{\alpha+\nu
k^{2}}\int\frac{d^{2}m}{(2\pi)^{2}}\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{k})$
$\displaystyle\frac{m^{2}k^{2}}{(\alpha+\nu m^{2})(\alpha+\nu
k^{2})(3\alpha+\nu k^{2}+\nu p^{2}+\nu m^{2})}$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)\left(\frac{1}{m^{2}}-\frac{1}{{k}^{2}}\right)(m_{2}k_{1}-m_{1}k_{2})^{2}.$
(43)

The general expressions (42) and (43) enable one to analyze the correction
$\delta F=F_{\mathrm{a}}+F_{\mathrm{b}}$ to the pair correlation function and
to establish its behavior in different limiting cases.

Let us integrate over $\bm{k}$ the product $(\alpha+\nu k^{2})\delta
F(\bm{k})$. One finds from Eq. (42)

$\displaystyle\int\frac{d^{2}k}{(2\pi)^{2}}(\alpha+\nu
k^{2})F_{\mathrm{a}}(\bm{k})=\epsilon^{2}\int\frac{d^{2}p\,d^{2}m}{(2\pi)^{4}}\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{p})$
$\displaystyle\frac{p^{2}}{(\alpha+\nu m^{2})(\alpha+\nu p^{2})(3\alpha+\nu
k^{2}+\nu p^{2}+\nu m^{2})}$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)(m_{2}p_{1}-m_{1}p_{2})^{2},\qquad$
(44)

where we passed to the integration over $\bm{p}$ and used symmetric properties
of the integrand. The second integral is equal to

$\displaystyle\int\frac{d^{2}k}{(2\pi)^{2}}(\alpha+\nu
k^{2})F_{\mathrm{b}}(\bm{k})=-\epsilon^{2}\int\frac{d^{2}k\,d^{2}m}{(2\pi)^{4}}\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{k})$
$\displaystyle\frac{1}{(\alpha+\nu m^{2})(\alpha+\nu k^{2})(3\alpha+\nu
k^{2}+\nu p^{2}+\nu m^{2})}$
$\displaystyle\left[k^{2}-\frac{k^{2}m^{2}}{p^{2}}\right]\left(\frac{1}{m^{2}}-\frac{1}{{k}^{2}}\right)(m_{2}k_{1}-m_{1}k_{2})^{2}.$

Here the first contribution, related to the term $k^{2}$ in the square
brackets, is cancelled by the integral (44), and the second contribution,
related to the term $k^{2}m^{2}/p^{2}$, is zero since the integrand changes
its sign under the permutation $\bm{m}\leftrightarrow\bm{k}$. We conclude that

$\int\frac{d^{2}k}{(2\pi)^{2}}(\alpha+\nu k^{2})\delta F(\bm{k})=0.$ (45)

The relation is a direct consequence of the enstrophy balance (7), since
integral (45) is the correction of the order of $\epsilon^{2}$ to the right
hand side of relation (7), and expression (10) shows that such correction
should be absent. Thus, our calculations are in agreement with the general
relation.

Another check can be extracted from the energy balance (6). It demonstrates
that the correction of the order of $\epsilon^{2}$ to the right-hand side of
relation (6) should be absent

$\int\frac{d^{2}k}{(2\pi)^{2}}\left(\frac{\alpha}{k^{2}}+\nu\right)\delta
F(\bm{k})=0.$ (46)

Indeed,

$\displaystyle\int\frac{d^{2}k}{(2\pi)^{2}}\left(\frac{\alpha}{k^{2}}+\nu\right)F_{\mathrm{a}}(\bm{k})=\epsilon^{2}\int\frac{d^{2}p\,d^{2}m}{(2\pi)^{4}}\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{p})$
$\displaystyle\frac{(-m^{2}/k^{2})}{(\alpha+\nu m^{2})(\alpha+\nu
p^{2})(3\alpha+\nu k^{2}+\nu p^{2}+\nu m^{2})}$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)(m_{2}p_{1}-m_{1}p_{2})^{2},\qquad$
(47)

and

$\displaystyle\int\frac{d^{2}k}{(2\pi)^{2}}\left(\frac{\alpha}{k^{2}}+\nu\right)F_{\mathrm{b}}(\bm{k})=-\epsilon^{2}\int\frac{d^{2}k\,d^{2}m}{(2\pi)^{4}}\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{k})$
$\displaystyle\frac{1}{(\alpha+\nu m^{2})(\alpha+\nu k^{2})(3\alpha+\nu
k^{2}+\nu p^{2}+\nu m^{2})}$
$\displaystyle\left[1-\frac{m^{2}}{p^{2}}\right]\left(\frac{1}{m^{2}}-\frac{1}{{k}^{2}}\right)(m_{2}k_{1}-m_{1}k_{2})^{2}.$

Here the first contribution, related to unity in the square brackets, is zero
since the integrand changes its sign under the permutation
$\bm{m}\leftrightarrow\bm{k}$, and the second contribution, related to the
term $m^{2}/p^{2}$ in the square brackets, is cancelled by Eq. (47).

Next, let us consider the limiting case $\alpha\ll\nu k_{f}^{2}$. Putting
$\alpha\to 0$ in the expressions (42) and (43), we obtain

$\displaystyle
F_{\mathrm{a}}(\bm{k})=\frac{\epsilon^{2}}{2\nu^{4}k^{2}}\int\frac{d^{2}m}{(2\pi)^{2}}\frac{\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{p})}{k^{2}+p^{2}+m^{2}}$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)^{2}(m_{2}k_{1}-m_{1}k_{2})^{2},$
(48)

and

$\displaystyle
F_{\mathrm{b}}(\bm{k})=-\frac{\epsilon^{2}}{\nu^{4}k^{2}}\int\frac{d^{2}m}{(2\pi)^{2}}\frac{\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{k})}{k^{2}+p^{2}+m^{2}}$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)\left(\frac{1}{\bm{m}^{2}}-\frac{1}{{\bm{k}}^{2}}\right)(m_{2}k_{1}-m_{1}k_{2})^{2}.$
(49)

One can rewrite the sum $F_{\mathrm{a}}(\bm{k})+F_{\mathrm{b}}(\bm{k})$ as
$A+B$, where

$\displaystyle
A=\frac{\epsilon^{2}}{\nu^{4}k^{2}}\int\frac{d^{2}m}{(2\pi)^{2}}\frac{\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{p})}{k^{2}+p^{2}+m^{2}}(m_{2}k_{1}-m_{1}k_{2})^{2}$
$\displaystyle\left[-\frac{1}{2m^{4}}+\frac{1}{2p^{4}}+\frac{1}{k^{2}m^{2}}-\frac{1}{k^{2}p^{2}}\right],$
(50)

and

$\displaystyle
B=-\frac{\epsilon^{2}}{\nu^{4}k^{2}}\int\frac{d^{2}m}{(2\pi)^{2}}\frac{\tilde{\varXi}(\bm{m})[\tilde{\varXi}(\bm{k})-\tilde{\varXi}(\bm{p})]}{k^{2}+p^{2}+m^{2}}$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)\left(\frac{1}{m^{2}}-\frac{1}{{k}^{2}}\right)(m_{2}k_{1}-m_{1}k_{2})^{2}.$
(51)

The quantity $A$ is equal to zero, it can be checked performing the change of
the integration variable $\bm{m}\to\bm{p}$ in Eq. (50). Thus, the correction
$\delta F$ is reduced to the value of $B$.

Figure 4: The integral (55) as a function of $k/k_{f}$.

There are no divergencies in integral (51). Therefore, we conclude that at
$k\sim k_{f}$

$\delta
F(\bm{k})=B\sim\frac{\epsilon^{2}}{\nu^{4}k_{f}^{6}}\sim\frac{\epsilon}{\nu^{3}k_{f}^{4}}F,$
(52)

where we used expression (31). Thus, the small parameter of the perturbation
series is the parameter $\beta_{\nu}$ (14). Note that in the limit $\alpha\to
0$ the relations (45) and (46) are reduced to

$\int d^{2}k\,B(\bm{k})=0,\quad\int d^{2}k\,k^{2}B(\bm{k})=0.$ (53)

These relations can be checked directly, using the symmetry properties of
integral (51). In accordance with Eq. (53), the function $B(k)$ should change
its sign twice when varying $k$.

Figure 5: The integral (61) as a function of $k/k_{f}$.

It is instructive to take a particular function $\varXi$, enabling to
calculate integrals for the corrections up to a number. We will exploit
function (11), which has Fourier transform (12). Then expression (51) is
reduced to

$\displaystyle\delta F(\bm{k})=B=-\frac{\epsilon^{2}}{\nu^{4}k_{f}^{6}}I,$
(54) $\displaystyle
I=\int_{0}^{\infty}\frac{dm}{m}\int_{-1}^{1}dx\sqrt{1-x^{2}}$
$\displaystyle\frac{k_{f}^{2}}{k^{2}+m^{2}-kmx}\exp\left(-\frac{m^{2}+k^{2}}{2k_{f}^{2}}\right)$
$\displaystyle\left[1-\exp\left(-\frac{m^{2}-2kmx}{2k_{f}^{2}}\right)\right]$
$\displaystyle\left(1-\frac{m^{2}}{k^{2}+m^{2}-2kmx}\right)\left(1-\frac{m^{2}}{{k}^{2}}\right).$
(55)

The integral converges and is, consequently, of the order of unity for $k\sim
k_{f}$. Thus, expression (54) confirms the general evaluation (52).

The integral (55) (found numerically) is plotted in Fig. 4 as a function of
$k/k_{f}$. In accordance with the general properties (53), the function
changes its sign twice when varying $k/k_{f}$. The integral has the following
asymptotics

$\displaystyle I\approx-\frac{\pi}{8}\ln(k_{f}/k),\quad k\ll k_{f},$
$\displaystyle I\approx-\frac{2^{9}\pi
k_{f}^{8}}{k^{8}}\exp\left(-\frac{k^{2}}{4k_{f}^{2}}\right),\quad k\gg k_{f},$

which can be found from Eq. (55). Note that the argument of the exponent is
two times smaller than in Eq. (12). It is a consequence of the fact that
$\delta F$ is formed by nonlinear confluence of two fluctuations.

Next, we consider the opposite limiting case $\alpha\gg\nu k_{f}^{2}$. Then
expressions (42) and (43) are reduced to

$\displaystyle
F_{\mathrm{a}}(\bm{k})=\frac{\epsilon^{2}}{6\alpha^{4}}\int\frac{d^{2}m}{(2\pi)^{2}}\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{p})$
$\displaystyle(p^{2}-m^{2})\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)(m_{2}k_{1}-m_{1}k_{2})^{2},$
(56)

and

$\displaystyle
F_{\mathrm{b}}(\bm{k})=-\frac{\epsilon^{2}}{3\alpha^{4}}\int\frac{d^{2}m}{(2\pi)^{2}}\tilde{\varXi}(\bm{m})\tilde{\varXi}(\bm{k})$
$\displaystyle\left(\frac{1}{m^{2}}-\frac{1}{p^{2}}\right)(k^{2}-m^{2})(m_{2}k_{1}-m_{1}k_{2})^{2}.$
(57)

Assuming $k\sim m\sim k_{f}$, we arrive at the estimate

$F_{\mathrm{a}}\sim
F_{\mathrm{b}}\sim\frac{\epsilon^{2}k_{f}^{2}}{\alpha^{4}}\sim\frac{\epsilon
k_{f}^{2}}{\alpha^{3}}F,$ (58)

where we used expression (31). Thus, the small parameter of the perturbation
series in the limiting case $\alpha\gg\nu k_{f}^{2}$ is the parameter
$\beta_{\alpha}$ (14).

Note that in the limit $\nu\to 0$, relations (45) and (46) are reduced to

$\int d^{2}k\,\delta F(\bm{k})=0,\quad\int d^{2}k\,k^{-2}\delta F(\bm{k})=0.$
(59)

These relations can be checked directly, using the symmetry properties of
expressions (56) and (57). In accordance with Eq. (59), the function $\delta
F(k)$ should change its sign twice when varying $k$.

Figure 6: Dependence of functions $C_{\eta}$ on parameter $k/k_{f}$ for
different values of $\eta=\nu k_{f}^{2}/\alpha$.

Again, it is instructive to calculate the correction for a particular function
$\varXi$. Using expression (12), one finds

$F_{\mathrm{a}}+F_{\mathrm{b}}=-\frac{\epsilon^{2}k_{f}^{2}}{3\alpha^{4}}J,$
(60)

where

$\displaystyle
J=\frac{k^{2}}{k_{f}^{6}}\int_{0}^{\infty}dm\,m\int_{-1}^{1}dx\,\sqrt{1-x^{2}}$
$\displaystyle\exp\left(-\frac{k^{2}+m^{2}}{2k_{f}^{2}}\right)\left\\{2\frac{(k^{2}-m^{2})(k^{2}-2kmx)}{k^{2}+m^{2}-2kmx}\right.$
$\displaystyle\left.-\frac{(k^{2}-2kmx)^{2}}{k^{2}+m^{2}-2kmx}\exp\left(-\frac{m^{2}}{2k_{f}^{2}}+\frac{mkx}{k_{f}^{2}}\right)\right\\}.$
(61)

The integral converges, and thus $J\sim 1$ for $k\sim k_{f}$, confirming the
general evaluation (58).

The integral (61) (found numerically) is plotted in Fig. 5 as a function of
$k/k_{f}$. In accordance with the general properties (59), the function
changes its sign twice when varying $k/k_{f}$. The asymptotics of the integral
(61)

$\displaystyle J\approx-\frac{\pi}{4}\frac{k^{4}}{k_{f}^{4}},\quad k\ll
k_{f},$ $\displaystyle
J\approx-8\pi\exp\left(-\frac{k^{2}}{4k_{f}^{2}}\right),\quad k\gg k_{f},$

can be found directly from expression (61).

To handle an arbitrary ratio between parameters $\alpha$ and $\nu k_{f}^{2}$,
based on expressions (52) and (58), we introduce the following expression for
the correction

$\delta
F(\bm{k})=\dfrac{\beta_{\alpha}\beta_{\nu}}{\beta_{\alpha}+\beta_{\nu}}F(\bm{k})C_{\eta}(\bm{k}),$

where $\eta=\nu k_{f}^{2}/\alpha$. The function $C_{\eta}(\bm{k})$ can be
found numerically using the expressions (42-43) and (31). The results for
forcing with the spectrum (12) are shown in Fig. 6. It can be concluded that
the values of $C_{\eta}$ are numerically small, so the perturbation theory can
be expected to work up to values of $\beta_{\nu}$ or $\beta_{\alpha}$ of the
order of $10$.

Some words about the case of a narrow pumping spectrum $\varXi(\bm{k})$ in the
Fourier space. In the extreme case (38), both contributions $F_{a}$ and
$F_{b}$ are equal to zero. When the spectrum has a finite width $\delta_{f}\ll
k_{f}$, corrections $F_{a}$ and $F_{b}$ are no longer equal to zero, but
contain an additional smallness compared to the estimates discussed earlier in
this section, which is controlled by the ratio $\delta_{f}/k_{f}\ll 1$. It is
interesting, that the contribution $F_{b}$ is localized near $k_{f}$, while
the contribution $F_{a}$ determines the correction $\delta F(k)$ in the region
$|k-k_{f}|\gg\delta_{f}$.

###  IV.2 Correction to the second moment of vorticity

The second moment of vorticity $\langle\varpi^{2}\rangle$ is the integral over
the wave vector of its simultaneous pair correlation function. Therefore the
first correction to the moment is

$\delta\langle\varpi^{2}\rangle=\int\frac{d^{2}k}{(2\pi)^{2}}(F_{\mathrm{a}}+F_{\mathrm{b}}),$

where $F_{\mathrm{a}}$ and $F_{\mathrm{b}}$ are corrections to the pair
correlation function (42) and (43). If $\alpha\sim\nu k_{f}^{2}$ then we find
the estimate

$\delta\langle\varpi^{2}\rangle\sim\frac{\epsilon^{2}}{\nu^{4}k_{f}^{4}}\sim\frac{\epsilon^{2}k_{f}^{4}}{\alpha^{4}},$
(62)

following from Eqs. (52) and (58). However, if we are interested in the
limiting cases $\nu k_{f}^{2}\gg\alpha$ or $\nu k_{f}^{2}\ll\alpha$, we cannot
use expression (51) or expressions (56) and (57) for the corrections to the
pair correlation function since the integrals of the functions are equal to
zero, see Eqs. (53) and (59).

To calculate correction to the second moment, we will use general expressions
(42) and (43). In the expression for $F_{\mathrm{a}}$, we pass from the
integration over $\bm{k}$ to the integration over $\bm{p}$. In the expression
for $F_{\mathrm{b}}$, we first symmetrize the integrand in $\bm{k},\bm{m}$ and
then permute $\bm{k}\leftrightarrow\bm{p}$. Collecting all terms, one finds

$\displaystyle\delta\langle\varpi^{2}\rangle=-\frac{\epsilon^{2}\alpha}{2\nu^{5}}\int\frac{d^{2}p\,d^{2}m}{(2\pi)^{4}}$
$\displaystyle\frac{\tilde{\varXi}(\bm{p})}{(p^{2}+\alpha/\nu)^{2}}\
\frac{\tilde{\varXi}(\bm{m})}{(m^{2}+\alpha/\nu)^{2}}\
(p_{1}m_{2}-p_{2}m_{1})^{2}$
$\displaystyle\frac{(m^{2}-p^{2})^{2}}{m^{2}p^{2}(\bm{m}+\bm{p})^{2}[\alpha/\nu+(\bm{m}+\bm{p})^{2}]}$
$\displaystyle\frac{[(\bm{m}+\bm{p})^{2}-m^{2}][(\bm{m}+\bm{p})^{2}-p^{2}]}{3\alpha/\nu+p^{2}+m^{2}+(\bm{m}+\bm{p})^{2}}.$
(63)

The expression is correct for an arbitrary function $\tilde{\varXi}$.

Now, let us consider the limiting case $\alpha\ll\nu k_{f}^{2}$. There is
small factor $\alpha/(\nu k_{f}^{2})$ in integral (63), but one should be
careful in estimating the correction. At small $\alpha$, there is an infrared
contribution to integral (63) coming from $p,m\sim\sqrt{\alpha/\nu}$. The
contribution behaves $\propto\alpha^{-1}$ and, consequently, can compensate
the small factor $\alpha$ in the integral. Since the infrared contribution to
integral (63) comes from $p,m\sim\sqrt{\alpha/\nu}\ll k_{f}$, we can
substitute factors $\tilde{\varXi}$ by constants $\tilde{\varXi}(\bm{0})$.
Then the infrared contribution to integral (63) appears to be proportional to

$\displaystyle\int d^{2}m\,d^{2}p\,d^{2}q\,\delta(\bm{m}+\bm{p}+\bm{q})$
$\displaystyle\frac{(p_{1}m_{2}-p_{2}m_{1})^{2}}{m^{2}p^{2}q^{2}}\frac{1}{3\alpha/\nu+p^{2}+m^{2}+q^{2}}$
$\displaystyle\left[\frac{a}{hc}-\frac{2}{c}+\frac{h}{ac}-\frac{a}{h^{2}}+\frac{1}{h}\right.$
$\displaystyle\left.+\frac{1}{a}-\frac{h}{a^{2}}+\frac{c}{h^{2}}-\frac{2c}{ah}+\frac{c}{a^{2}}\right],$

where $a=\alpha/\nu+p^{2}$, $h=\alpha/\nu+m^{2}$, $c=\alpha/\nu+q^{2}$. The
first factor in the integrand is invariant under permutations of
$\bm{p},\bm{m},\bm{q}$. Such permutation means the corresponding permutation
of $a,h,c$. Performing the permutations we find that the second factor in the
integrand (in the square brackets) is cancelled to zero. We conclude that the
leading infrared contribution to integral (63) is zero.

However, there is a subleading infrared contribution to integral (63)
determining the principal value of the correction to
$\langle\varpi^{2}\rangle$ due to the logarithmic divergence. To extract the
principal value, one has to use the next term of the expansion
$\tilde{\varXi}$ in the wave vectors:

$\tilde{\varXi}(p)\tilde{\varXi}(m)\to-\frac{2\pi^{2}C}{k_{f}^{6}}(p^{2}+m^{2}).$
(64)

Here $C$ is the constant of the order of unity, it is equal to unity $C=1$ for
the particular function (12). Substituting expression (64) into integral (63),
putting $\alpha\to 0$ in the integrand and taking the integral over the
angles, one finds

$\displaystyle\delta\langle\varpi^{2}\rangle=\frac{C\epsilon^{2}\alpha}{4\nu^{5}k_{f}^{6}}\int_{0}^{\infty}\frac{dp\,dm}{pm}\frac{(p^{2}-m^{2})^{2}}{(p^{2}+m^{2})^{2}}$
$\displaystyle\frac{1}{b^{4}}\left\\{12+\frac{b^{4}+4b^{2}-32}{\sqrt{4-b^{2}}}-\frac{b^{4}+2b^{2}-4}{\sqrt{1-b^{2}}}\right\\},$

where $b=2pm(p^{2}+m^{2})^{-1}$. Passing to the variables $b$,
$S=p^{2}+m^{2}$, we find

$\displaystyle\delta\langle\varpi^{2}\rangle=\frac{C\epsilon^{2}\alpha}{4\nu^{5}k_{f}^{6}}\int_{0}^{\infty}dS\int_{0}^{1}db\frac{\sqrt{1-b^{2}}}{Sb}$
$\displaystyle\frac{1}{b^{4}}\left\\{12+\frac{b^{4}+4b^{2}-32}{\sqrt{4-b^{2}}}-\frac{b^{4}+2b^{2}-4}{\sqrt{1-b^{2}}}\right\\}.$
(65)

The integral diverges logarithmically at small $S$ and $b$.

The logarithmic divergency is cut from below by the omitted terms with
$\alpha$ in the integrand and from above by the wave vector $k_{f}$. Thus we
arrive at the following regions of integration in Eq. (65)

$\frac{\alpha}{\nu}<S<k_{f}^{2},\quad\left(\frac{\alpha}{\nu
S}\right)^{1/2}<b<1.$

Keeping the main logarithmic term, one finds

$\displaystyle\delta\langle\varpi^{2}\rangle=-\frac{C\epsilon^{2}\alpha}{2^{10}\nu^{5}k_{f}^{6}}\left(\ln\frac{\nu
k_{f}^{2}}{\alpha}\right)^{2}.$ (66)

Therefore we arrive at the estimate (up to the logarithmic factor)

$\delta\langle\varpi^{2}\rangle\sim-\frac{\alpha\epsilon^{2}}{\nu^{5}k_{f}^{6}}\sim-\frac{\alpha}{\nu
k_{f}^{2}}\frac{\epsilon}{\nu^{3}k_{f}^{4}}\langle\varpi^{2}\rangle_{0},$ (67)

where $\langle\varpi^{2}\rangle_{0}\sim\epsilon/\nu$ is the zero order
contribution to the second moment. We conclude that besides the factor
$\beta_{\nu}$ (14), controlling the perturbation series in the limit $\nu
k_{f}^{2}\gg\alpha$, the correction to the second moment contains the
additional small factor $\alpha/(\nu k_{f}^{2})$.

Next, we consider the opposite limiting case $\alpha\gg\nu k_{f}^{2}$. Then
expression (63) can be rewritten as

$\displaystyle-\frac{\epsilon^{2}\nu}{3\alpha^{5}}\int\frac{dp\,dm}{(2\pi)^{2}}pm(m^{2}-p^{2})^{2}{\tilde{\varXi}(p)}{\tilde{\varXi}(m)}$
$\displaystyle\int\frac{d\psi}{2\pi}\frac{\sin^{2}\psi}{m^{2}+p^{2}+2pm\cos\psi}$
$\displaystyle\left[p^{2}m^{2}+2(p^{2}+m^{2})pm\cos\psi+4p^{2}m^{2}\cos^{2}\psi\right],$

where $\psi$ is the angle between $\bm{p}$ and $\bm{m}$, and we assumed $p\sim
m\sim k_{f}$. Taking the integral over $\psi$, we find

$\displaystyle\delta\langle\varpi^{2}\rangle=-\frac{\epsilon^{2}\nu}{12\alpha^{5}}\int\frac{dp\,dm}{(2\pi)^{2}}{\tilde{\varXi}(p)}{\tilde{\varXi}(m)}$
$\displaystyle pm(m^{2}-p^{2})^{2}(m^{2}+p^{2}-|m^{2}-p^{2}|).$ (68)

The integrand here is positive and has no singularities. Therefore we arrive
at the estimate

$\displaystyle\delta\langle\varpi^{2}\rangle\sim-\frac{\epsilon^{2}\nu
k_{f}^{6}}{\alpha^{5}}\sim-\frac{\epsilon\nu
k_{f}^{4}}{\alpha^{4}}\langle\varpi^{2}\rangle_{0},$ (69)

where we have taken into account the estimate
$\langle\varpi^{2}\rangle_{0}\sim\epsilon k_{f}^{2}/\alpha$. For the
particular function (12) we find from Eq. (68)

$\delta\langle\varpi^{2}\rangle=-\frac{4\epsilon^{2}\nu
k_{f}^{6}}{3\alpha^{5}}.$

The result confirms estimate (69). We conclude that besides the factor
$\beta_{\alpha}$ (14), controlling the perturbation series in the limit
$\alpha\gg\nu k_{f}^{2}$, the correction to the second moment contains the
additional small factor $\nu k_{f}^{2}/\alpha$.

Thus, we have established estimates for the correction to the second moment of
vorticity $\delta\langle\varpi^{2}\rangle$ in both limiting cases, see Eqs.
(67) and (69). They are anomalously small in accordance with the general
properties following from the energy balance and the enstrophy balance. It is
interesting that estimates (67) and (69) contain both kinetic coefficients,
$\nu$ and $\alpha$, whereas naively one may think that one of the kinetic
coefficients would drop.

$-t$$\bm{k}$$\bm{q}$$\bm{p}$$\bm{x}$$\bm{y}$$\bm{r}$ Figure 7: One of the
first order contributions to the triple correlation function of vorticity.

##  V Triple correlation function of vorticity

In zero approximation, when one neglects the nonlinear interaction, the
contribution to the triple correlation function of vorticity is equal to zero.
Therefore the correlation function directly determines properties of the
nonlinear interaction. Here we consider the first contribution to the triple
correlation function of vorticity. It is determined by the first term in the
expansion over the third-order term (17) in the effective action in the
functional integral determining the triple correlation function of vorticity.

The first contribution to the triple correlation function of vorticity is
determined by three similar diagrams, one of which is shown in Fig. 7. The
corresponding analytical expression for the simultaneous triple correlation
function is

$\displaystyle\langle\varpi(\bm{x})\varpi(\bm{y})\varpi(\bm{r})\rangle=\int_{0}^{\infty}dt\int\frac{d^{2}k\,d^{2}q\,d^{2}p}{(2\pi)^{4}}$
$\displaystyle\delta(\bm{k}+\bm{q}+\bm{p})\exp[i\bm{k}\bm{x}+i\bm{q}\bm{y}+i\bm{p}\bm{r}]$
$\displaystyle\left\\{\left(\frac{1}{q^{2}}-\frac{1}{k^{2}}\right)(q_{2}k_{1}-q_{1}k_{2})F(t,\bm{k})F(t,\bm{q})G(t,\bm{p})\right.$
$\displaystyle+\left(\frac{1}{k^{2}}-\frac{1}{p^{2}}\right)(k_{2}p_{1}-k_{1}p_{2})F(t,\bm{k})F(t,\bm{p})G(t,\bm{q})$
$\displaystyle\left.+\left(\frac{1}{p^{2}}-\frac{1}{q^{2}}\right)(p_{2}q_{1}-p_{1}q_{2})F(t,\bm{p})F(t,\bm{q})G(t,\bm{k})\right\\}.$
(70)

Substituting here expression (32) and taking the integral over time $t$, one
finds

$\displaystyle\langle\varpi(\bm{x})\varpi(\bm{y})\varpi(\bm{r})\rangle=\int\frac{d^{2}p\,d^{2}q\,d^{2}k}{(2\pi)^{4}}$
$\displaystyle\frac{\exp(i\bm{k}\bm{x}+i\bm{q}\bm{y}+i\bm{p}\bm{r})}{\nu(k^{2}+q^{2}+p^{2})+3\alpha}\delta(\bm{k}+\bm{q}+\bm{p})$
$\displaystyle(q_{2}k_{1}-q_{1}k_{2})\left\\{\left(\frac{1}{q^{2}}-\frac{1}{k^{2}}\right)F(\bm{q})F(\bm{k})\right.$
$\displaystyle+\left(\frac{1}{k^{2}}-\frac{1}{p^{2}}\right)F(\bm{k})F(\bm{p})$
$\displaystyle\left.+\left(\frac{1}{p^{2}}-\frac{1}{q^{2}}\right)F(\bm{p})F(\bm{q})\right\\}.$

The expression can be rewritten as

$\displaystyle\langle\varpi(\bm{x})\varpi(\bm{y})\varpi(\bm{r})\rangle=\varPhi(\bm{x},\bm{y},\bm{r})$
$\displaystyle+\varPhi(\bm{y},\bm{r},\bm{x})+\varPhi(\bm{r},\bm{x},\bm{y}),$
(71)

where

$\displaystyle\varPhi(\bm{x},\bm{y},\bm{0})=\int\frac{d^{2}q\,d^{2}k}{(2\pi)^{4}}\frac{\exp(i\bm{k}\bm{x}+i\bm{q}\bm{y})}{\nu(k^{2}+q^{2}+p^{2})+3\alpha}$
$\displaystyle
F(\bm{q})F(\bm{k})(q_{2}k_{1}-q_{1}k_{2})\left(\frac{1}{q^{2}}-\frac{1}{k^{2}}\right).$
(72)

To restore $\varPhi(\bm{x},\bm{y},\bm{r})$ from
$\varPhi(\bm{x},\bm{y},\bm{0})$ one should substitute
$\bm{x}\to\bm{x}-\bm{r}$, $\bm{y}\to\bm{y}-\bm{r}$ into
$\varPhi(\bm{x},\bm{y},\bm{0})$.

Let us analyze the limiting case $\alpha\ll\nu k_{f}^{2}$. Then for the
separations between points $\bm{x},\bm{y},\bm{r}$ of the order of $k_{f}^{-1}$
we find from Eqs. (71) and (72)

$\langle\varpi\varpi\varpi\rangle\sim\epsilon^{2}\nu^{-3}k_{f}^{-2},$ (73)

where we substituted $F\sim\epsilon/(\nu k_{f}^{2})$. Thus, we find the
estimate

$\frac{\langle\varpi\varpi\varpi\rangle^{2}}{(\langle\varpi^{2}\rangle_{0})^{3}}\sim\frac{\epsilon}{\nu^{3}k_{f}^{4}}.$
(74)

The quantity in the right-hand side of Eq. (74) is the parameter $\beta_{\nu}$
(14). We conclude that the perturbation parameter controls the third order
correlation function in the limit $\alpha\ll\nu k_{f}^{2}$.

In the opposite limiting case $\alpha\gg\nu k_{f}^{2}$, we find for the
separations $\sim k_{f}^{-1}$ from Eqs. (71) and (72)

$\langle\varpi\varpi\varpi\rangle\sim\epsilon^{2}k_{f}^{4}\alpha^{-3},$ (75)

where we substituted $F\sim\epsilon/\alpha$. Thus, we obtain

$\frac{\langle\varpi\varpi\varpi\rangle^{2}}{(\langle\varpi^{2}\rangle_{0})^{3}}\sim\frac{\epsilon
k_{f}^{2}}{\alpha^{3}}.$ (76)

The quantity in the right-hand side of Eq. (76) is the parameter
$\beta_{\alpha}$ (14). We conclude that the perturbation parameter controls
the third order correlation function in the limit $\alpha\gg\nu k_{f}^{2}$.

In the limit $\alpha\gg\nu k_{f}^{2}$, one can explicitly find
$\varPhi(\bm{x},\bm{y},\bm{0})$ for the particular function (12), since in
this case integral (72) is reduced to a Gaussian one. The result is

$\displaystyle\varPhi(\bm{x},\bm{y},\bm{0})=\frac{\epsilon^{2}k_{f}^{4}}{3\alpha^{3}}\left(\frac{\partial}{\partial
x_{1}}\frac{\partial}{\partial y_{2}}-\frac{\partial}{\partial
y_{1}}\frac{\partial}{\partial x_{2}}\right)$
$\displaystyle\left[(x^{2}-y^{2})\exp\left(-\frac{k_{f}^{2}x^{2}}{2}-\frac{k_{f}^{2}y^{2}}{2}\right)\right].$
(77)

These expression is in agreement with the general estimate (75).

We know that the third moment of vorticity is zero, see Section III. Therefore
the triple correlation function should be small for the separations between
the points $\bm{x},\bm{y},\bm{r}$ much smaller than $k_{f}^{-1}$. To analyze
integral (72) in this case, we introduce the parametrization

$\displaystyle k_{1}=k\cos\varphi,\quad k_{2}=k\sin\varphi,\quad
q_{1}=q\cos(\varphi+\psi),$ $\displaystyle q_{2}=q\sin(\varphi+\psi),\quad
k_{1}q_{2}-k_{2}q_{1}=kq\sin\psi.$

Taking the integral over $\varphi$, one obtains

$\displaystyle\varPhi(\bm{x},\bm{y},\bm{0})=\int\frac{dq\,dk\,d\psi}{(2\pi)^{3}}F(q)F(k)$
$\displaystyle\frac{(k^{2}-q^{2})J_{0}(a)\sin\psi}{2\nu(k^{2}+q^{2}+kq\cos\psi)+3\alpha},$
(78)

where

$\displaystyle a^{2}=k^{2}x^{2}+q^{2}y^{2}+2kq\bm{x}\bm{y}\cos\psi$
$\displaystyle+2kq(x_{1}y_{2}-x_{2}y_{1})\sin\psi.$ (79)

The quantity $a^{2}$ is positively defined.

Expanding $J_{0}$ in Eq. (78) in $a^{2}$ and performing the integration we
find a series over the separations. One can easily find that the terms of the
expansion of zero, second, fourth and sixth orders in the separations produce
zero contributions to
$\langle\varpi(\bm{x})\varpi(\bm{y})\varpi(\bm{r})\rangle$ in accordance with
Eq. (71). Thus, the first non-vanishing contribution to
$\langle\varpi(\bm{x})\varpi(\bm{y})\varpi(\bm{r})\rangle$ is of the eighth
order in the separations. The property can be directly checked for the
explicit expression (77).

If $\alpha\ll\nu k_{f}^{2}$ then there is an interval of the separations
between points
$k_{f}^{-1}\ll|\bm{x}-\bm{y}|,|\bm{x}-\bm{r}|,|\bm{r}-\bm{y}|\ll\sqrt{\nu/\alpha}$
where the triple correlation function of vorticity decays according to a
universal power law. To find the law, one should note that for such distances
the integral (72) is gained at the wave vectors $k,q$ of the order of the
inverse separations, $\sqrt{\alpha/\nu}\ll k,q\ll k_{f}$. Therefore, one can
put $\alpha\to 0$ and expand the product $F(\bm{k})F(\bm{q})$ in $k^{2}$,
$q^{2}$, since zero-order term of the expansion (that is a constant) leads to
the function $\varPhi(\bm{x},\bm{y},\bm{0})$ producing zero contribution to
the triple correlation function (71). Using equation (64), we obtain

$\displaystyle\langle\varpi(\bm{x})\varpi(\bm{y})\varpi(\bm{r})\rangle=-C\frac{\epsilon^{2}}{(\nu
k_{f}^{2})^{3}}$
$\displaystyle\left[(x_{1}-r_{1})(y_{2}-r_{2})-(x_{2}-r_{2})(y_{1}-r_{1})\right]$
$\displaystyle\left\\{\frac{4}{R^{6}}\left[\frac{(\bm{x}-\bm{r})(2\bm{y}-\bm{x}-\bm{r})}{(\bm{x}-\bm{r})^{2}}-\right.\right.$
$\displaystyle\left.\frac{(\bm{y}-\bm{r})(2\bm{x}-\bm{y}-\bm{r})}{(\bm{y}-\bm{r})^{2}}-\frac{(\bm{x}-\bm{y})(2\bm{r}-\bm{y}-\bm{x})}{(\bm{y}-\bm{x})^{2}}\right]$
$\displaystyle+\frac{1}{R^{4}}\left[\frac{(\bm{x}-\bm{r})(2\bm{y}-\bm{x}-\bm{r})}{(\bm{x}-\bm{r})^{4}}-\right.$
(80)
$\displaystyle\left.\left.\frac{(\bm{y}-\bm{r})(2\bm{x}-\bm{y}-\bm{r})}{(\bm{y}-\bm{r})^{4}}-\frac{(\bm{x}-\bm{y})(2\bm{r}-\bm{y}-\bm{x})}{(\bm{y}-\bm{x})^{4}}\right]\right\\},$

where

$R^{2}=(\bm{x}-\bm{r})^{2}+(\bm{y}-\bm{r})^{2}+(\bm{x}-\bm{y})^{2}.$

The constant $C$ is equal to unity for the forcing with spectrum (12).

Expression (80) is invariant under permutations of $\bm{x},\bm{y},\bm{r}$.
Note that it is equal to zero when the points $\bm{x},\bm{y},\bm{r}$ are
vertices of isosceles triangle. It is a consequence of the general symmetry:
any odd correlation function of $\varpi$ should change its sign when reflected
relative to an arbitrary axis. One can find a simple expression by putting
$\bm{r}\to 0$ and assuming $x\ll y$. In this case, expression (80) takes the
form:

$\displaystyle\langle\varpi(\bm{x})\varpi(\bm{y})\varpi(\bm{0})\rangle=\frac{\epsilon^{2}}{4(\nu
k_{f}^{2})^{3}}\frac{\sin(2\varphi)}{x^{2}y^{2}},$ (81)

where $\varphi$ is the angle between the vectors $\bm{x}$ and $\bm{y}$.

When separations are much larger than the length $\sqrt{\nu/\alpha}$, the
triple correlation function
$\langle\varpi(\bm{x})\varpi(\bm{y})\varpi(\bm{r})\rangle$ decays
exponentially. The argument of the exponent is determined by the same length
$\sqrt{\nu/\alpha}$. Note that this behavior is universal, being independent
of the pumping details.

##  VI Direct Numerical Simulations

In this section, we present the results of direct numerical simulations (DNS),
which were performed to validate analytical expressions that underlie our
analysis. The DNS results were obtained by integrating the Navier-Stokes
equation in the vorticity formulation (3) using the GeophysicalFlows.jl
pseudospectral code Constantinou21 , fully dealiazed by the two-thirds rule.
This code can be executed on the GPU, resulting in high computational
performance. In our simulations, we used a periodic square domain of size
$L=2\pi$.

Figure 8: Correction to the one-dimensional enstrophy spectrum. The DNS
results (blue line) are consistent with the numerical integration of
expressions (42) and (43) (markers).

We performed several runs to illustrate various features of obtained results.
In the first run, our goal was to verify expressions (42) and (43), which
describe corrections $F_{a}$ and $F_{b}$ to the pair correlation function of
vorticity. We set the parameters $\nu=1.6\times 10^{-3}$ and $\alpha=1$, so
that $\alpha=\nu k_{f}^{2}$. The random forcing was isotropic in space,
shortly correlated in time and had a covariance spectrum (12) with $k_{f}=25$.
The grid resolution was $512^{2}$, which made it possible to collect the
necessary statistics. Since $F_{a},\,F_{b}\propto\epsilon^{2}$, we attempted
to find the largest possible value of the parameter $\epsilon$ at which
perturbation theory could still be applied. This will make it easier to
distinguish the small correction $F_{a}+F_{b}$ from the background noise of
the main contribution $F$, see expression (31). We found that
$\epsilon=1.6\times 10^{-2}$ corresponding to the values of
$\beta_{\nu}=\beta_{\alpha}=10$ leads to satisfactory results.

The initial state was a state of rest, and we first conducted simulations for
some time to reach a statistical steady-state, observed by the saturation of
the total kinetic energy. The integration time step was $\Delta t=0.005$,
which was also the correlation time of the exciting force. Once stationary, we
saved the vorticity field and ran several simulations in parallel with
different realizations of the random force, using the saved data as the
initial condition. Then we stored data every $200$ integration steps and
gathered a total of about $0.6\times 10^{6}$ snapshots of the vorticity field.

Next, to compare DNS with analytical findings, we calculated the time-averaged
one-dimensional enstrophy spectrum $k(F+\delta F)/(2\pi)$ and subtracted from
it the contribution corresponding to the linear approximation (31). The
results are shown in Fig. 8 and they agree with the answer obtained by
numerical integration of expressions (42) and (43), which confirms their
correctness. The slight noise noticeable in the DNS results is associated with
a moderate volume of statistics.

Figure 9: Enstrophy spectrum under narrowband pumping (82) with $k_{f}=50$ and
$\delta_{f}=1.5$. The solid line shows DNS results, the dashed line
corresponds to the linear approximation (31), and the markers are obtained by
numerical integration of expression (42).

In the second run, our purpose was to considered the random forcing, which has
a narrow spectrum in $k$-space. Pumping with such a spatial structure is
commonly utilized in the numerical modeling of turbulent problems Chertkov2007
; Laurie2014 ; Frishman2018 ; Parfenyev2022  and, as mentioned in Section
IV.1, might result in some peculiarities. To illustrate the scenario, we
performed DNS with

$\tilde{\varXi}(\bm{k})=\frac{\sqrt{2\pi}}{\delta_{f}k_{f}}\exp\left[-(k-k_{f})^{2}/2\delta_{f}^{2}\right],$
(82)

where we set $k_{f}=50$ and $\delta_{f}=1.5$. The other parameters were
$\alpha=0.01$, $\nu=1.6\times 10^{-5}$, $\epsilon=2.56\times 10^{-11}$ and the
grid resolution was $1024^{2}$. The integration time step was $\Delta t=0.01$
and after reaching a steady-state, we saved the vorticity field every $100$
integration steps, and collected $10^{4}$ snapshots.

The results for the time-averaged one-dimensional enstrophy spectrum
$k(F+\delta F)/(2\pi)$ are presented in Fig. 9. The linear approximation (31)
well describes the main peak near $k=k_{f}$, and the correction $F_{a}$
accounts for the wide plateau on which it is located. Note that in this case
DNS does not require the collection of large statistics, since correction
$F_{a}$ is separated from the main contribution in the Fourier space. The
contribution $F_{b}$ is localized near $k_{f}$ and is not visible against the
background of the main peak. These results confirm the qualitative analysis
carried out in Subsection IV.1, and considering narrowband pumping allows us
to double-check the correctness of the expression (42) with a few statistics
in hand.

Figure 10: Triple correlation function of vorticity: comparison of analytical
expression (77) with DNS.

In the third run, we focused on triple correlation function of vorticity (71).
Analytical calculations yielded a concise answer (77) in the limiting case
$\nu k_{f}^{2}/\alpha\ll 1$ and when the pumping covariance spectrum had the
form (12). We reproduced these conditions in DNS with parameters
$\alpha=1000$, $\nu=1.6\times 10^{-3}$, $k_{f}=25$, $\epsilon=160$. The grid
resolution was $512^{2}$ and the integration time step was set to $\Delta
t=10^{-5}$. We started simulations from the state of rest, and once
stationary, we saved the vorticity field and ran six simulations in parallel
with different random realizations of forcing. We gathered a total of about
$0.6\times 10^{6}$ snapshots of the vorticity field, storing the data every
$100$ integration steps.

Next, we calculated a triple correlation function
$\langle\omega_{1}\omega_{2}\omega_{3}\rangle$, where the values of vorticity
were taken in the vertices of right triangle with legs $a$ and $b$. The
results are presented in Fig. 10. The solid lines show the analytical
calculations using expression (77), and the markers represent the DNS results
averaged over all collected data. The error bars indicate the maximum and
minimum values obtained from six parallel simulations. Note that to increase
statistics we used the symmetries of the problem: rotations of the triangle do
not change the correlation function, and reflections relative to straight
lines passing through the legs change its sign. Overall, it can be concluded
that the analytical and numerical results are in good agreement with each
other, which confirms the correctness of the developed perturbation method.

##  VII Conclusion

We examined properties of the perturbation theory for two-dimensional fluid
flow in the framework of the model where the pumping force, exciting the flow,
is random in space and shortly correlated in time. The model admits advanced
analytical calculations revealing qualitative properties of the perturbation
theory. The peculiarity of two-dimensional turbulence is existence of two
dissipation mechanisms: bottom friction and viscosity. There are two
dimensionless parameters (14) controlling the perturbation series for weak
bottom friction and weak viscosity, correspondingly. We performed the
calculations of the first correction to the pair correlation function of
vorticity and the main contribution to the triple correlation function of
vorticity, which enable one to evaluate the applicability region of the
perturbation series.

There are some interesting features revealed by the calculations. Fourier
transform of the correction to the pair correlation function changes its sign
twice. It is the consequence of the integral relations (46) and (45),
reflecting the energy balance and the enstrophy balance. The third moment of
vorticity is zero in our model thanks to the symmetry reasoning, see Section
III. At small separations between the points the triple correlation function
of vorticity behaves as the eighth degree of the separations. If the bottom
friction is weak then there is the universal region of distances where the
triple correlation function scales as the power $-4$ of the distances.

Surprisingly, the corrections to the second moment of vorticity do not obey
the general properties of the perturbation series. In both limit cases, $\nu
k_{f}^{2}\gg\alpha$ and $\nu k_{f}^{2}\ll\alpha$, the second moment is much
less than the estimate, related to the corresponding small perturbation
parameters (14), see Subsection IV.2. Such behavior is explained by the energy
and the enstrophy balances, leading to disappearing the main contributions to
the correction to the second moment, see Subsection IV.2. Thus one should be
careful with evaluating the parameter controlling the perturbation series via
the second moment of vorticity.

To confirm our analytical results, we performed direct numerical simulations
of two-dimensional turbulence excited by an external force shortly correlated
in time. The simulations were performed in a square box with periodic boundary
conditions. The results of the simulations are in excellent agreement with our
analytical predictions both for corrections to the pair correlation function
and for the triple correlation function. We tested also the pumping with a
narrow spectrum in $\bm{k}$-space. The results are in agreement with the
analytical calculations as well.

One expects, that turbulence in a $2d$ film is excited where both parameters
(14) are large. The character of the state depends, generally, on the ratio
$\alpha/(\nu k_{f}^{2})$. However, large-scale fluctuations (with scales
larger than the pumping length) are more sensitive to bottom friction whereas
small-scale fluctuations (with scales smaller than the pumping length) are
more sensitive to viscosity. Thus, the description of the turbulent state can
be complicated, being dependent on the ratio $\nu k_{f}^{2}/\alpha$. The
problem needs a special investigation.

The passive regime of the flow fluctuations where the perturbation series is
applicable can be realized inside the coherent vortices appearing via the
inverse cascade in two-dimensional turbulence KL15 ; kolokolov2016structure ;
kolokolov2016velocity ; frishman17 . The reason is that the fluctuations are
suppressed by the coherent vortex flow, which behaves locally as a shear flow.
To establish the criterion of the applicability of the approach one should
generalize the scheme developed in the present work to the case of the two-
dimensional flow excited in the strong external shear flow. First steps in
this direction were made in the papers KLT23a ; KLT23b .

###### Acknowledgements.

I.V.K. acknowledges support from the Ministry of Science and Higher Education
of the Russian Federation. The work of V.V.L. and V.M.P. was performed in the
Laboratory “Modern Hydrodynamics” created in frames of Grant 075-15-2022-1099
of the Ministry of Science and Higher Education of the Russian Federation in
Landau Institute for Theoretical Physics of RAS and was supported by Grant
23-72-30006 of Russian Science Foundation. V.M.P. also acknowledges support
from the Foundation for the Advancement of Theoretical Physics and Mathematics
“BASIS”, Project No. 22-1-3-24-1. Simulations were performed on the cluster of
the Landau Institute for Theoretical Physics of RAS.

## References

  * (1) J. Sommeria, J . Fluid Mech. 170, 139-168 (1986).
  * (2) P. Tabeling, Phys. Rep. 362, 1-62 (2002).
  * (3) H.Xia, M.Shats, and G.Falkovich, Phys. Fluids 21, 125101 (2009).
  * (4) A. V. Orlov, M. Yu. Brazhnikov, A.A.Levchenko, Pis’ma v ZhETF 107, 166 (2018) [JETP Letters 107, 157 (2018)].
  * (5) D. Molenaar, H.J.H. Clercx, G.J.F. van Heijst, Physica D 196, 329 (2004).
  * (6) P. K. Mishra, J. Herault, S. Fauve, and M. K. Verma, Phys. Rev. E 91, 053005 (2015).
  * (7) A. N. Doludenko, S. V. Fortova, I. V. Kolokolov, and V. V. Lebedev, Annals of Physics, 447, 169072 (2022).
  * (8) V.M. Parfenyev, S.S. Vergeles, Physics of Fluids 33, 115128 (2021).
  * (9) H. Kellay, W.I. Goldburg, Reports on Progress in Physics 65, 845 (2002).
  * (10) S.V. Yablonskii, N.M. Kurbatov, and V.M. Parfenyev, Phys. Rev. E 95, 012707 (2017).
  * (11) I.V. Kolokolov and V.V. Lebedev, Pis’ma v ZhETF, 101, 181 (2015) [JETP Letters, 101, 164 (2015)].
  * (12) I.V. Kolokolov and V.V. Lebedev, Phys. Rev. E, 93, 033104 (2016).
  * (13) I. V. Kolokolov and V. V. Lebedev, J. of Fluid Mechanics, 809, R2 (2016).
  * (14) A.Frishman, J.Laurie, G.Falkovich, Phys. Rev. Fluids 2, 032602 (2017).
  * (15) H. W. Wyld, Ann. Phys. (N.Y.) 14, 134 (1961).
  * (16) P. C. Martin, E. D. Siggia, and H. A. Rose, Phys. Rev. A 8, 423 (1973).
  * (17) S. Weinberg, The Quantum Theory of Fields, Vol. I (Foundations) Cambridge University Press, Cambridge, U. K., (1995).
  * (18) J. Hertz, Y. Roudi, P. Sollich, J. Phys. A: Math. Theor. 50, 033001 (2017).
  * (19) G. C. Wick, Phys. Rev. 80, 268 (1950).
  * (20) I.V.Kolokolov, V.V.Lebedev, M.M.Tumakova, Pis’ma v ZhETF 117, 127 (2023); [JETP Letters, 117, 122 (2023).]
  * (21) I.V.Kolokolov, V.V.Lebedev, M.M.Tumakova, ZhETF 163, 881 (2023); [JETP 136, 785 (2023)].
  * (22) N. C. Constantinou, G. L. Wagner, L. Siegelman, B. C. Pearson, and A. Paloczy, Journal of Open Source Software 6, 3053 (2021).
  * (23) M. Chertkov, C. Connaughton, I. Kolokolov, and V. Lebedev, Phys. Rev. Lett. 99, 084501 (2007).
  * (24) J. Laurie, G. Boffetta, G. Falkovich, I. Kolokolov, and V. Lebedev, Phys. Rev. Lett. 113, 254503 (2014).
  * (25) A. Frishman and C. Herbert, Phys. Rev. Lett. 120, 204505 (2018).
  * (26) V. Parfenyev, Phys. Rev. E 106, 025102 (2022).