1. Introduction
The inverse energy cascade of two-dimensional turbulence, with small-scale forcing and weak damping, and in a finite domain, leads to a condensation in which the energy is overwhelmingly concentrated in the lowest available modes of the system. In a square domain with periodic boundary, the mean state takes the form of a vortex dipole (Chertkov et al. Reference Chertkov, Connaughton, Kolokolov and Lebedev2007; Laurie et al. Reference Laurie, Boffetta, Falkovich, Kolokolov and Lebedev2014; Frishman Reference Frishman2017), with a well-defined mean radial vorticity profile. By considering the energy and angular momentum of the mean state, and considering fluctuations at leading order only, Laurie et al. (Reference Laurie, Boffetta, Falkovich, Kolokolov and Lebedev2014, hereafter L14) showed that the radial structure of the individual vortices in the dipole follows a simple power law outside a viscous core, taking the form
\begin{equation} {{\bar{\omega}}} =(r{\bar{u}})_r/r= U/r, \quad{\bar{u}}=U \end{equation}
in a polar coordinate system centred on each vortex. Here, 
$\omega$ is vorticity, 
$u$ is the azimuthal velocity, and an overbar denotes a suitable average (in the case of L14, a long-time and azimuthal average). In forced and frictionally dissipated flow, with equilibrated energy 
$ {\mathcal {E}}=\varepsilon /2\alpha$, the constant 
$U$ was predicted to satisfy 
$U^2=3\varepsilon /\alpha$, where 
$\varepsilon$ is the energy input rate, and 
$\alpha$ is the friction coefficient.
The form (1.1a,b) may be taken to hold at intermediate radial distances 
$ {r_c}< r< r_d$, outside a (hyper)viscous core of radius 
$ {r_c}$, and on scales smaller than the dipole separation 
$r_d$. Within the viscous core, 
$ {{\bar {\omega }}}$ may not increase without bound, and a simple solid rotation has been proposed, while at distances approaching 
$r_d$, the flow must assume a hyperbolic structure (Kolokolov & Lebedev Reference Kolokolov and Lebedev2016). Within the range of validity, the form (1.1a,b) and the numerical value of 
$U$ are both well supported by the numerical experiments of L14 and later ones at higher resolution (Frishman & Herbert Reference Frishman and Herbert2018).
In a different context, recent advances in the understanding of geophysical flows have established the ubiquity of inhomogeneous mixing on a background gradient of planetary vorticity. Under suitable conditions, the potential vorticity, the scalar product of absolute vorticity and stratification, can organize into an approximate staircase profile where mixing is confined to distinct latitudinal bands separated by near discontinuous jumps (Dritschel & McIntyre Reference Dritschel and McIntyre2008; Dunkerton & Scott Reference Dunkerton and Scott2008; Scott & Dritschel Reference Scott and Dritschel2012; Scott & Tissier Reference Scott and Tissier2012). The relation between potential vorticity and zonal velocity implies well-defined eastward jets at the jump locations. Similar inhomogeneous mixing has been identified in other contexts, notably vertical mixing on a background density gradient (Dritschel & McIntyre Reference Dritschel and McIntyre2008) and more recently in tokamak flows of neutral plasmas and other plasma flows (Dif-Pradalier et al. Reference Dif-Pradalier, Hornung, Garbet, Ghendrih, Grandgirard, Latu and Sarazin2017).
The development of such zonal mean flows can be considered in some sense as an instability of the background vorticity gradient, which is ‘an unstable equilibrium in the presence of mixing by waves and instabilities’ (Dunkerton & Scott Reference Dunkerton and Scott2008). The instability relies on the fact that Rossby waves propagating on eastward jets have critical layers on the jet flanks, where further mixing is favoured, sharpening the jet cores (Dritschel & McIntyre Reference Dritschel and McIntyre2008). A natural question then arises as to whether the mean vorticity gradient of (1.1a,b), embedded in a background turbulent flow, might also exhibit a susceptibility to inhomogeneous mixing, or similar self-organization of the mean vorticity. The numerical experiments of Dritschel, Qi & Marston (Reference Dritschel, Qi and Marston2015) on freely decaying turbulence on the sphere, in which late-time quadrupolar vortices were found to develop such step-like vorticity profiles, suggest that this might be the case.
In the context of the vorticity profile (1.1a,b), a vortex Rossby wave on a central core may likewise be expected to possess a critical layer at larger radii. By way of illustration, consider a condensate core in uniform solid body rotation such that 
$ {\bar {u}}=Ur/ {r_c}$ for 
$r< {r_c}$, continuously joined to the profile (1.1a,b) outside the core radius 
$ {r_c}$, as suggested by Kolokolov & Lebedev (Reference Kolokolov and Lebedev2016). The vorticity inside the core is 
$\omega =2U/ {r_c}$. Neglecting the vorticity gradient outside 
$ {r_c}$, waves on the core edge satisfy the dispersion relation
\begin{equation} \sigma={-}U(m-1)/{r_c}, \end{equation}
where 
$m$ is the azimuthal angular wavenumber. The angular phase speed of such waves, 
$\sigma /m$, would then match the angular speed 
$U/r$ of the condensate at radii 
$r=m {r_c}/(m-1)$, implying critical layers outside the core. This simplified argument neglects the effect of the vorticity gradient outside 
$ {r_c}$ on the frequencies of waves on the core. It also neglects all wave-like fluctuations on the mean vorticity profile outside 
$ {r_c}$, which may also be expected to possess their own critical layers. A complete description would require the calculation of quasi-modes on the full profile (1.1a,b) (e.g. Schecter et al. Reference Schecter, Dubin, Cass, Driscoll, Lansky and O'Neil2000; Turner Reference Turner2014). The point here is only that such critical layers may exist generally on such a mean flow, and that where they do, inhomogeneous mixing would be expected.
The development of a staircase profile in the geophysical context involves independent changes to the total energy and angular momentum of the system (Dunkerton & Scott Reference Dunkerton and Scott2008; Dunkerton Reference Dunkerton2018). As shown below, angular momentum preserving rearrangements of (1.1a,b) into a staircase profile involve an increase in total energy. Such changes are permitted in forced turbulent flows, in which total energy increases at a constant rate, or in which energy can be transferred between the mean and fluctuation fields. The extent to which these arguments are relevant to actual condensate flows, however, is limited; as will be shown in the numerical experiments below, the mean flow often develops a profile that is non-monotonic in radius. The process through which this occurs may be tentatively linked to the advection of high vorticity values from the inner edge of a critical layer region that does not have an analogue in the geophysical flows of the kind referred to here. The construction is nonetheless instructive, and the examples shown in § 2 are included more for the purposes of illustration than as quantitative predictions.
In this paper, we provide preliminary numerical evidence that spatially local self-organization of the condensate mean flow will occur in a forced system. Two sets of experiments are presented. In the first, an initially laminar vortex dipole, possessing the radial structure (1.1a,b), is subjected to turbulent fluctuations by a stochastic small-scale forcing. This allows the examination of turbulent mixing on a condensate-like vortex, without the need for the long spin-up times required to establish the condensate from rest. In the second set of experiments, the flow is started from rest and the condensate emerges out of the turbulent inverse energy cascade as in previous studies. Here, however, a time dependence is introduced into the energy input rate that allows for a more rapid establishment of the condensate and an attenuation of the turbulent fluctuations. The two sets of experiments may be considered as complementary, lying in some sense on either side of the traditional approach of forcing with weak, constant energy input over very long times.
The remainder of the paper is organized as follows. In § 2, we present some examples of idealized vorticity rearrangements, by analogy with the geophysical case, that illustrate the limiting forms that a vorticity staircase on the condensate may assume, and give an indication of the accompanying changes in angular momentum and energy that the profiles imply. In § 3, we outline the system of equations and implementation of the numerical experiments, recall under what conditions the condensate is expected, and summarize the configuration and parameter choices of the two sets of experiments. In § 4, we present the results of each set of experiments, illustrating in each case the tendency toward self-organization of the background condensate profile. Conclusions are given in § 5.
2. Staircase perturbations
From an ‘ideal’ condensate profile of the form given by (1.1a,b), we consider possible axisymmetric perturbations that may develop from a rearrangement of vorticity followed by azimuthal averaging. We restrict attention to a range of radii 
$r$ satisfying 
$ {r_c}< r< r_d$, and ignore details of the inner core and outer hyperbolic regions.
One of the simplest perturbations consists of an advective, downgradient rearrangement of vorticity over a finite interval 
$[ {r_0}, {r_1}]$, followed by an azimuthal average with uniform final vorticity on 
$[ {r_0}, {r_1}]$, as shown in figure 1(a):
\begin{equation} \omega=1/{r_m}, \quad r\in[{r_0},{r_1}]. \end{equation}
Here, 
$ {r_m}$ is the point at which the vorticity value of the original profile equals that of the mixed profile.
Figure 1. Illustrative staircase profiles on the ideal condensate: perturbed vorticity (red) and velocity (blue) with the unperturbed profile (1.1a,b) (dotted).
If the mixing leading to such a rearrangement is purely advective, then the area integrated vorticity is conserved. This ensures that the circulation at 
$r> {r_1}$, and hence the azimuthal velocity outside the mixing region, remains unchanged. In an infinite domain, this is the only kind of rearrangement permissible, without an infinite change in total energy. Direct integration of the vorticity difference over the annular area then yields the constraint
\begin{equation} 2{r_m}={r_0}+{r_1}. \end{equation} The blue curve in figure 1(a) shows the velocity profile corresponding to the mixed vorticity profile, calculated by direct integration of 
$r\omega = (ru)_r$. From this, it is immediately clear that the changes in both angular momentum and energy resulting from the perturbation must be negative. The actual change in angular momentum may be evaluated from
\begin{equation} \Delta M = \int_{{r_0}}^{{r_1}} ( m(r) - {\bar{m}}(r) ) r\,{\rm d}r, \end{equation}
where 
$m=ru$, i.e. 
$\omega =m_r/r$, and 
$ {\bar {m}}$ denotes the unperturbed profile. Integrating yields the simple result 
$\Delta M = -\delta ^3/12$, where 
$\delta = {r_1}- {r_0}$ is the step width.
A second simple case can be obtained from a non-advective upgradient rearrangement of vorticity, again across an interval 
$[ {r_0}, {r_1}]$, in which vorticity is increased on the interval 
$[ {r_0}, {r_m}]$ to its value at 
$ {r_0}$, and decreased on the interval 
$[ {r_m}, {r_1}]$ to its value at 
$ {r_1}$ (figure 1b):
\begin{equation} \omega= \begin{cases} 1/{r_0}, & r\in[{r_0},{r_m}], \\ 1/{r_1}, & r\in[{r_m},{r_1}]. \end{cases} \end{equation}
The vorticity is continuous at 
$ {r_0}$ and 
$ {r_1}$ but has a jump across 
$ {r_m}$. Such a profile cannot be obtained by advective mixing of the original vorticity, but must involve creation and destruction of mean vorticity through the action of eddy stresses 
$\overline {v'\omega '}$, primes denoting departures from the azimuthal mean. The circulation requirement that the mean flow outside 
$ {r_1}$ is unchanged now implies
\begin{equation} {r_m}^2={r_0}{r_1}. \end{equation} The jump in vorticity at 
$ {r_m}$ is associated with a local maximum in the mean velocity (blue curve), an analogue of the zonal jet profile of the potential vorticity staircase. From this mean velocity, it is clear that the changes in both energy and angular momentum are now positive. Such a profile could be achieved only if there is a transfer of both energy and angular momentum from the fluctuations to the mean flow. Integration of 
$m(r)- {\bar {m}}(r)$ separately on the intervals 
$[ {r_0}, {r_m}]$ and 
$[ {r_m}, {r_1}]$, and making use of (2.5), now gives 
$\Delta M = \delta ^3/24$.
Following Dunkerton & Scott (Reference Dunkerton and Scott2008), we can modify both the above perturbations in such a way that the angular momentum is unchanged from the original profile, 
$\Delta M=0$, through the addition of mixed zones on either side of the interval 
$[ {r_0}, {r_1}]$. In both cases, the vorticity is made uniform on intervals 
$[ {r_-}, {r_0}]$ and 
$[ {r_1}, {r_+}]$. In the first case, the vorticity values on these intervals are such that the vorticity is continuous at 
$ {r_-}$ and 
$ {r_+}$, increasing the jumps across 
$ {r_0}$ and 
$ {r_1}$, while in the second case, the values are such that the vorticity is continuous at 
$ {r_0}$ and 
$ {r_1}$, introducing new jumps at 
$ {r_-}$ and 
$ {r_+}$. See figures 1(c) and 1(d), respectively.
Again, the requirement that the change in circulation across the full interval 
$[ {r_-}, {r_+}]$ is zero provides one constraint on the two new parameters 
$ {r_-}$ and 
$ {r_+}$. The requirement that 
$\Delta M=0$ provides a second constraint, uniquely determining 
$ {r_-}$ and 
$ {r_+}$ for a given 
$ {r_0}$ and 
$ {r_1}$. The accompanying changes in total energy are of the same sign as the changes in the cases of figures 1(a) and 1(b), but an order of magnitude smaller.
The profiles shown here are purely for illustration of the type of mixing that might occur, in analogy with the generation of the potential vorticity staircase on the 
$\beta$-plane. They may be regarded as limiting cases, and the degree to which they are realized in practice will depend on the strength of the eddy fluctuations that accomplish the mixing. In numerical simulations at finite resolution, small-scale diffusion may additionally be expected to smooth out the jumps in vorticity between the mixed regions. Moreover, mixing that leads to non-monotonic profiles may be possible to a greater extent here than in the 
$\beta$-plane case, where a non-monotonic potential vorticity profile is known to be shear-unstable.
To illustrate a simple mixing process that might give rise to a profile resembling those just described, a simple experiment can be constructed in which an initial vorticity profile (1.1a,b) is subjected to an externally imposed straining flow.
To avoid the singular vorticity at 
$r=0$, we add a uniform core inside a radius 
$ {r_c}$ and take the initial vorticity to be
\begin{equation} {{\bar{\omega}}} = \begin{cases} 2U/{r_c}=2{{\varOmega_c}}, & r<{r_c},\\ U/r, & r>{r_c}. \end{cases} \end{equation}
The corresponding velocity is then 
$ {\bar {u}}=rU/ {r_c}$ for 
$r< {r_c}$, and 
$ {\bar {u}}=U$ for 
$r> {r_c}$.
The profile is forced by specifying an external straining flow 
$ {{\boldsymbol {u}}_{s}}$, with associated streamfunction 
$ {\psi _{s}}$, similar to the experiments of Turner (Reference Turner2014), so that the evolution of the vorticity is governed by the vorticity equation
\begin{equation} \omega_t + J(\psi + {\psi_{s}},\omega) =0, \quad\text{with } \psi=\nabla^{{-}2}\omega. \end{equation}The external streamfunction is defined to be
\begin{equation} {\psi_{s}}= \gamma(t)\,r^2 \cos(2\theta+\phi(t)), \end{equation}
where 
$\theta$ is the azimuthal angle, 
$\phi$ is a time-varying phase, and 
$\gamma$ is the amplitude. This straining flow is applied in the region of interest and modified at large 
$r$ to satisfy a no-normal-flow condition at the domain boundary. The domain boundary is set at large enough 
$r$ and has been verified to have negligible effect on the dynamics of interest, owing to the fact that vorticity gradients are very weak at large 
$r$.
The amplitude is grown gradually and then reduced again according to
\begin{equation} \gamma=\gamma_0 \exp({-(t-t_0)^2/\tau^2}). \end{equation}
The numerical values have been chosen as 
$t_0=250/ {{\varOmega _c}}$, 
$\tau =50/ {{\varOmega _c}}$, and peak amplitude 
$\gamma _0=0.8 {{\varOmega _c}}$, again for illustration only.
The phase 
$\phi (t)$ in (2.8) may be chosen to excite different modes of perturbation to the condensate profile. Taking a uniform rotation 
$\phi =\varOmega t$, the strain field can be made to project onto vortex Rossby waves that exist on the vorticity gradient of the initial vortex. Although the dispersion relation for these waves cannot be calculated analytically, specific values of 
$\varOmega$ can be expected to excite waves whose angular phase speeds match the local angular velocity of the vortex profile at some value of 
$r$. The resulting critical layer mixing will differ in details but can be expected to share similarities with that studied extensively in linear shear flows (e.g. Stewartson Reference Stewartson1978, among others).
Three examples are shown in figures 2(a,d), 2(b,e) and 2(c,f). As 
$\varOmega$ increases, the external strain matches the condensate angular velocity at decreasing radii. In all cases, the critical layer structure is clear in the early stage of the evolution (
$t=275/ {{\varOmega _c}}$, just past the peak forcing amplitude), with the advection of vorticity following closed streamlines in a finite range of radii. Beyond 
$t=250/ {{\varOmega _c}}$, the straining flow reduces in amplitude, and further forcing of the critical layer diminishes. The flow then symmetrizes azimuthally; by late time (
$t=1000/ {{\varOmega _c}}$), the wavenumber two structure has been mixed and vorticity homogenized across the original critical layer. Notably, there is a thin band of elevated vorticity values that persists at the outer edge of the original critical layer and appears to be dynamically stable, at least on the time scales considered here.
Figure 2. Snapshots of the vorticity field at (a–c) 
$t=275/\varOmega _c$ and (d–f) 
$t=1000/\varOmega _c$, at various 
$\varOmega$. Colour saturates white for 
$r< r_c$.
These examples are for illustration only, and to point out the analogy with mixing in geophysical flows. Details differ: for example, the thin band of elevated vorticity is not typically found in critical layer mixing on a linear background vorticity gradient, and seems due to the nonlinear variation of 
$\omega _r$ across the critical layer and, in particular, to the rapidly increasing 
$\omega$ on the inner edge. Similar elevated vorticity regions and non-monotonicity in radius turn out to be a fairly common feature of the fully turbulent simulations below.
3. Direct numerical simulation
3.1. Numerical model and configuration
We now examine the development of corresponding annular structure on a condensate, or condensate-like vorticity profile, in response to general small-scale stochastic forcing, typical of studies of forced-dissipative turbulence. Two sets of experiments are considered.
In the first set of experiments, an initial laminar vortex dipole is prescribed at 
$t=0$, and forcing and damping are then introduced to generate turbulent fluctuations while retaining the same total energy of the system as time evolves. This approach avoids the long spin-up stage needed to develop the condensate from an initial state of rest, and permits higher-resolution calculations and an exploration of a broader range of forcing parameters.
In the second set of experiments, the flow is started from rest and forcing applied that leads to the spontaneous generation of a condensate dipole in a background turbulent field. Crucially, the energy input rate used in this second set is time-varying, and allows for a rapid condensate development with diminishing fluctuation energy in time.
The two sets of experiments are complementary, and demonstrate the robustness of the tendency for annular self-organization to develop on the condensate profile.
In both cases, the equations solved are
\begin{equation} \omega_t + J(\psi,\omega) ={-}\alpha\omega + f + \nu\,\nabla^4\omega, \quad\text{with }\psi=\nabla^{{-}2}\omega, \end{equation}
where 
$\alpha$ is a coefficient of linear friction, and 
$f$ is a stochastic forcing, 
$\delta$-correlated in time, that inputs energy at rate 
$\varepsilon$ in a narrow band centred on wavenumber 
$ {k_f}$. A weak hyperdiffusion with coefficient 
$\nu$ is included for numerical stability. The equations are solved in a square periodic domain of length 
$L=2{\rm \pi}$. Simulations are performed at a resolutions from 
$N=512$ to 
$N=2048$ grid points in each direction using a standard pseudo-spectral method with de-aliasing via a spectral filter.
3.2. Initial conditions and parameters for simulation set A
In the first set of experiments (set A), the initial condition is given by
\begin{equation} \omega_0=\sum_{i=0}^{1} \frac{({-}1)^i U}{|{{\boldsymbol{x}}}-{{\boldsymbol{x}}}_i|} + \text{periodic extension}, \end{equation}
where 
$ {{\boldsymbol {x}}}_0 = ({\rm \pi} /2,{\rm \pi} /2)$ and 
$ {{\boldsymbol {x}}}_1 = (-{\rm \pi} /2,-{\rm \pi} /2)$ in a domain of size 
$L=2{\rm \pi}$.
We note that the root mean square (r.m.s.) velocity of the flow defined by (3.2) satisfies
\begin{equation} U^2\approx2.2 {U_{rms}}^2, \end{equation}
which is comparable but distinct from the theoretical prediction 
$U^2=3 {U_{rms}}^2$ put forward in L14. The value 
$2.2$ was obtained by evaluating 
$ {U_{rms}}$ numerically for a given 
$U$, and was found to be insensitive to the details in the condensate core. The difference between this value and the theoretical prediction may indicate that the profile (1.1a,b) does not hold at large radii, comparable to the dipole separation (i.e. in the hyperbolic region). The energy input rate and frictional damping rate are chosen to satisfy 
$\alpha =\varepsilon / {U_{rms}}^2$ so as to ensure that the total energy remains constant in time.
In L14, the establishment of the condensate is controlled by the parameter
\begin{equation} \delta= \alpha (L^2/\varepsilon)^{1/3}, \end{equation}
the ratio of frictional time scale to a large-scale eddy time scale. The condensate emerges for 
$\delta \ll 1$, that is, when the energy of the condensate significantly dominates the energy of the turbulent fluctuations. Values of 
$\delta$ down to 
$10^{-3}$ were considered in Frishman & Herbert (Reference Frishman and Herbert2018). Here, we consider energy input rates 
$\varepsilon =U^3 \beta /L$, for 
$\beta$ a dimensionless amplitude with 
$\beta =2^{ {n_f}} 2{\rm \pi} \times 10^{-7}$, where 
$ {n_f}=-2,\ldots,5$ is an integral forcing parameter. The values of 
$ {n_f}$ used here give 
$\delta$ in the range
\begin{equation} \delta=(\varepsilon L)^{2/3}/{U_{rms}}^2 \in [0.064,2.6]\times 10^{{-}3}, \end{equation}
where we have used 
$\alpha =\varepsilon / {U_{rms}}^2$.
A related parameter can be constructed by considering the ratio of an enstrophy based time scale 
$\eta ^{-1/3}$, where 
$\eta = {k_f}^2\varepsilon$ is the enstrophy input rate, with a time scale constructed from the mean flow and the forcing scale, 
$1/ {k_f} U$. We thus define
\begin{equation} \zeta=\frac{\eta^{1/3}}{{k_f} U} = \frac{{U_{rms}}^2}{U^2}\,\frac{\delta^2}{L{k_f}}, \end{equation}
and from its definition, we expect that 
$\zeta$ might influence the level of small-scale mixing of the mean flow by the forced vorticity fluctuations. This interpretation is supported by considering a rescaling of the length scales of the problem. Rescaling length by a factor 
$\mu$, so that 
$r\rightarrow \mu r$, corresponds to a rescaling of the forcing scale as 
$ {k_f}\rightarrow \mu ^{-1} {k_f}$. The radial vorticity profile then scales as 
$\omega =U/r\rightarrow \mu ^{-1}\omega$, implying that the relevant vortex turnaround time scales as 
$\tau \sim \omega ^{-1}\rightarrow \mu \tau$. The scaling of the turnaround time in turn suggests a scaling of both friction and energy input rate as 
$\alpha \rightarrow \mu ^{-1}\alpha$ and 
$\varepsilon \rightarrow \mu ^{-1}\varepsilon$, so that 
$ {U_{rms}}$ is fixed also, as is usual in the case of constant forcing with friction. Since 
$\zeta$ can also be written as
\begin{equation} \zeta=\left(\frac{\varepsilon}{{k_f} U^3}\right)^{1/3}, \end{equation}
it follows that 
$\zeta$ is invariant under the scaling by 
$\mu$. In this rescaling, forcing at larger 
$ {k_f}$ essentially leads to mixing at smaller distances 
$r$ from the vortex core, and on a shorter time scale. This prediction is supported by the numerical results reported below.
Simulations have been carried out with three different values of 
$ {k_f}$: 
$ {k_f}=16$, 32, 64. Following the above scaling, the grid resolution is kept proportional to 
$ {k_f}$, with 
$N=512$, 1024, 2048 grid points in each direction, respectively. Similarly, the integration time is set proportional to 
$ {k_f}^{-1}$. Defining a domain-scale turnaround time based on the initial condition, 
$\tau _0=L/U$, we integrate to 
$t=T$, where 
$T=1600\tau _0$, 
$800\tau _0$, 
$400\tau _0$, for 
$ {k_f}=16$, 32, 64. The full set of parameter values is summarized in table 1
Table 1. Parameter values used in simulation set A.
3.3. Forcing time dependence and parameters for simulation set B
In the second set of simulations (set B), the vorticity field is initially zero and subjected to the stochastic forcing 
$f$. Here, we induce the rapid development of the condensate by using a time-varying energy input rate 
$\varepsilon (t)$, which is initially large, to force strong turbulent fluctuations, but which then decreases as the condensate emerges, to allow the condensate mean flow to dominate the turbulent fluctuations in a shorter time than would be achievable with constant 
$\epsilon$. Further, the time dependence of 
$\varepsilon (t)$ is chosen to mimic the damping of constant energy input by frictional dissipation.
When 
$\varepsilon =\varepsilon _0$ is constant, frictional dissipation at rate 
$\alpha$ gives a total energy that grows in time as
\begin{equation} {\mathcal{E}}(t)=\frac{\varepsilon_0}{2\alpha}\,( 1 - {\rm e}^{{-}2\alpha t} ), \end{equation}
tending towards the steady state value 
$\varepsilon _0/2\alpha$. To mimic this energy growth with time-varying forcing, we consider a frictionless system with 
$\alpha =0$. Energy then grows simply according to
\begin{equation} \dot{{\mathcal{E}}} = \varepsilon(t). \end{equation}By choosing
\begin{equation} \varepsilon(t)=\varepsilon_0\,{\rm e}^{{-}t/{\tau_f}}, \end{equation}
we recover the energy growth of (3.8) provided that 
$ {\tau _f}^{-1}$ is identified with the 
$2\alpha$ of the constant 
$\varepsilon =\varepsilon _0$ case. The total energy in the limit 
$t\rightarrow \infty$ is 
$\varepsilon_0 {\tau _f}$. One advantage of this time dependence is that an initially rapid energy input allows turbulent fluctuations, through the inverse cascade, to reach the domain scale relatively quickly, generating a condensate. At later times, as the condensate develops, the turbulent energy fluctuations decrease and the energy in the condensate dominates.
In the simulations presented below, we have chosen values for 
$\varepsilon _0$ and 
$ {\tau _f}$ that give a similar 
$ {U_{rms}}$ to that of the simulations in set A. We thus take a range 
$\varepsilon _0=2,4,8,\ldots$, with corresponding 
$ {\tau _f}=320,160,80,\ldots$. All simulations are run to time 
$T=4 {\tau _f}$, which gives a final energy input rate of approximately 
$\varepsilon _T=0.0183\varepsilon _0$.
Following the scaling argument above, we expect more rapid mixing at larger 
$ {k_f}$, and therefore take generally smaller values of 
$ {\tau _f}$ and larger values of 
$\varepsilon _0$ at larger 
$ {k_f}$. Finally, the equivalent parameters 
$\delta$ and 
$\zeta$ can be considered based on the final energy input rate 
$\varepsilon _T$ and the actual 
$ {U_{rms}}$ to give an approximate characterization of the forcing at the end of each simulation. Values of the various parameter choices for this set of simulations are summarized in table 2.
Table 2. Parameter values used in simulation set B.
4. Results
4.1. Forcing of an initially laminar vortex dipole
Snapshots of the vorticity evolution at the final time, 
$t=T$, are shown in figure 3 for a variety of forcing wavenumbers 
$ {k_f}=16,32,64$ and forcing strengths (weak–intermediate forcing in figures 3(a–c), stronger forcing in figures 3(d–f)). In all cases, a clear annular banding of the initially smooth vorticity gradient has occurred by the final time, though with considerable variation both across different cases and even between the positive and negative vortices of a single case. The banding may take the form of a uniform step-like but broadly monotonic distribution (e.g. the negative (blue) vortices in figures 3a,e), or a more distinct non-monotonic distribution (e.g. the positive (orange) vortex in figure 3(b) or the negative vortex in figure 3(d). The radial scale of the banding decreases with increasing 
$ {k_f}$ and is harder to make out in figures 3(c,f), 
$ {k_f}=64$, which we examine in more detail below. The vortices remain broadly axisymmetric, with stronger fluctuations visible in the more strongly forced cases, as expected.
Figure 3. Snapshots of the full vorticity field at 
$t/T=1$ for (a,d) 
$ {k_f}=16$, (b,e) 
$ {k_f}=32$, and (c,f) 
$ {k_f}=64$, and (a–c) forcing index 
$ {n_f}=0,1,2$ (intermediate forcing) and (d–f) forcing index 
$ {n_f}=3,4,5$ (strong forcing). Orange denotes positive vorticity, blue denotes negative vorticity and the colour scale saturates to white at value 
$|\omega |=U/ {r_c}$ with 
$ {r_c}=0.1$, i.e. at a distance 
$ {r_c}$ from the vortex centre in the case of the exact condensate profile.
The fluctuation field can be examined more clearly by removing the mean flow. Here, the mean flow is defined instantaneously as an azimuthal average centred on each vortex, thus retaining the banded radial structure observed in the full fields. The fluctuation vorticity, defined as the full vorticity minus this mean flow, is shown in figure 4 for the same cases as in figure 3. The colour saturates to white at vorticity values proportional to 
$ {k_f}$, and in all cases much lower values than for the full fields in figure 3. Across each level of forcing (intermediate in figures 4a–c, strong in figures 4d–f), the intensity of the fluctuation field appears approximately uniform in the scaled colour range, supporting the scaling discussed § 3.2, namely that the flow is approximately invariant under 
$( {k_f},\varepsilon,\omega )\rightarrow \mu ^{-1} ( {k_f},\varepsilon,\omega )$.
Figure 4. As figure 3, but for the fluctuation vorticity field. The colour saturates to white at vorticity values 
$|\omega |=U/2$ for (a,d) 
$ {k_f}=16$, 
$|\omega |=U$ for (b,e) 
$ {k_f}=32$ and 
$|\omega |=2U$ for (c,f) 
$ {k_f}=64$.
It is important to note that the r.m.s. vorticity fluctuations equilibrate at very early times in each calculation. Fields plotted at times as early as 
$T/10$ (not shown) are qualitatively indistinguishable from those of figure 4. Energy spectra are shown in figure 5 for the strong forcing cases, with lines plotted for times 
$t=T/2$ and 
$t=T$ being almost indistinguishable.
Figure 5. Energy spectra for (a) 
$ {k_f}=16$, 
$n_f=3$, (b) 
$ {k_f}=32$, 
$n_f=4$ and (c) 
$ {k_f}=64$, 
$n_f=5$, plotted at 
$t=T/2$ (solid) and 
$t=T$ (dotted). Straight lines have slope 
$-3$ and are for reference only.
Figure 6 shows profiles of the azimuthally averaged vorticity field 
$ {{\bar {\omega }}}$ centred on each vortex for the full range of simulations, together with an enlargement of each vortex region for the selected cases shown by the blue lines in figures 6(a–c). Organization of the initial vorticity profile (dotted) is clearly visible in multiple bands surrounding each vortex. In some places, the turbulent fluctuations have flattened the vorticity over an annular region, while in other places, the azimuthal mean vorticity profile has developed a significant non-monotonicity in the radial direction. The processes giving rise to the annular structures are highly variable, but often can be clearly linked to cat's-eye-like structures, similar to those found in critical layers, as in figure 6(e), for example. The scale of the annular regions appears to scale approximately with 
$ {k_f}^{-1}$, in both the profiles and the detailed vorticity field – noting that the range of 
$r$ in figures 6(a–c), and of 
$(x,y)$ in figures 6(d–f), is also rescaled with 
$ {k_f}^{-1}$ – consistent with the scaling arguments given in § 3.2. Indeed, there is no other preferred length scale in the system, other than the core radius 
$r_c$, which is smaller and does not appear to have a significant effect. (Here, 
$r_c$ is set by the numerical truncation, which is also scaled with 
$ {k_f}^{-1}$ and so fixed relative to the plot range. While it is possible that the scale of the annular regions could be controlled by 
$r_c$, it seems much more likely that it is controlled by the forcing scale and the induced turbulent fluctuations.)
Figure 6. Azimuthally averaged vorticity 
$ {{\bar {\omega }}}(r)$ for (a) 
$ {k_f}=16$, (b) 
$ {k_f}=32$ and (c) 
$ {k_f}=64$, and a range of forcing strengths 
$n_f$, with lines corresponding to different 
$n_f$ offset in the vertical for clarity; dotted lines correspond to the initial 
$ {{\bar {\omega }}}(r)$. Detailed view of the full vorticity field for the cases in (a–c) indicated by the blue lines: (d) 
$ {k_f}=16$, 
$n_f=2$, (e) 
$ {k_f}=32$, 
$n_f=3$ and (f) 
$ {k_f}=64$, 
$n_f=4$. Magenta lines in (a–c) correspond to the cases shown in figure 3.
Self-organization into annular structures is less distinct at weaker forcing (bottom two curves in each of figures 6a–c), but is still present: there is clear flattening of the initial vorticity profile over some regions, e.g. around 
$r=0.8$–
$1.6$ in the case 
$ {n_f}=-2$ in figure 6(a), or on either side of 
$r=0.25$ in the case 
$ {n_f}=0$ in figure 6(c). However, large overshoots are mainly found for stronger forcing, suggesting that their generation requires a certain threshold of the fluctuation vorticity field.
The annular distributions shown are typical of other times in the simulations. To illustrate the temporal evolution, we construct the relative vorticity difference
\begin{equation} \Delta {{\bar{\omega}}}(r,t) = \frac{{{\bar{\omega}}}(r,t)-{{\bar{\omega}}}(r,0)}{{{\bar{\omega}}}(r,0)}. \end{equation}
This is shown in figure 7 for the same cases as in figure 3, over the full time of the integrations from 
$t=0$ to 
$t=T$. The radial range is again proportional to 
$ {k_f}^{-1}$. Alternating bands of positive and negative azimuthal vorticity develop in all cases from early time, and persist with varying degrees of transience throughout the duration of the simulation. In some cases, the bands are remarkably steady over the time interval considered (e.g. figure 7(b), top). In general, there appears to be greater transience in the more strongly forced cases (figures 7d–f), which might be expected on account of stronger turbulent fluctuations. In these cases, the relative azimuthal vorticity differences are stronger than for weaker forcing, consistent with the stronger non-monotonicity in the upper versus lower magenta lines in figure 6. The degree of variance that might be expected in an ensemble can be inferred to some extent from the differences between the positive (upper) and negative (lower) vortices of each case. Various of the lower-resolution calculations have been run beyond 
$t=T$ (not shown) and indicate that the patterns found here persist on longer time scales, with occasional regime-like behaviour with more or less well-defined annular structures persisting over prolonged periods.
Figure 7. Relative vorticity difference 
$\Delta {{\bar {\omega }}}(r,t)$ for (a–c) 
$ {k_f}=16,32,64$ and 
$ {n_f}=0,1,2$, respectively, and (d–f) 
$ {k_f}=16,32,64$ and 
$ {n_f}=3,4,5$, respectively. Orange and blue denote positive and negative differences, with colour saturating white for 
$\Delta {{\bar {\omega }}}=0.2$.
To quantify the dependence of the annular structure on forcing strength and strength of the turbulent fluctuations, we calculate the total, radially integrated annular vorticity difference 
$Z$, as well as the r.m.s. of the fluctuation vorticity field (that shown in figure 4) 
$F$:
\begin{equation} Z = \int_0^{r_0} |\Delta{{\bar{\omega}}}| \,{\rm d}r \quad\text{and}\quad F = \int_\mathcal{D} |\omega'| \,{\rm d}\kern 0.06em x\,{\rm d}y. \end{equation}
In 
$Z$, we have omitted the area weighting, since 
$\Delta \omega$ itself contains a normalization by 
$\omega _0\propto 1/r$; 
$Z$ therefore gives an area-weighted measure of the azimuthal vorticity difference over the range of integration. Here, 
$r_0$ has been taken to be 
$0.7\times 32/ {k_f}$, i.e. 0.7 of the domain shown in figure 7.
The measures are summarized in figure 8 for all values of 
$ {k_f}$ and 
$ {n_f}$ considered, time-averaged over the last half of each simulation. There is considerable variance in 
$Z$, but a general approximately linear dependence on 
$\zeta$ holds across the range of forcing strengths considered. For given 
$\zeta$, there is little systematic variation with 
$ {k_f}$.
Figure 8. (a) Integrated azimuthal vorticity difference 
$Z(\zeta )$ and (b) integrated fluctuation vorticity 
$F(\zeta )$ (log scales) for different energy input rates and forcing scales. Values have been time-averaged over the last half of each simulation. Values of 
$2^{ {n_f}}/ {k_f}$ range from 
$1/64$ (leftmost points) to 
$1/2$ (rightmost points), corresponding to 
$\zeta \propto (2^{n_f}/k_f)^{1/3}$ between 0.0016 and 0.00368 (see table 1).
The variation of the fluctuation field has less variance and shows bigger departures from a simple scaling with 
$\zeta$. The generally weaker 
$F$ for 
$ {k_f}=16$ might plausibly be attributed to insufficient scale separation between the forcing scale and the vortex scale, and the fact that the simple scaling, which neglects many aspects of the problem, is not satisfied. On the other hand, there is good agreement between the 
$ {k_f}=32$ and 
$ {k_f}=64$ cases. At these these forcing values, there is a clear monotonic increase in fluctuation strength with 
$\zeta$, though it is not a linear or even power-law dependence over the range considered.
4.2. Annular self-organization in a spontaneously developing condensate
The second set of simulations approaches the final condensate from the opposite direction. Instead of starting with a pure laminar condensate profile and adding turbulent fluctuations, the vortex dipole condenses self-consistently from the turbulent inverse energy cascade. The process is sped up by allowing for the time-dependent energy input rate described in § 3.3, in which forcing of the turbulent fluctuations decays exponentially in the process.
The energy, enstrophy and vorticity maxima for the simulations set summarized in table 2 are shown in figure 9. In all cases, the energy grows in close agreement with (3.8). The enstrophy (red) shows an early-time dependence on 
$ {\tau _f}$, since shorter 
$ {\tau _f}$ corresponds to higher 
$\varepsilon _0$ and stronger initial vorticity fluctuations. The variation with 
$ {\tau _f}$ decreases with time as fluctuations die out and the condensate regime is established. The vorticity maxima are related more closely to the condensate core, rather than the fluctuations, and consequently show no systematic dependence on 
$ {\tau _f}$ but a clear scaling with resolution (proportional to 
$ {k_f}$), which controls the size of the condensate core. In all cases, an approximate equilibrium is being approached towards the end of the simulation time 
$t=4 {\tau _f}$.
Figure 9. Normalized energy (black), enstrophy (red) and vorticity maximum (blue), where 
$\mathcal {E}_0=\varepsilon _0\tau _f$ for (a) 
$ {k_f}=16$ and (b) 
$ {k_f}=32$. Solid, dotted and dashed lines indicate 
$ {\tau _f}=40,80,160$ in (a), and 
$ {\tau _f}=20,40,80$ in (b), respectively.
Snapshots at times 
$t/ {\tau _f}=0.25,1,4$ of the full vorticity field are shown in figure 10 for the cases 
$ {k_f}=32$ and 
$ {\tau _f}=20$ (figures 10a–c) and 
$ {\tau _f}=80$ (figures 10d–f). They illustrate the emergence of the condensate from an initially turbulent flow comprising a distribution of vortices across area; see Burgess & Scott (Reference Burgess and Scott2017) for a description of this forced vortex population at earlier times, before the energy reaches the domain scale. Here, at 
$t/ {\tau _f}=0.25$, energy at the domain scale is becoming significant, with the emerging condensate dipole just visible in the case 
$ {\tau _f}=80$. The condensate is more clearly defined by 
$t/ {\tau _f}=1$, at which time there is still considerable energy in the turbulent fluctuations. By 
$t/ {\tau _f}=4$, the energy in turbulent fluctuations is significantly reduced. Note that at both 
$t/ {\tau _f}=1$ and 
$t/ {\tau _f}=4$, a distinct banding is visible on the condensate dipole, more pronounced at the later time.
Figure 10. Snapshots of the full vorticity field for 
$ {k_f}=32$ at (a) 
$t=0.25 {\tau _f}$, 
$ {\tau _f}=20$, (b) 
$t= {\tau _f}$, 
$ {\tau _f}=20$, (c) 
$t=4 {\tau _f}$, 
$ {\tau _f}=20$, (d) 
$t=0.25 {\tau _f}$, 
$ {\tau _f}=80$, (e) 
$t= {\tau _f}$, 
$ {\tau _f}=80$ and (f) 
$t=4 {\tau _f}$, 
$ {\tau _f}=80$.
As in the previous subsection, the turbulent fluctuations can be separated by defining an azimuthal mean profile centred on each vortex at each time, and subtracting this mean profile from the full vorticity field. The fluctuation field at 
$t/ {\tau _f}=4$ is shown in figure 11 for the cases 
$ {k_f}=32$ and 
$ {\tau _f}=20,40,80$; figures 11(a,c) correspond to figures 10(c,f). The colour scale is the same across all cases, saturating to white at a vorticity value one-sixth that of figure 10, and shows a stronger fluctuation field at smaller 
$ {\tau _f}$ (larger 
$\varepsilon _0$). Also just discernible in figures 11(a–c) is a qualitative difference in the structure of the fluctuations in a square-lattice-like domain separating the individual dipole vortices. This difference was even more apparent in the experiments of set A shown in figure 4. Their existence in the set B experiments confirms that they are not artefacts of the initial condition but emerge as inherent parts of the condensate flow. Presumably, they are related to the hyperbolic regions discussed in Kolokolov & Lebedev (Reference Kolokolov and Lebedev2016), but they appear to exist in relatively narrow regions of the domain along the separatrices of the dipole.
Figure 11. Snapshots of the fluctuation vorticity field at 
$t=4 {\tau _f}$ for 
$ {k_f}=32$ and (a) 
$ {\tau _f}=20$, (b) 
$ {\tau _f}=40$ and (c) 
$ {\tau _f}=80$. The colour saturates to white at vorticity values one-sixth of those in figures 10(c,f).
Energy spectra for the same cases 
$ {k_f}=32$ and 
$ {\tau _f}=20,40,80$ are shown in figure 12, with lines indicating times 
$t/ {\tau _f}=0.25,1,4$. They confirm the above picture, namely that at 
$t/ {\tau _f}=0.25$, energy has reached the domain scale but the 
$k^{-3}$ profile of the condensate has not fully developed. At scales smaller than the forcing scale, there is a clear decrease in the energy level with increasing 
$ {\tau _f}$, and in all cases a distinct break between the large- and small-scale distributions.
Figure 12. Energy spectra for (a) 
$ {\tau _f}=20$, (b) 
$ {\tau _f}=40$ and (c) 
$ {\tau _f}=80$, plotted at 
$t= {\tau _f}/4$ (solid), 
$t= {\tau _f}$ (dotted) and 
$t=4 {\tau _f}$ (dashed). Straight lines have slope 
$-3$ and are for reference only.
Finally, figure 13 shows the azimuthal mean vorticity for each pair of the condensate dipole, and for all simulations in this set. Individual lines indicate the instantaneous 
$ {{\bar {\omega }}}(r)$ equally spaced between 
$t/ {\tau _f}=2$ and 
$t/ {\tau _f}=4$, offset in the vertical for clarity. There is a significant degree of variation from one vortex to another, although the basic 
$1/r$ dependence is clear in all cases. On top of this general 
$1/r$ dependence, however, discrete regions of annular banded fluctuations again appear spontaneously, mostly persisting in time over the duration of the simulation, though with some transient behaviour (e.g. figure 13(d) right, figure 13(f) left). Comparing figures 13(a–c) and 13(d–f), the width of the annular regions appears to scale approximately with 
$ {k_f}^{-1}$, though again with considerable variation across individual cases. There is no systematic indication that the annular structure weakens as the time scale 
$ {\tau _f}$ is increased.
Figure 13. Azimuthally averaged vorticity 
$ {{\bar {\omega }}}(r)$ for (a–c) 
$ {k_f}=16$, 
$ {\tau _f}=40,80,160$, respectively, and (d–f) 
$ {k_f}=32$, 
$ {\tau _f}=20,40,80$, respectively. Individual lines represent equally spaced times from 
$t=2 {\tau _f}$ to 
$t=4 {\tau _f}$, offset in the vertical for clarity (time increasing upwards). The range of 
$r$ values is proportional to 
$ {k_f}^{-1}$. The left and right plots in each panel denote the two vortices of each dipole.
5. Conclusions
While the condensate profile (1.1a,b) put forward by L14 may be obtained under very-long-time averaging, a distinct radial structure appears over shorter time scales, in which the background vorticity gradient is reduced over well-defined annular intervals. The structure appears to be generated by a self-organization of the background vorticity gradient of the long-time average profile by turbulent fluctuations, in a manner analogous to the inhomogeneous mixing of the background potential vorticity gradient on a rotating planet: rearrangement of the mean vorticity by the turbulent fluctuations gives rise to annular regions over which the vorticity gradient is reduced, and can even change sign, giving rise to a non-monotonic profile. The appearance of the annular regions is robustly observed both in simulations in which the flow begins in a pre-existing, ideal, laminar condensate profile before being subject to mixing by turbulent fluctuations, and in simulations in which the condensate develops spontaneously from the forced inverse energy cascade. Notably, there was little discernible difference between the two sets of simulations, in terms of the fluctuation vorticity field, energy spectra and the development of the annular zones.
Although they persist on the time scales of the current simulations, transience in the annular structures on longer time scales may obscure their observation in very-long-time averages. Probably for this reason, they were not detected in the numerical simulations of L14 and the subsequent studies (Frishman Reference Frishman2017; Frishman & Herbert Reference Frishman and Herbert2018), particularly at the larger values of 
$\delta$ considered in those works (for which transience is greater). It is speculated that the ‘zero modes’ of Frishman & Herbert (Reference Frishman and Herbert2018) may be related to the annular structures found here, and that the vorticity staircase construction might provide a useful framework for their further analysis.
In 
$\beta$-plane turbulence, the emergence of the potential vorticity staircase depends on the dimensionless parameter 
$ {L_{Rh}}/ {L_{\varepsilon }}$, where 
$ {L_{Rh}}=\sqrt {U/\beta }$ is the Rhines scale, with characteristic velocity 
$U$, and 
$ {L_{\varepsilon }}=(\varepsilon /\beta ^3)^{1/5}$. Strong staircases are found for 
$ {L_{Rh}}/ {L_{\varepsilon }}\gtrsim 6$ (Scott & Dritschel Reference Scott and Dritschel2012). Corresponding local length scales for the condensate can be defined analogously by
\begin{equation} {L_{Rh}}=\sqrt{U/{{\bar{\omega}}}_r}=r \quad \text{and} \quad {L_{\varepsilon}}=(\varepsilon/{{\bar{\omega}}}_r^3)^{1/5}, \end{equation}which results in
\begin{equation} {L_{Rh}}/{L_{\varepsilon}} = (U^3/\varepsilon r)^{1/5}. \end{equation}
At a given radius, this is proportional to 
$\delta ^{-3/10}$. At 
$r<1$, values 
$\delta \sim 10^{-3}$ give 
$ {L_{Rh}}/ {L_{\varepsilon }}\gtrsim 10$, which is well within the strong staircase regime, where inhomogeneous mixing of the background vorticity might be expected.
The annular width of the annular regions observed scales approximately with the forcing scale, suggesting that they arise in response to turbulent fluctuations at the scale of the forcing itself, and therefore that the inverse energy cascade plays a limited role in organizing the fluctuation vorticity before it mixes the background gradient. This means that at very small forcing scale, the scale of the annular regions will become increasingly small. However, a small forcing scale requires a correspondingly higher numerical resolution and weaker (hyper)diffusion (Kolokolov & Lebedev Reference Kolokolov and Lebedev2016). Consequently the (hyper)viscous core of the condensate vortex will be correspondingly smaller, and the vorticity values correspondingly higher. Thus a self-similarity exists in which higher 
$ {k_f}$ corresponds to a smaller and more intense vortex core, and narrower annular zones at smaller radii, developing over correspondingly shorter time scales. This self-similarity is encapsulated by the observation of similar dynamics for similar values of the parameter 
$\zeta$ defined by (3.6). However, because the annular zones are generated on a scale proportional to 
$k_f^{-1}$, they will become manifest only once the condensate dipole has reached a certain intensity and scale.
Finally, the results of the simulation set A suggest that the degree of self-organization scales approximately with the parameter 
$\zeta$, as shown in figure 8(a). This suggests that for sufficiently weak forcing (small 
$\zeta$), the annular zones may eventually cease to be observed, and the condensate may approach more closely to its theoretical form. To attain such a regime in a flow that is forced from a state of rest and in which the condensate develops spontaneously would, however, require weak forcing over extremely long time scales, far beyond what is currently computationally feasible. For now, therefore, it remains an open question whether a perfect universal condensate in the sense of (1.1a,b) might eventually be attained with small enough energy input rates.
Declaration of interests
The author reports no conflict of interest.
Open access
