# Velocity profiles of cyclones and anticyclones in a rotating turbulent flow

Vladimir M. Parfenyev<EMAIL_ADDRESS>Ivan A. Vointsev  Alyona O. Skoba
Sergey S. Vergeles<EMAIL_ADDRESS>Landau Institute for Theoretical Physics,
Russian Academy of Sciences,
1-A Akademika Semenova av., 142432 Chernogolovka, Russia  National Research
University Higher School of Economics, Faculty of Physics,
Myasnitskaya 20, 101000 Moscow, Russia

###### Abstract

Strong rotation makes an underlying turbulent flow quasi-two-dimensional that
leads to the upscale energy transfer. Recent numerical simulations show that
under certain conditions, the energy is accumulated at the largest scales of
the system, forming coherent vortex structures known as condensates. We
analytically describe the interaction of a strong condensate with weak small-
scale turbulent pulsations and obtain an equation that allows us to determine
the radial velocity profile $U(r)$ of a coherent vortex. When external
rotation is fast, the velocity profiles of cyclones and anticyclones are
identical to each other and are well described by the dependence
$U(r)\propto\pm r\ln(R/r)$, where $R$ is the transverse size of the vortex. As
the external rotation decreases, this symmetry disappears: the maximum
velocity in cyclones is greater and the position of the maximum is closer to
the axis of the vortex in comparison with anticyclones. Besides, our analysis
shows that the size $R$ of the anticyclone cannot exceed a certain critical
value, which depends on the Rossby and Reynolds numbers. The maximum size of
the cyclones is limited only by the system size under the same conditions. Our
predictions are based on the linear evolution of turbulent pulsations on the
background of the coherent vortex flow and are accompanied by estimates
following from the nonlinear Navier-Stokes equation.

##  I Introduction

It is well known that rotation plays a crucial role in determining the
properties of underlying turbulence in three-dimensional flows, see, e.g., the
recent review [1]. When rotation is strong, the velocity field tends to become
quasi-two-dimensional varying weakly along the rotation axis following the
Taylor-Proudman theorem [2, 3]. For this reason, it is often said that strong
rotation effectively makes the flow two-dimensional and opens the way for an
inverse energy cascade (from small to large scales), which is characteristic
of two-dimensional turbulence [4]. If the damping of the large-scale motion in
the system is small enough, then the inverse cascade can lead to the formation
of large-scale columnar vortices that are stable over many periods of their
rotation (the so-called condensate).

The spontaneous formation of long-living columnar vortices aligned in the
direction of rotation in three-dimensional systems has been observed both in
experiments [5, 6, 7, 8, 9, 10, 11] and in numerical simulations [12, 13, 14,
15, 16, 17]. Remarkably, the distribution of these vortices is asymmetric with
the predominance of cyclones (co-rotating with the external rotation) over
anticyclones. The qualitative explanation is that anticyclones weaken the
external rotation due to their own rotation which prevents the development of
the inverse energy cascade inside them. Quantitatively, the cyclone-
anticyclone asymmetry is studied by statistical methods using various
correlation functions [10, 18, 19] that involve spatial averaging over the
entire observation area.

This approach is justified, since in most of the works mentioned above, the
large-scale flow is not strictly stationary and at large times required for
collecting statistics, the configuration of cyclones and anticyclones changes.
The only exception among the experiments known to us was reported in Ref. 7,
where truly stationary vortices were observed and their spatial arrangement
was fixed by inhomogeneity in the permanent forcing of the flow. However, the
results of recent numerical simulations demonstrate that with increasing
external rotation, columnar vortices become very long-lived [17]. Besides,
their transverse size increases; it exceeds the forcing scale and becomes
comparable to the system size $L$. The velocity amplitude $U\sim
L\sqrt{\epsilon/\nu}$ inside the vortex is the result of a balance between the
forcing power $\epsilon$ per unit mass and the viscous friction inside the
condensate, where $\nu$ is the kinematic viscosity of a fluid. In this state,
new questions arise related to the structure of an individual vortex. How does
the azimuthal velocity of the vortex change depending on the distance to its
axis? How does the vortex velocity profile depend on the global parameters of
the problem? Is there any difference between cyclones and anticyclones?

The first attempt to address these questions was made in Ref. 20, where the
authors studied the limiting case of fast rotation (Rossby number
$\mathrm{Ro}_{\scriptscriptstyle R}=\sqrt{\epsilon/\nu}/(2\Omega_{0})\ll 1$,
where $\Omega_{0}$ is the angular velocity of the external rotation) and
small-scale forcing ($k_{f}L\gg 1/\mathrm{Ro}_{\scriptscriptstyle R}$, where
$k_{f}$ is the characteristic wavenumber of the forcing). In this limiting
case, the power balance between the forcing and the viscous friction is local
in space. The forcing first excites small-scale inertial waves and then their
kinetic energy is transferred to the large-scale vortex flow. This leads to
the universal velocity profile $U(r)\sim\pm\sqrt{\epsilon/\nu}\,r\ln(R/r)$,
which describes both cyclones and anticyclones. Here $r$ is the distance to
the vortex axis and $R\lesssim L$ is the transverse vortex size. The absence
of cyclone-anticyclone asymmetry in the limit of small Rossby numbers is in
agreement with the results of numerical simulations [17]. Note that the local
power balance between the forcing and the viscous dissipation inside the
coherent flow can be achieved in two-dimensional flow as well, resulting in
the same linear-logarithmic velocity profile observed in Ref. 21. As for
experimental attempts to measure the velocity profiles of coherent vortices,
we would like to mention works 5, 6, although the condensates in them were
clearly underdeveloped (large-scale vortices live for a relatively short time
and their size is small compared to the system size).

Here we consistently study the structure of the vortex condensates assuming
that the Rossby number is moderate and analyze the cyclone-anticyclone
asymmetry in the developed condensate. For this purpose, we extend the
approach proposed in Ref. 20 to the case of moderately low Rossby numbers. We
consider the turbulent three-dimensional pulsations on the background of the
coherent vortex flow in the linear regime, assuming that the rapid distortion
theory is applicable. The turbulent pulsations produce the Reynolds stress
which maintains the coherent vortex. Then we obtain the mean velocity profile
$U(r)$ and analyze the difference between cyclones and anticyclones. Using the
estimates following from the full nonlinear Navier-Stokes equation, we derive
the criteria when the coherent vortices can exist.

Actually, the cyclone-anticyclone asymmetry can be observed even at a small
Rossby number if the system size is large enough. Inside a coherent vortex,
the external rotation is superimposed by self-rotation
$U(r)/r\sim\pm\sqrt{\epsilon/\nu}\ln(R/r)$ induced by the vortex. The total
local rotation rate is equal to $\Omega=\Omega_{0}+U(r)/r$, so the local
Rossby number $\mathrm{Ro}(r)=r\partial_{r}(U/r)/(2\Omega)$ formally grows
indefinitely with approaching the axis of anticyclone. However, the linear-
logarithmic profile is not applicable below some core radius $r_{c}$, where
the local power balance approximation is no longer valid and the profile
changes to a rigid body rotation [20]. Thus, one can estimate the maximum
possible radius of an anticyclone from the condition $2\ln(R_{max}/r_{c})\sim
1/\mathrm{Ro}_{\scriptscriptstyle R}$, which ensures that the local Rossby
number inside the anticyclone core remains below unity to prevent the
development of direct energy cascade (characteristic of a three-dimensional
turbulent flow with weak rotation). The restriction $R<R_{max}$ for
anticyclones means that the cyclone-anticyclone asymmetry occurs in the case
$L\gg R_{\max}$, since in cyclones the self-rotation increases the total
rotation rate, so there is no similar limitation for them.

##  II Momentum balance

We consider a three-dimensional incompressible turbulent fluid that rotates
around the $Z$-axis with a constant angular velocity $\bm{\Omega}_{0}$. In the
rotating frame, the system is described by the Navier-Stokes equation with an
additional term describing the Coriolis force and with an external force
$\bm{f}$ random in time and space with zero mean (pumping or forcing), which
continuously excites small-scale turbulent pulsations

$\partial_{t}\bm{v}+(\bm{v}\nabla)\bm{v}+2\left[\bm{\Omega}_{0},\bm{v}\right]=-\nabla
p+\nu\Delta\bm{v}+\bm{f}.$ (1)

Equation (1) has to be supplemented by the incompressibility condition
$\mathrm{div}\,\bm{v}=0$. Here $\bm{v}$ is the fluid velocity, $p$ is the
effective pressure divided by the fluid mass density, which differs from the
physical pressure by addition of the potentials produced by the centrifugal
and gravitational forces, and $\nu$ is the kinematic viscosity coefficient. We
assume that the force $\bm{f}$ produces the energy flux per unit mass
$\epsilon$, has a correlation length in space equal to $1/k_{f}\ll L$, where
$L$ is the system size, and is shortly correlated in time. The Reynolds number
characterizing the pumping should be large,
${\mathrm{Re}}_{f}=\epsilon^{1/3}/\nu k_{f}^{4/3}\gg 1$, which is a
prerequisite for the overall flow to be turbulent.

Figure 1: Schematic of a coherent columnar vortex and the introduced
cylindrical $\\{{r},\varphi,z\\}$ and Cartesian $\\{\xi,\eta,z\\}$ coordinate
systems.

Next, we assume that a coherent columnar large-scale vortex has already formed
in our system. We introduce a cylindrical coordinate system
$\\{{r},\varphi,z\\}$ so that the vortex axis coincides with the direction of
the $Z$-axis, see Fig. 1. In the vortex flow region, the velocity field can be
represented as the sum of the mean flow $\bm{U}$ and turbulent pulsations
$\bm{u}$:

$\bm{v}=\bm{u}+\bm{U},\quad\bm{u}=\\{u_{{r}},u_{\varphi},u_{z}\\},\quad\bm{U}=\\{0,U,0\\},$
(2)

where, by definition, the average of the fluctuations is equal to zero,
$\langle\bm{u}\rangle=0$. Substituting this representation into equation (1),
averaging the result over turbulent pulsations, and projecting the equation
onto ${\bf e}_{\varphi}$, we obtain

$\partial_{t}U=-\left(\partial_{{r}}+\frac{2}{{r}}\right)\Pi^{\varphi{r}},\quad\Pi^{\varphi{r}}=\langle
u^{\varphi}u^{{r}}\rangle-\nu\Sigma,$ (3)

where $\Pi^{ij}$ is the Reynolds stress tensor collected with the viscous
term. Here we also introduced the local large-scale shear rate
$\Sigma={r}\,\partial_{{r}}(U/{r})$ which characterizes the difference between
local rotation in the vortex and rigid body rotation.

Figure 2: Reynolds stress tensor, for definiteness $\Sigma(r)<0$ is chosen.
The red arrow depicts the turbulent part of the $\xi$-component of the
momentum flux in $\eta$-direction. Virtual boundary is absent for the vortex
flow, it is added for comparison with the case of turbulent boundary layer
flow in a channel.

Note that the differential relationship between the shear rate profile
$\Sigma({r})$ and the velocity profile $U({r})$ can be rewritten in an
integral form:

$U({r})=-{r}\int\limits_{r}^{R}\frac{\Sigma({r}^{\prime})}{{r}^{\prime}}\mathrm{d}{r}^{\prime},$
(4)

where the distance $R$ plays a role of an integration constant. If
$\Sigma({r})$ has permanent sign, then $U({r})$ has the opposite sign inside
the tube of radius $R$, at ${r}<R$, and changes its sign outside at ${r}>R$.
It is reasonable to assume that beyond the vortex the mean flow is absent or
at least uncorrelated with the flow inside the vortex. Then the region
${r}\gtrsim R$ should be treated as the boundary of the vortex, where the
axial symmetry is violated.

In what follows, we consider a statistically stationary situation, therefore
$\partial_{t}U=0$ and the moment of forces $-{r}^{2}\Pi^{\varphi{r}}({r})$
acting on fluid in a circle of radius ${r}$ per unit length is also equal to
zero. As a result, we arrive at the equation $\Pi^{\varphi{r}}=0$ or

$\langle u^{\varphi}u^{{r}}\rangle-\nu\Sigma=0,$ (5)

which will further allow us to determine the velocity profile of a coherent
vortex. It is important to note that equation (5) requires the inequality

$\Sigma\langle u^{\varphi}u^{{r}}\rangle>0$ (6)

to be satisfied throughout the vortex flow region.

To discuss the meaning of expression (6), let us introduce the local Cartesian
coordinate system $\\{\xi,\eta,z\\}$, whose origin moves around the axis of
the vortex along a circle of radius $r$ with an angular velocity $U({r})/{r}$,
see Fig. 1. Its axes ${\bf e}_{\xi}$, ${\bf e}_{\eta}$ rotate with the same
angular velocity. In this coordinate system, the directions of unit vectors
${\bf e}_{r}$, ${\bf e}_{\varphi}$ remain constant, and the distance to the
vortex axis is given by $|{r}-\eta|$. The local assignment of the coordinate
system means that we are only interested in distances $\eta\ll r$.

Consider, for definiteness, a vortex with $\Sigma(r)<0$. Then the distribution
of the mean flow velocity near the origin of the introduced Cartesian
coordinate system will be as shown in Fig. 2. Inequality (6) means that the
turbulent part of the $\xi$-component of the momentum flux $\langle
u^{\varphi}u^{{r}}\rangle$ must be directed to the vortex axis. This
essentially distinguishes the situation in a vortex from, for example, a
turbulent flow along a stationary boundary. Within the geometry of Fig. 2, the
plane $\eta=0$ can correspond to the boundary, then the fluid flow will be
directed along the $\xi$-axis. In this case, the turbulent momentum flux must
be directed towards the boundary, see, e.g., Ref. 22, that is, in the opposite
direction. In the case of a vortex, inequality (6) is fulfilled due to the
local rotation of the fluid $\Omega$, which is the sum of the external
rotation $\Omega_{0}$ and the local addition $U({r})/{r}$,

$\Omega({r})=\Omega_{0}+U({r})/{r},$ (7)

see also Ref. 23. If $U({r})$ has the same sign as $\Omega_{0}$, then the
vortex is called a cyclone. According to relation (4), cyclone corresponds to
different signs of $\Sigma$ and $\Omega$. An anticyclone is a vortex in which
$U$ and $\Omega_{0}$ have opposite signs. The vortex is statistically stable
if $\Omega(r)$ does not change its sign inside the vortex. Thus, the signs of
$\Omega$ and $\Sigma$ coincide in an anticyclone.

##  III Energy balance

Figure 3: Numerical data for dimensionless tangent Reynolds stress
$F(\mathrm{Ro},\gamma)=\Sigma\langle u^{\varphi}u^{{r}}\rangle/\epsilon$ at
(a) small and (b) large positive Rossby numbers
$\mathrm{Ro}=\Sigma/(2\Omega)>0$, which correspond to the anticyclone. In Fig.
3a the dashed line shows that the curve $q=1$ coincides with the curve $F=0$,
see equation (12). In Fig. 3b the function $F(\mathrm{Ro},\gamma)$ ceases to
depend on the parameter $\mathrm{Ro}$ above the dashed line following Ref. 24.

The sign of the Reynolds stress tensor component $\langle
u^{\varphi}u^{{r}}\rangle$ also determines the direction of the energy flow
between the large-scale mean current ${\bm{U}}$ and small-scale turbulent
pulsations. Let us analyze the energy balance based on the Navier-Stokes
equation (1). To do this, we multiply this equation by $\bm{v}$ and carry out
averaging over time:

$\partial_{t}\langle{T}^{kin}\rangle+\mathop{\mathrm{div}}\langle{T}^{kin}\bm{v}\rangle=-\partial_{i}\langle
v^{i}\,p\rangle+\nu\,\langle v^{i}\Delta v^{i}\rangle+\epsilon,$ (8)

where we have introduced the kinetic energy density ${T}^{kin}=\bm{v}^{2}/2$.

First, we consider the viscous term

$\displaystyle\langle v^{i}\Delta v^{i}\rangle=U^{i}\Delta U^{i}+\langle
u^{i}\Delta u^{i}\rangle=$ (9) $\displaystyle
U(\partial_{{r}}+2/{r})\Sigma+\frac{\Delta\langle{\bm{u}}^{2}\rangle}{2}-q\epsilon/\nu,$

where $q=\nu\langle\partial_{i}u^{k}\partial_{i}u^{k}\rangle/\epsilon$ is the
dimensionless rate of viscous dissipation of the kinetic energy of turbulent
pulsations. We assume that the local approximation is valid, which means that
the vortex size $R$ is much larger than the pumping scale $1/k_{f}$ and,
therefore, the term $\Delta\langle{\bm{u}}^{2}\rangle/2$ in equation (9) is
small compared to $q\epsilon/\nu$ as $1/(k_{f}R)^{2}$ and can be dropped.

Next, since we are analyzing the stationary regime, which implies the axial
symmetry of the vortex, the kinetic energy flux is directed only along the
radial direction,
$\langle{T}^{kin}{\bm{v}}\rangle=\langle{T}^{kin}u^{r}\rangle{\bf e}_{r}$. The
average amplitude is equal to

$\langle{T}^{kin}u^{r}\rangle=U\langle
u^{r}u^{\varphi}\rangle+\frac{1}{2}\langle{\bm{u}}^{2}u^{r}\rangle.$ (10)

As turbulent pulsations are assumed to be weak, the last term in expression
(10) can also be neglected.

Finally, we analyze the term with pressure. Due to the axial symmetry of the
vortex and its homogeneity along its axis, only the component $\langle
v^{r}p\rangle=\langle u^{r}p\rangle$ should be considered. Next, we set
$\langle u^{r}p\rangle=0$ following Ref. 25, where the absence of the
pressure-velocity correlation was adopted as a hypothesis based on empirical
observations [26]. As a result, we arrive to energy balance equation for the
statistically stable vortex flow, see below, which will be verified
numerically in an independent way in Section IV and thereby substantiate the
validity of this assumption in our case.

Collecting all the results together and equating the rate of change in the
average kinetic energy density to zero, we transform equation (8) to the form

$\frac{1}{{r}}\,\partial_{{r}}\big{(}{r}U\langle
u^{{r}}u^{\varphi}\rangle\big{)}=\nu
U(\partial_{{r}}+2/{r})\Sigma+\epsilon(1-q).$ (11)

Now we use $\langle u^{{r}}u^{\varphi}\rangle=\nu\Sigma$ which follows from
relation (5), and obtain the final expression for the energy balance in the
system

$\epsilon(1-q)=\Sigma\langle u^{{r}}u^{\varphi}\rangle\qquad\mathrm{or}\qquad
q=1-\frac{\Sigma}{\epsilon}\langle u^{r}u^{\varphi}\rangle.$ (12)

Thus, $\Sigma\langle u^{{r}}u^{\varphi}\rangle=\nu\Sigma^{2}$ is the power
density that is transferred from turbulent pulsations to the large-scale
coherent vortex. This power compensates for the viscous damping of the large-
scale flow, and it should be positive throughout the vortex flow region, if
only the coherent vortex exists, see equation (6). Let us stress again that
relation (12) takes place only if the vortex is large-scale, i.e. $k_{f}R\gg
1$.

Based on equation (12), one can estimate $\langle
u^{r}u^{\varphi}\rangle\sim\epsilon/\Sigma$ and then, according to expression
(5), one finds $\Sigma\sim\sqrt{\epsilon/\nu}$. For further consideration, it
is convenient to use these estimates to introduce the Rossby number for the
large-scale flow $\mathrm{Ro}_{\scriptscriptstyle
R}=\sqrt{\epsilon/\nu}/(2\Omega_{0})$. Together with the Reynolds number
$\mathrm{Re}_{f}=\epsilon^{1/3}/\nu k_{f}^{4/3}$ characterizing the pumping,
these two dimensionless parameters determine the behavior of the system. Note
that if one takes the standard definition for small-scale Rossby number
$\mathrm{Ro}_{f}=\epsilon^{1/3}k_{f}^{2/3}/2\Omega_{0}$, then the large-scale
Rossby number $\mathrm{Ro}_{\scriptscriptstyle
R}=\mathrm{Ro}_{f}\sqrt{\mathrm{Re}_{f}}$.

##  IV Statistics of turbulent pulsations

Equation (5) contains a component $\langle u^{\varphi}u^{{r}}\rangle$ of the
Reynolds stress tensor, which can be calculated if we assume that the gradient
of the velocity ${\bm{u}}$ in the pulsations is small in comparison with the
gradient of the mean current $\Sigma$. Then the evolution of velocity
pulsations can be described in the linear approximation. In this
approximation, equation (1) written in the Cartesian coordinate system
$\\{\xi,\eta,z\\}$, which was introduced above, takes the form

$(\partial_{t}-\Sigma\,\eta\,\partial_{\xi})\bm{u}-u^{\eta}\,\Sigma\,{\bm{e}}_{\xi}+2\,\Omega\left[{\bm{e}}_{z},\bm{u}\right]=-\nabla\tilde{p}+\nu\Delta\bm{u}+\bm{f},$
(13)

where the local angular velocity $\Omega$ of fluid rotation was defined in
expression (7). The random external force is assumed to be statistically
isotropic and shortly correlated in time. Its pair correlation function is

$\displaystyle\langle
f^{i}_{\bm{k}}(t_{1})f^{j}_{\bm{q}}(t_{2})\rangle=(2\pi)^{3}\left(\delta^{ij}-\frac{k^{i}k^{j}}{{\bm{k}}^{2}}\right)\times$
(14) $\displaystyle\epsilon\delta(t_{1}-t_{2})\delta(\bm{k}+\bm{q})\chi(k)$

in the Fourier space. The spatial profile $\chi(k)$ tends to zero at $k\to 0$
and $k\to\infty$, so the force does not excite very large and very small
eddies; it reaches maximum at $k\sim k_{f}$. For further calculations, we
choose the Gaussian spatial profile

$\chi(k)=\dfrac{16\pi^{3/2}}{3k_{f}^{5}}{k}^{2}\exp\left[-\dfrac{{k}^{2}}{k_{f}^{2}}\right].$
(15)

Figure 4: Numerical data for dimensionless tangent Reynolds stress
$F(\mathrm{Ro},\gamma)=\Sigma\langle u^{\varphi}u^{{r}}\rangle/\epsilon$ at
negative Rossby numbers $\mathrm{Ro}=\Sigma/(2\Omega)<0$, which correspond to
the cyclone. The dashed line shows that the curve $q=1$ coincides with the
curve $F=0$, see equation (12).

The problem (13) can be solved analytically [27] or numerically if the
function ${\bm{f}}(t,{\bm{r}})$ is known and the parameters $\Omega$ and
$\Sigma$ do not depend on time, see Appendix A. The statistical properties of
${\bm{u}}$ are determined by two dimensionless parameters: the local Rossby
number $\mathrm{Ro}$ for the large-scale flow and the local parameter $\gamma$

$\mathrm{Ro}=\frac{\Sigma}{2\Omega},\qquad\gamma=\frac{2\nu
k_{f}^{2}}{|\Sigma|}.$ (16)

${\mathrm{Ro}}$ is the local version of the global Rossby number
${\mathrm{Ro}}_{\scriptscriptstyle R}$. If one admits
$\Sigma\sim\sqrt{\epsilon/\nu}$, then $\gamma\sim({\mathrm{Re}}_{f})^{-2/3}\ll
1$.

For the component $\langle u^{r}u^{\varphi}\rangle$ of the Reynolds stress
tensor, we obtain

$\langle u^{r}u^{\varphi}\rangle=-\langle
u^{\eta}u^{\xi}\rangle=\frac{\epsilon}{\Sigma}F(\mathrm{Ro},\gamma),$ (17)

where the function $F(\mathrm{Ro},\gamma)$ should be evaluated numerically
based on equation (45) with $i=\eta$, $j=\xi$. The results are shown in Fig. 3
and in Fig. 4 for positive (anticyclone) and negative (cyclone) Rossby numbers
$\mathrm{Ro}$ respectively. Note that in both cases $F$ has
$\mathrm{Ro}$-independent limit at small Rossby numbers [20].

To verify our hypothesis concerning the absence of correlation between
pressure and velocity we independently compute the viscous dissipation rate
$q$ of the kinetic energy of turbulent pulsations, see equation (46).
According to energy balance (12), one finds

$q=1-F,$ (18)

where $F>0$ due to equation (6). As both $q$ and $F$ are analytical functions
of $\mathrm{Ro}$ and $\gamma$, equality (18) should also hold for negative
values of $F$, although $F<0$ cannot correspond to the vortex interior. Our
numerical calculations confirm equation (18). In particular, the dashed lines
in Fig. 3a and in Fig. 4 demonstrate that the curve $q=1$ coincides with the
curve $F=0$. Also note, that the dimensionless function $F$ obeys the
restriction $F<1$, since the viscous dissipation rate $q$ is always positive.

When Rossby number $\operatorname{Ro}$ and parameter $\gamma$ are small, the
function $F$ approaches the unity, see Ref. 20, Eq.(26). One can find the
asymptotic behavior

$q=1-F\approx q_{1/3}\gamma^{1/3},$ (19)

where the value of $q_{1/3}\approx 1.8$ was calculated numerically. Note that
the value of the constant $q_{1/3}$ depends on the specific choice of the
spatial pumping profile $\chi(k)$, see equation (15).

##  V Vortex velocity profile

Figure 5: Dependence of the cyclone velocity profile on the parameters
$\mathrm{Re}_{f}$ (left panel, $\operatorname{Ro}_{R}=0.5$) and
$\operatorname{Ro}_{R}$ (right panel, $\mathrm{Re}_{f}=20$). (a,b) Azimuthal
velocity of the cyclone $U$ and (c,d) parameter $\Sigma=r\partial_{r}(U/r)$ as
a function of the distance $r/R$ to the vortex axis,
$\Sigma_{0}=\sqrt{\epsilon/\nu}$. (e,f) Variation of the local parameters
$\operatorname{Ro}$ and $\gamma$ inside the cyclone. The markers indicate the
points corresponding to $r/R=1$. The area of curves downward from the marker
corresponds to the movement to the vortex axis, in the opposite direction – to
its periphery (shown with dashed color lines in the left panel). The black
dashed line corresponds to $F/\operatorname{Ro}=1/2$ (one-half takes into
account the anisotropy of the velocity statistics) and above the line the
linear analysis of turbulent pulsations becomes inapplicable.

In this section we analyze the radial profile of the mean velocity $U(r)$ in
the coherent vortex. First, let us rewrite equation (5) in the dimensionless
form. For this purpose, we introduce dimensionless variables: local shear rate
$\sigma=\Sigma/\Sigma_{0}$, local rotation rate $\omega=2\Omega/\Sigma_{0}$
and the distance $\rho=r/R$ to the vortex axis, where
$\Sigma_{0}=\sqrt{\epsilon/\nu}$. Then we obtain

$\displaystyle\sigma^{2}={F\left(\frac{\sigma}{\omega},\frac{2}{({\mathrm{Re}}_{f})^{3/2}|\sigma|}\right)},\qquad\sigma=\rho\partial_{\rho}\omega/2,$
(20)

and this essentially nonlinear first-order differential equation involving
unknown function $\omega(\rho)$ should be supplemented by the boundary
condition

$\omega\big{|}_{\rho=1}=\frac{1}{\mathrm{Ro}_{\scriptscriptstyle R}},$ (21)

which means that $U(R)=0$ according to expression (4). Note that equation (20)
contains only one dimensionless parameter ${\mathrm{Re}}_{f}$, although the
vortex structure generally depends on ${\mathrm{Re}}_{f}$ and
$\mathrm{Ro}_{\scriptscriptstyle R}$. Thus, all solutions of equation (20)
belong to a one-parameter family of curves, and the choice of a particular
curve is determined by ${\mathrm{Re}}_{f}$. The value of
$\mathrm{Ro}_{\scriptscriptstyle R}$ defines the position on the curve
corresponding to the boundary or ’initial’ condition (21), see Figs. 5,6 and
discussion below.

###  V.1 Cyclones

In the case of a cyclone, the negative root of equation (20) should be chosen,

$\displaystyle\sigma=-\sqrt{F\left(\frac{\sigma}{\omega},\frac{2}{({\mathrm{Re}}_{f})^{3/2}\cdot|\sigma|}\right)},\quad\sigma=\rho\partial_{\rho}\omega/2.$
(22)

We then numerically solve this differential-algebraic equation supplemented
with the initial condition (21). The results are presented in Fig. 5.

The left panel shows how the radial velocity profile depends on the forcing
Reynolds number $Re_{f}$ at a fixed global Rossby number
$\operatorname{Ro}_{R}=0.5$. Fig. 5a demonstrates that with an increase in the
Reynolds number, the velocity in the vortex increases, and the radial profile
is quite close to the linear-logarithmic one. This corresponds to the almost
complete independence of the local shear rate $\Sigma(r)$ on the distance $r$
to the vortex axis, see Fig. 5c and equation (4). The dashed lines in Fig. 5a
show the dependence $U=-\Sigma^{*}r\ln R/r$, where $\Sigma^{*}=\Sigma(R)$. The
angular velocity $U(r)/r$ associated with the self-rotation of the vortex
increases when approaching the cyclone axis; therefore, the absolute value of
the local Rossby number $\operatorname{Ro}$ decreases. This is illustrated in
Fig. 5e, where we show the phase portrait of equation (22) on the background
of a color map corresponding to the values of $F(\operatorname{Ro},\gamma)$
obtained in Section IV. By changing the Reynolds number $\mathrm{Re}_{f}$, we
move from one curve to another, and the marker corresponds to the coordinate
$r=R$. The movement towards the axis of the cyclone in this diagram
corresponds to the movement downward along the curves.

The right panel shows how the behavior of the system changes if we vary the
global Rossby number $\operatorname{Ro}_{R}$ at a fixed value of the Reynolds
number $\mathrm{Re}_{f}=20$. In this case, we remain on the same phase curve,
but the position of the initial point $r=R$ changes, see Fig. 5f. The larger
the Rossby number $\operatorname{Ro}_{R}$, the higher we are on this curve and
the more significant the difference between the radial velocity profile and
the linear-logarithmic prediction, since the local parameter $\Sigma(r)$
begins to depend on $r$, see Figs. 5b and 5d. Note that for the values
$\operatorname{Ro}_{R}=1$ and $\mathrm{Re}_{f}=20$, the value
$\operatorname{Ro}_{f}\approx 0.2\lesssim 1$, which means that the fluid is
still rotating rapidly.

In the limit ${\mathrm{Ro}}_{R}\ll 1$ and $\mathrm{Re}_{f}\gg 1$ one can solve
equation (22) analytically using the asymptotic expression (19). The result is
a linear-logarithmic velocity profile $U=-\tilde{\sigma}\Sigma_{0}r\ln R/r$,
where the parameter $\tilde{\sigma}$ is the minimum root of the equation

$\tilde{\sigma}=-\sqrt{1-\dfrac{2^{1/3}q_{1/3}}{Re_{f}^{1/2}|\tilde{\sigma}|^{1/3}}}.$
(23)

In particular, for $\mathrm{Re}_{f}=20$, one finds
$\tilde{\sigma}\approx-0.64$, which is in reasonable agreement with the
corresponding result presented in Fig. 5d for $\operatorname{Ro}_{R}=0.1$.

The evolution of turbulent pulsations cannot be considered in the local
approximation (adopted in Section IV) near the axis of the vortex; therefore,
our analysis should be limited to the region $r>r_{c}$, where $r_{c}$ defines
the radius of the vortex core. Inside the core, inertial waves do not have
time to transfer their energy to the large-scale flow before they leave the
core. The group velocity of the waves is $v_{g}\sim\Omega(r_{c})/k_{f}$ and
the time needed for the energy transfer is $\sim 1/\Sigma(r_{c})$. Thus, the
core radius is determined by the equation [20]

${\mathrm{Ro}}(r_{c})\cdot k_{f}r_{c}\sim 1.$ (24)

Inside the core, the inertial waves are incoherently affected by inhomogeneous
large-scale velocity field during their propagation in space and thus they do
not produce considerable Reynolds stress $\langle u^{r}u^{\varphi}\rangle$.
Therefore, one expects approximately rigid-body rotation of the fluid and the
velocity profile has a linear dependence on $r$ in the region $r<r_{c}$.

Besides, the developed theory is not applicable in the region where the
evolution of turbulent pulsations goes beyond the linear approximation. This
occurs near the curve $F=0$, since in this case there is practically no energy
transfer from fluctuations to the coherent vortex, see equation (12). To
establish the applicability condition of our consideration, we must compare
the nonlinear interaction rate with the shear rate $|\Sigma|$. The nonlinear
interaction rate $\Gamma_{\\!nl}$ stems from the nonlinear term
$({\bm{u}}\nabla){\bm{u}}$ omitted in equation (13). Since the Rossby number
is less than unity, the rate is diminished by the inertia wave oscillations
and its upper estimation is
$\Gamma_{\\!nl}\lesssim\langle(\nabla{\bm{u}})^{2}\rangle/2\Omega$, see, e.g.,
Ref. 28. Some overestimation is due to the fact that the term
$\langle(\nabla{\bm{u}})^{2}\rangle$ takes into account, inter alia, the
spatial derivative normal to the ${\bm{u}}$-direction, whereas this derivative
is absent in $({\bm{u}}\nabla){\bm{u}}$. The difference affects the result
since the largest component is $u^{\xi}$ and the correlation length reaches a
maximum in $\xi$-direction in the presence of a shear flow, see, e.g., Ref.
29. Next, in the linear approximation
$\langle(\nabla{\bm{u}})^{2}\rangle\sim(1-F)\epsilon/\nu$ and using equation
(26), we can rewrite the inequality $\Gamma_{\\!nl}\lesssim|\Sigma|$ as

$F/\mathrm{Ro}\gtrsim 1.$ (25)

The corresponding line is shown in Figs. 5e and 5f (the factor $1/2$ instead
of unity takes into account the aforementioned overestimation). Above this
line, the applicability condition (25) is violated, which means the
development of three-dimensional turbulence. In this regime, the Reynolds
stress tensor component $\langle u^{r}u^{\varphi}\rangle$ has the opposite
sign, which makes the existence of coherent vortices impossible. Thus, the
$\textquoteleft R{\text{\textquoteright}}$-point corresponding to the cyclone
radius $r=R$ and indicated with the marker should be located below the black
dashed curve in Figs. 5e,f. This imposes restrictions on the global flow
parameters ${\mathrm{Re}}_{f},{\mathrm{Ro}}_{\scriptscriptstyle R}$ at which
the cyclone can exist.

###  V.2 Anticyclones

Figure 6: Dependence of the anticyclone velocity profile on the parameters
$\mathrm{Re}_{f}$ (left panel, $\operatorname{Ro}_{R}=0.5$) and
$\operatorname{Ro}_{R}$ (right panel, $\mathrm{Re}_{f}=20$). (a,b) Azimuthal
velocity of the anticyclone $U$ and (c,d) parameter
$\Sigma=r\partial_{r}(U/r)$ as a function of the distance $r/R$ to the vortex
axis, $\Sigma_{0}=\sqrt{\epsilon/\nu}$. (e,f) Variation of the local
parameters $\operatorname{Ro}$ and $\gamma$ inside the vortex. The markers
indicate the points corresponding to $r/R=1$. The area of curves right-up from
the marker corresponds to the movement to the vortex axis, in the opposite
direction – to its periphery (shown with dashed color lines in the left
panel). The black dashed line corresponds to $F/\operatorname{Ro}=1/2$ (one-
half takes into account the anisotropy of the velocity statistics) and above
the line the linear analysis of turbulent pulsations becomes inapplicable.

In the case of an anticyclone, the positive root of equation (20) should be
chosen,

$\displaystyle\sigma=\sqrt{F\left(\frac{\sigma}{\omega},\frac{2}{({\mathrm{Re}}_{f})^{3/2}\sigma}\right)},\quad\sigma=\rho\partial_{\rho}\omega/2,$
(26)

and a similar analysis of solutions to this equation is presented in Fig. 6.
In comparison with cyclones under the same values of global parameters
$\mathrm{Re}_{f}$ and $\operatorname{Ro}_{R}$, it should be noted that the
maximum value of the velocity is lower in absolute value and the position of
the maximum is farther from the axis of the vortex.

The self-rotation of the anticyclone $U(r)/r$ is directed opposite to the
external rotation $\Omega_{0}$, and now the local Rossby number
$\operatorname{Ro}$ increases when approaching the vortex axis. Accordingly,
the movement along the phase curves shown in Figs. 6e and 6f occurs in the
opposite direction compared to cyclones. This leads to a stronger dependence
of the local shear rate $\Sigma(r)$ on the distance to the vortex axis $r$ at
the same values of the global parameters $\mathrm{Re}_{f}$ and
$\operatorname{Ro}_{R}$, which is expressed in a noticeable difference between
the velocity profile $U(r)$ and the linear-logarithmic prediction, see Figs.
6a-d. Good agreement is observed only for a small value of the Rossby number
$\operatorname{Ro}_{R}=0.1$. In the limit $\operatorname{Ro}_{R}\ll 1$ and
$\mathrm{Re}_{f}\gg 1$ the profiles of cyclones and anticyclones coincide up
to sign, so the previous asymptotic expression (23) can be used.

When moving along the phase curves, we inevitably come to the region
$F\lesssim\operatorname{Ro}$, where the nonlinear interaction of fluctuations
becomes essential, which should lead to the destruction of the anticyclone due
to the direct energy cascade, see equation (25). The only way to prevent this
is to require that the size $r_{c}$ of the vortex core (24) be large enough so
that the point $r=r_{c}$ lies outside the region. As the shear rate $\Sigma$
vanishes inside the core, we should impose one more requirement to avoid the
development of three-dimensional turbulence inside the core. It is suppressed
if the rotation inside the core is large enough,
${\mathrm{Ro}}_{\epsilon}=(\epsilon k_{f}^{2})^{1/3}/(2\Omega(r_{c}))<1$
(number ${\mathrm{Ro}}_{\epsilon}$ is the local version of
${\mathrm{Ro}}_{f}$). Then the external force ${\bm{f}}$ excites primarily
inertia waves, which then spread away from the core. The requirement is always
satisfied if the condition (25) is met, since it can be rewritten in a form
$F/(\gamma\cdot{\mathrm{Ro}}^{3})>1$ at $r=r_{c}$. Thus, it is necessary to
take care only about the relation (25), which should be fulfilled at
$r=r_{c}$.

Next, it follows from equation (26) that the ratio between the vortex and the
core radii is

$\frac{R}{r_{c}}=\exp\left(\int_{R}^{r_{c}}\frac{\mathrm{d}\ln(\gamma\cdot{\mathrm{Ro}})}{2{\mathrm{Ro}}}\right),$
(27)

where the integral should be taken along the phase curve. The location of
$\textquoteleft\\!R\text{\textquoteright}$-point on the phase plane is
determined by the global parameters ${\mathrm{Ro}}_{\scriptscriptstyle R}$ and
${\mathrm{Re}}_{f}$, see the marked points in Figs. 6e,f. Now suppose for
simplicity that $\textquoteleft r_{c}\text{\textquoteright}$-point corresponds
to $\mathrm{Ro}(r_{c})\ll 1$. Then the parameter $\gamma$ remains
approximately constant along the phase curve and we can perform integration in
equation (27)

$2\ln\frac{R}{r_{c}}=\frac{1}{{\mathrm{Ro}}(R)}-\frac{1}{{\mathrm{Ro}}(r_{c})}<\frac{1}{{\mathrm{Ro}}(R)}.$
(28)

Condition (21) implies $\mathrm{Ro}(R)=\sigma(R)\cdot\operatorname{Ro}_{R}$,
and in the limit $\operatorname{Ro}_{R}\ll 1$ and $\mathrm{Re}_{f}\gg 1$ one
finds $\sigma(R)=1$, see Ref. 20. The criterion in this form was mentioned in
the Introduction.

In the general case, one needs to find numerically the value of
$\rho_{nl}=r_{nl}/R$ that corresponds to the intersection between the black
dashed line (25) and the phase curve shown in Figs. 6e,f. Note that this value
will depend only on the global flow parameters $\mathrm{Re}_{f}$ and
$\operatorname{Ro}_{R}$, i.e.
$\rho_{nl}=\rho_{nl}(\mathrm{Re}_{f},\operatorname{Ro}_{R})$. Next, as we
discussed above, $r_{c}\geq r_{nl}$ and using equation (24) we arrive at

$R<R_{max}=\dfrac{1}{k_{f}}\max_{\rho_{nl}\leq\rho\leq
1}\dfrac{1}{\rho\operatorname{Ro}(\rho)}.$ (29)

The value of $R_{max}$ should be calculated numerically for specific values of
$\operatorname{Ro}_{R}$ and $\mathrm{Re}_{f}$. For the parameters under study,
the maximum is reached at point $\rho=\rho_{nl}$, which corresponds to
$r_{c}=r_{nl}$. In particular, for $\mathrm{Re}_{f}=20$ and
$\operatorname{Ro}_{R}=0.5$ one finds $k_{f}R_{max}\sim 5$ and for
$\mathrm{Re}_{f}=20$ and $\operatorname{Ro}_{R}=0.1$ we obtain
$k_{f}R_{max}\sim 2.6\times 10^{3}$. If the size of the system $L$
substantially exceeds $R_{max}$ then we expect a difference between the sizes
of cyclones and anticyclones in the condensate, since the size of cyclones is
limited only by the system size.

##  VI Concluding remarks

To summarize, we considered a coherent vortex in a condensate of turbulent
flow in a rotating fluid when Rossby number $\operatorname{Ro}_{R}\lesssim 1$
is small and Reynolds number $\mathrm{Re}_{f}\gtrsim 10$ is large. In the
limit $\operatorname{Ro}_{R}\ll 1$, the radial profiles of the mean azimuthal
velocity in cyclones and anticyclones are equal to each other up to a sign,
$U(r)=\mp\tilde{\sigma}\Sigma_{0}r\ln R/r$, where $R$ is the vortex radius,
$\Sigma_{0}=\sqrt{\epsilon/\nu}$, and $\tilde{\sigma}$ is determined by
equation (23). At finite Rossby numbers, we found some differences between
them. Qualitatively, the maximum of $|U(r)|$ shifts to the center inside
cyclones and to the periphery inside anticyclones, and the maximum value of
$|U(r)|$ is greater inside cyclones than inside anticyclones at the same
global parameters of the flow.

Coherent vortices can exist only if condition (25) is met at the point $r=R$,
which imposes restrictions on the allowable range of $\mathrm{Re}_{f}$ and
$\operatorname{Ro}_{R}$. The induced rotation $U(r)/r$ of the cyclone is added
to the external rotation $\Omega_{0}$, so the fulfillment of inequality (25)
at a point $r=R$ guarantees that it will be satisfied at any other point $r<R$
inside the vortex. In the case of anticyclones, the induced rotation
compensates the external rotation, so the movement along the phase curves in
Figs. 6e,f occurs in the opposite direction. This imposes an additional
condition on the existence of anticyclones – inequality (25) must be fulfilled
on the boundary of the vortex core, which size $r_{c}$ is mainly determined by
the forcing scale, see equation (24). This condition leads to the fact that
the size $R$ of the anticyclone is limited by $R_{max}$ that depends on the
global flow parameters $\mathrm{Re}_{f}$ and $\operatorname{Ro}_{R}$, see
equation (29). We predict that if the size of the system $L$ significantly
exceeds $R_{max}$, then the sizes of cyclones and anticyclones in the
condensate may also differ significantly. One should remember that expression
(25) is only an estimate, and the exact criterion must be obtained either by
numerical simulation or experimentally. Besides, our theory is applicable only
if the local approximation is valid for describing small-scale turbulent
pulsations, when ${\mathrm{Ro}}_{\scriptscriptstyle R}\cdot k_{f}L\gg 1$.

As for comparing the developed theory with experiments, we should emphasize
that our approach takes into account only bulk viscous dissipation without any
influence of boundaries. In experiments, the influence of the lower and upper
boundaries can be significant. Indeed, the thickness of the boundary layer is
$\delta\sim\sqrt{\nu/\Omega_{0}}$, and then the scale $r_{\alpha}$ at which
the bottom friction becomes comparable with the bulk viscosity is
$r_{\alpha}\sim(\nu/\Omega_{0})^{1/4}H^{1/2}$, where $H$ is the height of the
vessel. For typical values $H=50\,\text{cm}$,
$\nu=0.01\,\text{cm}^{2}/\text{s}$, $\Omega_{0}=3\,\text{s}^{-1}$, see Refs.
7, 11, 1, one finds $r_{\alpha}\sim 2\,\text{cm}$ that is of the order or
smaller than the radii of the observed columnar vortices. Thus, the
predictions of our theory can be directly applied only to the results of
numerical simulations with periodic boundary conditions, while the description
of the experiments requires further research.

###### Acknowledgements.

We are grateful to V. Lebedev and I. Kolokolov for valuable discussions. This
work was supported by the Russian Science Foundation, Grant No. 20-12-00383.
SV and IV were supported by Grant No. 19-1-2-46-1 of Foundation for the
Advancement of Theoretical Physics and Mathematics “BASIS.”

## Data AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Dynamics and statistical properties of turbulent pulsations

It is convenient to solve equation (13) in the Fourier space

$\left(\partial_{t}+\Sigma\,k_{\xi}\,\partial_{k_{\eta}}\right)\bm{u}_{\bm{k}}-u^{\eta}_{\bm{k}}\,\Sigma\,{\bm{e}}_{\xi}+2\,\Omega\left[{\bm{e}}_{z},\bm{u}_{\bm{k}}\right]=-i\bm{k}\tilde{p}_{\bm{k}}-\nu{\bm{k}}^{2}\bm{u}_{\bm{k}}+\bm{f}_{\bm{k}}$
(30)

using the method of characteristics. Then one has to solve the ordinary
differential equation

$\dot{\bm{u}}_{\bm{K}}=u^{\eta}_{\bm{K}}\,\Sigma\,{\bm{e}}_{\xi}-2\,\Omega\left[{\bm{e}}_{z},\bm{u}_{\bm{K}}\right]-i\bm{K}\tilde{p}_{\bm{K}}-\nu{\bm{K}}^{2}\bm{u}_{\bm{K}}+\bm{f}_{\bm{K}}$
(31)

along the characteristics

$\bm{K}(t)=\\{k_{\xi},\,K_{\eta}(t),\,k_{z}\\},\quad
K_{\eta}(t)=k_{\eta}+\Sigma\,t\,k_{\xi}.$ (32)

To exclude the pressure term, we note that
$\mathrm{d}\left(\bm{K},\,{\bm{u}}_{\bm{K}}\right)/\mathrm{d}t=0$, so the
projection of (31) onto ${\bm{K}}$ gives

$\tilde{p}_{\bm{K}}=\frac{2i}{{\bm{K}}^{2}}\biggl{(}\Omega
k_{\eta}u^{\xi}_{\bm{K}}-\left(\Omega+\Sigma\right)k_{\xi}u^{\eta}_{\bm{K}}\biggr{)}.$
(33)

Next, we substitute expression (33) into equation (31) and rewrite it in a
matrix form

$\dot{\bm{u}}_{\bm{K}}=\hat{C}\bm{u}_{\bm{K}}-\nu{\bm{K}}^{2}\bm{u}_{\bm{K}}+\bm{f}_{\bm{K}},\qquad\hat{C}=\left(\begin{array}[]{ccc}2\,\Omega\frac{k_{\xi}K_{\eta}}{K^{2}}&\Sigma\left(1-2\frac{k_{\xi}^{2}}{K^{2}}\right)+2\,\Omega\left(1-\frac{k_{\xi}^{2}}{K^{2}}\right)&0\\\
-2\,\Omega\left(1-\frac{K_{\eta}^{2}}{K^{2}}\right)&-(\Omega+\Sigma)\frac{2k_{\xi}K_{\eta}}{K^{2}}&0\\\
2\,\Omega\frac{K_{\eta}k_{z}}{K^{2}}&-(\Omega+\Sigma)\frac{2k_{\xi}k_{z}}{K^{2}}&0\end{array}\right),$
(34)

where ${\bm{u}}_{\bm{K}}$ is assumed to be a column vector. Due to the
incompressibility condition, the velocity vector ${\bm{u}}_{\bm{K}}$ depends
only on two independent scalar functions. In order to derive the equations on
these two functions, we follow Ref. [27] and we choose basis
$\left\\{\bm{e}_{1},\,\bm{e}_{2},\,\bm{e}_{3}\right\\}$

$\bm{e}_{3}=\frac{\bm{K}}{K},\quad\bm{e}_{1}=\frac{\left[{\bf
e}_{\eta},\bm{e}_{3}\right]}{\left|\left[{\bf
e}_{\eta},\bm{e}_{3}\right]\right|},\quad\bm{e}_{2}=\left[\bm{e}_{3},\bm{e}_{1}\right].$
(35)

The expansions of vectors $\bm{u}_{\bm{K}}$ and $\bm{f}_{\bm{K}}$ in this
basis are

$\bm{u}_{\bm{K}}=\sum_{\alpha}\psi_{\bm{K}}^{\alpha}\bm{e}_{\alpha},\quad\bm{f}_{\bm{K}}=\sum_{\alpha}\phi_{\bm{K}}^{\alpha}\bm{e}_{\alpha},\qquad\alpha,\beta=\\{1,2\\},$
(36)

but one should remember that the basis (35) is time dependent when
substituting the expansion (36) into the evolution equation (34). The result
of the corresponding substitution is

$\frac{d{}}{dt}\Psi_{\bm{K}}=\hat{H}\Psi_{\bm{K}}-\nu\mathbf{K}^{2}\Psi_{\bm{K}}+\Phi_{\bm{K}},\quad\Psi_{\bm{K}}=\begin{pmatrix}\psi^{1}_{\bm{K}}\\\
\psi^{2}_{\bm{K}}\\\
\end{pmatrix},\quad\Phi_{\bm{K}}=\begin{pmatrix}\phi^{1}_{\bm{K}}\\\
\phi^{2}_{\bm{K}}\\\
\end{pmatrix},\quad\hat{H}=\frac{1}{K^{2}}\begin{pmatrix}0&(\Sigma+2\,\Omega)Kk_{z}\\\
\\\ -2\,\Omega Kk_{z}&-\Sigma k_{\xi}K_{\eta}\end{pmatrix}.$ (37)

The pair correlation function of the external force (14) in the basis
$\left\\{\bm{e}_{1},\,\bm{e}_{2},\,\bm{e}_{3}\right\\}$ is equal to

$\langle\phi^{\alpha}_{\bm{k}}(t)\,\phi^{\beta}_{\bm{q}}(t^{\prime})\rangle=(2\pi)^{3}\delta({\bm{k}}+{\bm{q}})\big{(}\delta^{\alpha
1}\delta^{\beta 1}-\delta^{\alpha 2}\delta^{\beta
2}\big{)}\epsilon\delta(t-t^{\prime})\chi(k).$ (38)

The alternating sign in the tensor structure of equation (38) reflects the
fact that the basis (35) changes its orientation when ${\bm{k}}\to-{\bm{k}}$.

The solution of equation (37) can be represented in the form

$\Psi_{\bm{K}(t)}(t)=\int\limits_{-\infty}^{t}\mathrm{d}t^{\prime}\exp\left[-\nu\int\limits_{t^{\prime}}^{t}\mathrm{d}t^{\prime\prime}\,\bm{K}^{2}(t^{\prime\prime})\right]\hat{Q}(t,\,t^{\prime})\,\Phi_{\bm{K}(t^{\prime})}(t^{\prime}),$
(39)

where the evolution matrix $\hat{Q}(t,t^{\prime})$ obeys the equation

$\frac{d{}}{dt}\,\hat{Q}(t,t^{\prime})=\hat{H}(t)\hat{Q}(t,t^{\prime}),\qquad\hat{Q}(t^{\prime},t^{\prime})=\hat{1}.$
(40)

However, for numerical calculations it will be more convenient to set the
final moment and then the evolution matrix $\hat{Q}$ also obeys the backward
equation

$\frac{d}{dt^{\prime}}\,\hat{Q}(t,t^{\prime})=-\hat{Q}(t,t^{\prime})\,\hat{H}(t^{\prime}),\qquad\hat{Q}(t,\,t)=\hat{1}.$
(41)

Let us note that the matrix $\hat{H}$ (37) depends only on the direction of
wavevector ${\bm{k}}$, but not on its absolute value. Therefore, the same is
true for the evolution matrix $\hat{Q}$.

Now we are ready to proceed to the direct calculation of the Reynolds stress
tensor. According to the introduced definitions,

$\langle
u^{i}(t)u^{j}(t)\rangle=\int\dfrac{d^{3}\bm{k}}{(2\pi)^{3}}\int\dfrac{d^{3}\bm{s}}{(2\pi)^{3}}\langle
u^{i}_{\bm{K}}(t)u^{j}_{\bm{S}}(t)\rangle,\quad\langle
u^{i}_{\bm{K}}(t)u^{j}_{\bm{S}}(t)\rangle=\sum_{\alpha,\beta}\,\langle\psi^{\alpha}_{\bm{K}}(t)\psi^{\beta}_{\bm{S}}(t)\rangle
e^{i}_{\alpha}(\bm{K}(t))e^{j}_{\beta}(\bm{S}(t)),$ (42)

and using expression (39) we obtain

$\displaystyle\langle\psi^{\alpha}_{\bm{K}}(t)\psi^{\beta}_{\bm{S}}(t)\rangle=\int\limits_{-\infty}^{t}dt_{1}dt_{2}\exp\left[-\nu\left(\int\limits_{t_{1}}^{t}dt^{\prime}_{1}\,\bm{K}^{2}(t^{\prime}_{1})+\int\limits_{t_{2}}^{t}dt^{\prime}_{2}\,\bm{S}^{2}(t^{\prime}_{2})\right)\right]\times$
(43)
$\displaystyle\sum_{\sigma,\rho}\,Q^{\alpha\sigma}\bigl{(}\bm{K}(t),\bm{K}(t_{1})\bigr{)}Q^{\beta\rho}\bigl{(}\bm{S}(t),\bm{S}(t_{2})\bigr{)}\langle\phi^{\sigma}_{\bm{K}}(t_{1})\phi^{\rho}_{\bm{S}}(t_{2})\rangle.$

Since the statistics of the turbulent flow is assumed to be stationary in
time, we can set $t=0$ in equation (43). Further calculations, taking into
account relation (38), are rather straightforward and we omit them for
brevity. We present the final results using the introduced parameters (16).
Besides, we also introduce the dimensionless time
$\tau=-|\Sigma|t^{\prime}>0$, where $t^{\prime}$ appears in equation (41), and
$\varsigma=-\mathop{\mathrm{sign}}\Sigma$, so that $\varsigma=1$ and
$\varsigma=-1$ correspond to cyclone and anticyclone respectively. We also
define

$\displaystyle\lambda(\theta,\,\varphi,\,\tau,\,\varsigma)=1+\varsigma\tau\sin
2\varphi\sin^{2}\theta+(\tau\cos\varphi\sin\theta)^{2},\qquad\Lambda(\theta,\,\varphi,\,\tau,\,\varsigma)=\int\limits_{0}^{\tau}d\tau^{\prime}\lambda(\theta,\,\varphi,\,\tau^{\prime},\,\varsigma),$
(44)

and then the Reynolds stress tensor components are equal to

$\langle
u^{i}u^{j}\rangle=\dfrac{\epsilon}{\pi|\Sigma|}\,\int\limits_{0}^{\pi/2}d\theta\int\limits_{0}^{\pi}d\varphi\int\limits_{0}^{+\infty}d\tau\,\frac{\lambda\,\sin\theta}{\left(\lambda+\gamma\Lambda\right)^{-5/2}}\sum_{\alpha,\beta,\delta=1,2}e_{\alpha}^{i}\,e_{\beta}^{j}\,Q^{\alpha\delta}\,Q^{\beta\delta},$
(45)

where we used the Gaussian spatial profile (15) for the pumping. In a similar
way, for the dimensionless viscous dissipation rate we obtain

$q=\dfrac{5}{2\pi}\int\limits_{0}^{\pi/2}d\theta\int\limits_{0}^{\pi}d\varphi\int\limits_{0}^{+\infty}d\tau\,\frac{\sin\theta}{\left(\lambda+\gamma\Lambda\right)^{-7/2}}\sum_{\alpha,\beta=1,2}\big{(}Q^{\alpha\beta}\big{)}^{2}.$
(46)

The evolution matrix $Q^{\alpha\beta}$ should be found by solving equation
(41), which takes the form

$\dfrac{\partial{}}{\partial{\tau}}\,\hat{Q}(\tau)=\hat{Q}(\tau)\begin{pmatrix}0&-\varsigma\left(1+\dfrac{1}{\mathrm{Ro}}\right)\dfrac{\cos\theta}{\sqrt{\lambda(\tau)}}\\\\[5.0pt]
\dfrac{\varsigma}{\mathrm{Ro}}\dfrac{\cos\theta}{\sqrt{\lambda(\tau)}}&\dfrac{\cos\varphi\sin^{2}\theta(\tau\cos\varphi+\varsigma\sin\varphi)}{\lambda(\tau)}\end{pmatrix},\qquad\hat{Q}(0)=\hat{1}.$
(47)

## References

  * Godeferd and Moisy [2015] F. S. Godeferd and F. Moisy, “Structure and dynamics of rotating turbulence: a review of recent experimental and numerical results,” Applied Mechanics Reviews 67 (2015).
  * Proudman [1916] J. Proudman, “On the motion of solids in a liquid possessing vorticity,” Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 92, 408–424 (1916).
  * Taylor [1917] G. I. Taylor, “Motion of solids in fluids when the flow is not irrotational,” Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 93, 99–113 (1917).
  * Boffetta and Ecke [2012] G. Boffetta and R. E. Ecke, “Two-dimensional turbulence,” Annual Review of Fluid Mechanics 44, 427–451 (2012).
  * McEwan [1976] A. McEwan, “Angular momentum diffusion and the initiation of cyclones,” Nature 260, 126–128 (1976).
  * Hopfinger, Browand, and Gagne [1982] E. Hopfinger, F. Browand, and Y. Gagne, “Turbulence and waves in a rotating tank,” Journal of Fluid Mechanics 125, 505–534 (1982).
  * Ruppert-Felsot  _et al._ [2005] J. E. Ruppert-Felsot, O. Praud, E. Sharon, and H. L. Swinney, “Extraction of coherent structures in a rotating turbulent flow experiment,” Physical Review E 72, 016311 (2005).
  * Staplehurst, Davidson, and Dalziel [2008] P. Staplehurst, P. Davidson, and S. Dalziel, “Structure formation in homogeneous freely decaying rotating turbulence,” Journal of Fluid Mechanics 598, 81 (2008).
  * Moisy  _et al._ [2011] F. Moisy, C. Morize, M. Rabaud, and J. Sommeria, “Decay laws, anisotropy and cyclone-anticyclone asymmetry in decaying rotating turbulence,” Journal of Fluid Mechanics 666, 5 (2011).
  * Gallet  _et al._ [2014] B. Gallet, A. Campagne, P.-P. Cortet, and F. Moisy, “Scale-dependent cyclone-anticyclone asymmetry in a forced rotating turbulence experiment,” Physics of Fluids 26, 035108 (2014).
  * Boffetta  _et al._ [2020] G. Boffetta, F. Toselli, M. Manfrin, and S. Musacchio, “Cyclone–anticyclone asymmetry in rotating thin fluid layers,” Journal of Turbulence , 1–12 (2020).
  * Bartello, Métais, and Lesieur [1994] P. Bartello, O. Métais, and M. Lesieur, “Coherent structures in rotating three-dimensional turbulence,” Journal of Fluid Mechanics 273, 1–29 (1994).
  * Yeung and Zhou [1998] P. Yeung and Y. Zhou, “Numerical study of rotating turbulence with external forcing,” Physics of Fluids 10, 2895–2909 (1998).
  * Smith and Waleffe [1999] L. M. Smith and F. Waleffe, “Transfer of energy to two-dimensional large scales in forced, rotating three-dimensional turbulence,” Physics of Fluids 11, 1608–1622 (1999).
  * Yoshimatsu, Midorikawa, and Kaneda [2011] K. Yoshimatsu, M. Midorikawa, and Y. Kaneda, “Columnar eddy formation in freely decaying homogeneous rotating turbulence,” Journal of Fluid Mechanics 677, 154 (2011).
  * Biferale  _et al._ [2016] L. Biferale, F. Bonaccorso, I. M. Mazzitelli, M. A. van Hinsberg, A. S. Lanotte, S. Musacchio, P. Perlekar, and F. Toschi, “Coherent structures and extreme events in rotating multiphase turbulent flows,” Physical Review X 6, 041036 (2016).
  * Seshasayanan and Alexakis [2018] K. Seshasayanan and A. Alexakis, “Condensates in rotating turbulent flows,” Journal of Fluid Mechanics 841, 434–462 (2018).
  * Deusebio  _et al._ [2014] E. Deusebio, G. Boffetta, E. Lindborg, and S. Musacchio, “Dimensional transition in rotating turbulence,” Physical Review E 90, 023005 (2014).
  * Naso [2015] A. Naso, “Cyclone-anticyclone asymmetry and alignment statistics in homogeneous rotating turbulence,” Physics of Fluids 27, 035108 (2015).
  * Kolokolov, Ogorodnikov, and Vergeles [2020] I. Kolokolov, L. Ogorodnikov, and S. Vergeles, “Structure of coherent columnar vortices in three-dimensional rotating turbulent flow,” Physical Review Fluids 5, 034604 (2020).
  * Doludenko  _et al._ [2021] A. Doludenko, S. Fortova, I. Kolokolov, and V. Lebedev, “Coherent vortex in a spatially restricted two-dimensional turbulent flow in absence of bottom friction,” Physics of Fluids 33, 011704 (2021).
  * Smits, McKeon, and Marusic [2011] A. J. Smits, B. J. McKeon, and I. Marusic, “High–Reynolds number wall turbulence,” Annual Review of Fluid Mechanics 43, 353–375 (2011).
  * Balbus [2017] S. A. Balbus, “When is high Reynolds number shear flow not turbulent?” Journal of Fluid Mechanics 824, 1–4 (2017).
  * Nazarenko, Kevlahan, and Dubrulle [2000] S. Nazarenko, N.-R. Kevlahan, and B. Dubrulle, “Nonlinear RDT theory of near-wall turbulence,” Physica D: Nonlinear Phenomena 139, 158–176 (2000).
  * Falkovich [2016] G. Falkovich, “Interaction between mean flow and turbulence in two dimensions,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472, 20160287 (2016).
  * Laurie  _et al._ [2014] J. Laurie, G. Boffetta, G. Falkovich, I. Kolokolov, and V. Lebedev, “Universal profile of the vortex condensate in two-dimensional turbulence,” Physical Review Letters 113, 254503 (2014).
  * Salhi and Cambon [1997] A. Salhi and C. Cambon, “An analysis of rotating shear flow using linear theory and DNS and LES results,” Journal of Fluid Mechanics 347, 171–195 (1997).
  * Galtier [2003] S. Galtier, “Weak inertial-wave turbulence theory,” Physical Review E 68, 015301 (2003).
  * Champagne, Harris, and Corrsin [1970] F. Champagne, V. Harris, and S. Corrsin, “Experiments on nearly homogeneous turbulent shear flow,” Journal of Fluid Mechanics 41, 81–139 (1970).