Quasi-two-dimensional nonlinear evolution of helical magnetorotational instability in a magnetized Taylor-Couette flow

Magnetorotational instability (MRI) is one of the fundamental processes in astrophysics, driving angular momentum transport and mass accretion in a wide variety of cosmic objects. Despite much theoretical/numerical and experimental efforts over the last decades, its saturation mechanism and amplitude, which sets the angular momentum transport rate, remains not well understood, especially in the limit of high resistivity, or small magnetic Prandtl numbers typical to interiors (dead zones) of protoplanetary disks, liquid cores of planets and liquid metals in laboratory. Using direct numerical simulations, in this paper we investigate the nonlinear development and saturation properties of the helical magnetorotational instability (HMRI) -- a relative of the standard MRI -- in a magnetized Taylor-Couette flow at very low magnetic Prandtl number (correspondingly at low magnetic Reynolds number) relevant to liquid metals. For simplicity, the ratio of azimuthal field to axial field is kept fixed. From the linear theory of HMRI, it is known that the Elsasser number, or interaction parameter determines its growth rate and plays a special role in the dynamics. We show that this parameter is also important in the nonlinear problem. By increasing its value, a sudden transition from weakly nonlinear, where the system is slightly above the linear stability threshold, to strongly nonlinear, or turbulent regime occurs. We calculate the azimuthal and axial energy spectra corresponding to these two regimes and show that they differ qualitatively. Remarkably, the nonlinear state remains in all cases nearly axisymmetric suggesting that HMRI turbulence is quasi two-dimensional in nature. Although the contribution of non-axisymmetric modes increases moderately with the Elsasser number, their total energy remains much smaller than that of the axisymmetric ones.


Introduction
The magnetorotational instability (MRI, [1,2]) is one of the most important processes in magnetized differentially rotating conducting fluids, which largely determines their dynamics and evolution. It is a powerful linear instability arising as a result of the combined effect of weak magnetic field and radially decreasing angular velocity. The MRI is believed to operate in a vast variety of cosmic objects, ranging from astrophysical disks and stars to liquid-metal cores of planets. While discovered as early as 1959 [1], the astrophysical significance of MRI was first recognized only three decades later by Balbus & Hawley [3], who demonstrated that weak magnetic fields can destabilize Keplerian accretion disks around such diverse objects as supermassive black holes, black holes in X-ray binaries and young stellar objects, which would otherwise be linearly stable according to Rayleigh's criterion. The linear growth of MRI, which taps into the free energy of differential rotation, eventually breaks down into magnetohydrodynamic (MHD) turbulence [4,5]. This turbulence transports angular momentum outward and, as a consequence, matter inward in the disk, yielding mass accretion rates onto the central object close to observationally inferred values. More recently, MRI has also been discussed for explaining angular momentum transport in the Sun, massive stars and neutron stars [6,7,8,9] and also in the context of the geodynamo [10,11].
The MRI was originally discovered theoretically in a classical Taylor-Couette (TC) flow of a conducting fluid between two concentric rotating cylinders threaded by an external vertical (along the common axis of cylinders) magnetic field [1,12] -a setup which is also best suited and most frequently used nowadays for experimental investigations. The conducting medium filling the gap between cylinders is usually a liquid metal (sodium, gallium) characterized by an extremely small ratio of viscosity to Ohmic resistivity, or magnetic Prandtl number, P m = ν/η ∼ 10 −6 − 10 −5 . By suitably adjusting the rotation rates of outer and inner cylinders, the TC flow profile can be made very close to the Keplerian one [13,14,15,16], offering a unique possibility to study the disk MRI problem in laboratory as well, which has been up to now mostly carried out both via analytical means and numerical simulations.
First experimental efforts to study MRI in the laboratory were made at the University of Maryland [17], at Princeton University [18,19], and at Helmholtz-Zentrum Dresden-Rossendorf (HZDR) [20,21,22]. The liquid sodium spherical Couette experiment in Maryland had produced coherent velocity/magnetic field fluctuations showing up in a parameter region reminiscent of MRI [17]. The liquid GaInSn TC experiments in Princeton was designed to investigate the standard version of MRI (SMRI) with only an axial magnetic field being imposed. In this case, the azimuthal magnetic field (which is an essential participant in the MRI process) must be produced from the axial field by induction effects, which are proportional to the magnetic Reynolds number Rm of the flow. Rm, in turn, is proportional to the hydrodynamic (HD) Reynolds number Re according to Rm = P m · Re. Therefore, in order to achieve Rm 1 necessary for SMRI to operate (see e.g., [23,24]), taking into account the very small values of P m, Re ∼ 10 5 − 10 6 is needed. This makes SMRI experiments with TC flows extremely challenging, although evidence for slow magneto-Coriolis waves [18] and for a free-Shercliff layer instability [19] have already been obtained. Within the DRESDYN project at HZDR [25], it is planned to set-up a large liquid sodium TC experiment, which will allow to broaden the parameter range in such a way as to make SMRI reachable.
On the theoretical side, SMRI in TC flow has been extensively studied both in the linear and nonlinear regimes, with more focus on the low-P m regime as it is typical in experiments. Still, the typical Reynolds numbers Re ∼ 10 3 used in these analyses are orders of magnitude smaller than experimental values Re 10 6 , which are extremely demanding from the computational point of view. The linear analysis identified critical magnetic Reynolds numbers for the onset of the instability and discussed possibilities for its experimental detection (e.g., [24,26,27,28]). In particular, it was shown that the linear growth rate is determined by magnetic Reynolds number and Lundquist number (i.e., ratio of the magnetic diffusion time to the Alfvén crossing time), which should be both larger than unity for SMRI to operate. The saturation properties of SMRI following its exponential growth was investigated both in the weakly nonlinear regime, i.e., near the instability threshold, [29,30,31,32] and in the fully nonlinear regime [33,34,35,36,37,38,39]. It was demonstrated that SMRI saturates on the resistive time scale by modifying (reducing) the background velocity shear responsible for it and strengthening the background field. The dependence of the saturation level on the imposed field as well as on the magnetic Reynolds number, Prandtl number and Lundquist number was explored in greater detail.
While SMRI has been explored quite extensively both in TC flows and astrophysical disk context since the 1990s, its relative -helical magnetorotational instability (HMRI) -has become a subject of active theoretical and experimental research only in the last decade. These studies were initiated by Hollerbach & Rüdiger [40], who realized that adding an azimuthal magnetic field to the axial field can destabilize highly resistive flows at much smaller field strengths and at several orders of magnitude smaller Reynolds numbers, reducing the critical value from Re ∼ 10 6 needed for SMRI with purely axial field to Re ∼ 10 3 . Because this instability takes place in the presence of helical magnetic field, composed of azimuthal and axial components, it was termed HMRI. Like the SMRI, the HMRI also draws free energy from the background shear flow and transports angular momentum outward. Subsequently, HMRI has been widely studied theoretically by means of linear modal stability analysis [41,42,43,28,44,45,46,47]. It was shown that HMRI represents a destabilized inertial wave [41] and its growth rate is governed by the Reynolds number and the Hartmann number, Ha, characterizing the magnetic field strength. As a result, HMRI can persist even at extremely small magnetic Prandtl numbers typical to liquid metals, in contrast to SMRI, which is generally suppressed under this condition. Another difference with SMRI is that HMRI is restricted to rotational profiles with comparably steep negative shear or extremely steep positive shear, smaller than so-called the lower Liu limit or larger than the upper Liu limit, respectively [41,44]. Although the transition between HMRI and SMRI is monotonic when decreasing the ratio of azimuthal to axial field [40], the mathematics of this connection turned out to be quite subtle, including the formation of a spectral exceptional point where the original inertial mode coalesces with the slow magnetocoriolis mode [28]. Moreover, the recent nonmodal analysis revealed a fundamental connection between the modal growth of HMRI and the nonmodal dynamics in the corresponding HD problem [48]. The relevance of HMRI for astrophysical Keplerian disks has also been considered: it seems to be a promising candidate to replace SMRI in weakly ionized parts of disks with small P m, for example in the "dead zones" of protoplanetary disks or the outer parts of accretion disks around black holes. However, this issue is not yet fully resolved and under debate, because under standard conditions Keplerian shear alone might not be steep enough for HMRI to work [41] unless other physical factors, such as specific electrical boundary conditions [49,50] or additional axial currents within the fluid are involved [45].
The remarkable feature of HMRI to survive at very small magnetic Prandtl numbers and to set in at moderate Reynolds numbers makes it an ideal playground for experimental studies with liquid metals. Indeed, shortly after its theoretical discovery, a series of specially designed TC experiments [20,21] provided the first experimental evidence of HMRI at the liquid metal facility PROMISE and reproduced the main results of the linear theory, such as the stability threshold, the wavenumber and frequency of the HMRI-wave. In the experiments, the saturation amplitude as well as the propagation speed of the HMRI wave were measured in detail as a function of the system parameters (Re, Ha, ratio of rotation frequencies of outer and inner cylinders, etc.) and some differences with the numerical results were pointed out (see e.g., figure 10 in [21]).
This prompted further theoretical studies of HMRI in TC flows both with and without endcaps, focusing more on its nonlinear development and saturation. However, even today this is still a less explored area than the linear HMRI. First axisymmetric numerical simulations explored the dependence of the established HMRI-wave amplitude and propagation speed on the Reynolds number and the rotation ratio of the inner and outer cylinders as well as the role of endcaps [51,52,53]. However, the parameter range adopted in these studies was much narrower than that in the above experiments, being near the linear instability threshold of HMRI, when its dynamics is weakly nonlinear (see also [32]). More recent axisymmetric simulations [54,55] looked into the intrinsic nonlinear dynamics of HMRI in an axially infinite TC domain, avoiding complications because of the endcaps [56]. Different types of nonlinear regimes were analyzed and the applicability of the generalized quasi-linear approximation was tested by tracing the dynamics of large-scale Fourier modes. Despite these efforts, which undoubtedly contributed to a better understanding of the nonlinear dynamics of HMRI, detailed physics of its saturation and sustenance is still missing, especially when comparison with experiments is concerned. Evidently, it is this nonlinear saturated state which is observed in experiments, while the initial transient phase is much harder to identify due to the limited sensitivity of the applied ultrasonic flow measurement techniques.
In some respect, more progress has been made recently on the nonlinear dynamics of a "sibling" of HMRI -the azimuthal magnetorotational instability (AMRI). This nonaxisymmetric instability, which emerges in the presence of an imposed purely azimuthal fields in TC liquid metal flows, had been first identified after HMRI [22,57] (see also [46] and references therein). Using numerical simulations, Guseva et al. [58,59,60] probed much broader ranges of Reynolds and Hartmann numbers than those done for HMRI so far and identified different regimes of nonlinear saturation, from supercritical Hopf bifurcation near the linear instability threshold up to a catastrophic transition to spatiotemporal defects, which are mediated by a subcritical subharmonic Hopf bifurcation, and ultimately to turbulence. The scaling of the angular momentum transport of AMRI in the nonlinear regime with respect to Reynolds, Hartmann and Prandtl numbers was also explored in these papers.
Motivated by the experimental results of [21] and by the recent progress on understanding the nonlinear dynamics of AMRI [58], in this paper we investigate the evolution of HMRI, from its linear exponential growth to nonlinear saturation in liquid metal TC flows at small magnetic Prandtl numbers using numerical simulations. We consider infinite cylinders with periodic boundary conditions in the axial direction in order to avoid complications because of the endcaps (Ekman pumping) and concentrate on the intrinsic dynamics of HMRI driven by the combination of an imposed helical magnetic field and differential rotation of the flow. Distinctly from the above-mentioned previous studies on nonlinear HMRI, we do not restrict ourselves to axisymmetric perturbations only and allow non-axisymmetric modes to naturally develop during dynamical evolution, so that we can quantify the degree of non-axisymmetry of the saturated state of HMRI depending on the system parameters. It is well known that HMRI is axisymmetric in the linear regime, i.e. the most unstable modes do not have any azimuthal variation [40,42,44,54], and it is therefore important to know if this feature is retained in the nonlinear regime, which in turn can shed light on the saturation mechanism. We show that different regimes of the saturation are realized in the flow and analyze the characteristics of these states, such as angular momentum transport and energy spectra. This study, aiming at understanding the basic nonlinear dynamics of HMRI -underlying mechanism and properties of its saturation -is intended to be guiding and preparatory for the upcoming liquid sodium TC experiments in the DRESDYN facility at HZDR, which will combine and enhance the previous experiments on HMRI [21], AMRI [22] and on the current-driven kink-type Tayler instability [61].
The paper is organized as follows. The basic equations are introduced in section 2. Direct numerical simulations of the nonlinear saturation and evolution of HMRI as well as the spectral characteristics of the nonlinear state are presented in section 3. Summary and conclusions are given in section 4.

Main equations
The basic equations of non-ideal MHD governing the motion of an incompressible conductive fluid with constant kinematic viscosity ν and Ohmic resistivity η are ∂u ∂t where ρ is the constant density, p is the thermal pressure, u is the velocity, B is the magnetic field and µ 0 is the magnetic permeability of vacuum. The equilibrium state is an axisymmetric cylindrical magnetized TC flow between two coaxial cylinders with the inner radius R i and the outer radius R o , rotating with the angular velocities Ω i and Ω o , respectively. The flow is threaded by an externally imposed helical magnetic field B 0 = (0, B 0φ (r), B 0z ) consisting of constant axial, B 0z , and radially varying current-free azimuthal, B 0φ (r) = βB 0z R i /r, components in cylindrical coordinates (r, φ, z), where β is the dimensionless parameter characterizing the helicity of the field. We ignore the effects of endcaps in order to gain insight into the basic/intrinsic nonlinear evolution of HMRI. The ratio of the radii a ≡ R i /R o = 0.5 and the height of the cylinders, L z = 10d, where d = R o − R i = R i is the gap width, are chosen as in PROMISE [21]. An unperturbed flow between the cylinders, as follows from Eqs.
(1)-(3), is a standard vertically uniform and axisymmetric TC flow, U 0 = (0, rΩ(r), 0), with the angular velocity given by pressure p 0 (r) and constant density ρ 0 , satisfying the hydrostatic balance equation ρ 0 rΩ 2 = dp 0 /dr. The rotation ratio of the cylinders is fixed to µ = Ω o /Ω i = 0.27, slightly larger than the Rayleigh line µ c = 0.25 at R i /R o = 0.5, ensuring that purely HD instabilities are excluded, so that the flow can become unstable solely due to the presence of the magnetic field. In the following, we introduce the nondimensional variables by using the gap width d as the unit of length, Ω −1 i as the unit of time, Ω i d as the unit of velocity and ρ 0 Ω 2 i d 2 as the unit of pressure and energy density. (Since in the present case d = R i , this velocity scale is in fact the rotational velocity of the inner cylinder, Ω i R i .) The magnetic field is normalized by the imposed axial field B 0z . The Reynolds number is defined in terms of the rotation rate of the inner cylinder Its value is chosen to be Re = 6000 throughout the paper, which is within the range used in the related experiments [21]. In this paper, we consider a highly resistive fluid with very small magnetic Prandtl number, P m = ν/η = 1.4 × 10 −6 , typical to the liquid metal alloy GaInSn, so that the corresponding magnetic Reynolds number of the flow is also small, Rm = Re · P m = 8.4 × 10 −3 . The strength of the imposed axial field is measured by the Hartmann number, whereas the strength of the azimuthal field relative to the vertical one is measured by the helicity β parameter introduced above. Everywhere below this parameter will be fixed to β = 4.53, as adopted in some of the experiments [21]. Previous linear analysis showed that HMRI is effective for relatively strong azimuthal fields, β 1 (e.g., [40,41,28,46]). Consider perturbations of the velocity, pressure and magnetic field about the equilibrium, In the limit Rm ≪ 1, which applies here, the magnetic field perturbation induced by the perturbed flow is much smaller than the imposed field and scales with Rm [62,63,41], so we change it to b ′ → Rm · b ′ , such that the new b ′ is comparable to the perturbed velocity u. Substituting this into Eqs.
(1)-(3), taking into account that the imposed magnetic field is current-free, and using the above normalization, we arrive at the following equations governing perturbations with arbitrary amplitude in nondimensional form to first order in Rm (primes henceforth will be omitted): where the terms proportional to Rm on the left hand side of induction equation (6), including the time-derivative, are usually neglected in the limit of small Rm → 0 (inductionless, or quasi-static approximation, see e.g., [63,51,54]). In that case, the advective derivative, or the inertial term involving the perturbed velocity, (u·∇)u, is the only nonlinearity that remains in these equations. Here we preferred to keep these small terms for completeness. In the chosen nondimensional units, the rotational frequency (4) of the background velocity becomes and the imposed field is B 0 = (0, β/r, 1). Since η, µ, β and Re are fixed, only Ha is a free parameter, entering Eqs. (5)-(7) through the Elsasser number, or the interaction parameter: which is composed of the imposed magnetic field, rotation and resistivity and is independent of viscosity. The interaction parameter has been shown to play a special role in the linear dynamics of the HMRI, as it determines its modal growth rate, the corresponding critical wavenumber and other characteristics [44,46,48]. Specifically, in the inviscid case (Re → ∞), the HMRI can exist in a certain range of Λ and disappears at sufficiently large or small values of this parameter. In the first limit, the magnetic field is too strong and stabilizes the flow, whereas in the second limit it is too weak to support the instability. At finite Re, this implies that for a fixed Ha, HMRI operates in a range of Re and for a fixed Re, in a range of Ha [42,43,47], which was also confirmed experimentally [64,21]. Since Λ mainly determines the linear dynamics of HMRI, it is natural to explore how the nonlinear dynamics and saturation properties of HMRI depend on this parameter, which is our main goal. This is the main reason why we decided to follow the dynamics as a function of Λ. We anticipate that, like in the linear theory, this parameter will be decisive in the nonlinear outcome of the HMRI. Previous related nonlinear studies [51,53] focused mostly on the variation of the saturated state's characteristics with the Reynolds number at fixed Hartmann number. Besides, for these parameters, HMRI was in the weakly nonlinear regime. In this connection, we mention that the interaction parameter also plays an important role in the dynamics of forced MHD turbulence with mean field in the small-P m regime [62].

Numerical method
We solve the basic Eqs. (1)-(3) using the pseudo-spectral code from [58]. It is based on the Fourier expansion in the axial z-and azimuthal φ-directions and the finitedifference method in the radial direction. The nonlinear terms are calculated using the pseudospectral method and are de-aliased using the 3/2-rule. Integration in time is done with a second-order scheme based on the implicit Crank-Nicolson method. Although the inductionless approximation holds in the present case of small P m, the code updates the induction Eq. (2) in time in a general manner without invoking this approximation. Further details of the code and its validation tests on HD and MHD problems in TC flows can be found in [58]. The cylindrical flow domain in nondimensional units spans the radial range (R i , R o ) = (1, 2) divided into N r = 100 grid points, azimuthal range (0, 2π) and axial range (0, L z ) = (0, 10) with, respectively, N φ = 80 and N z = 400 spatial Fourier modes. The radial grid points are placed at collocation points of Chebyshev polynomials (despite the code is based on finite differences in r) that ensures higher resolution near the cylinder surfaces. Thus, although the HMRI is thought to be predominantly axisymmetric at least in the linear regime, we still allow for non-axisymmetric modes to naturally develop in the nonlinear regime by encompassing azimuthal wavenumbers, m, in the range (−N φ /2, N φ /2). A resolution test we performed (not shown here) showed that the adopted resolution is well sufficient for the present problem, capturing a whole range of wavenumbers from dynamically important intermediate ones, corresponding to the most unstable HMRI modes, down to wavenumbers dominated by viscous dissipation, as indicated by the energy spectra calculations in section 3.3 below. The radial boundary conditions at the cylinders are standard no-slip for the velocity perturbation, u| r=R i ,Ro = 0, and insulating for the magnetic field perturbation, meaning that current does not penetrate into cylinders, that is, ∇ × b = 0 at r < R i and r > R o . These radial boundary conditions are implemented in the code (see details in [63,65,58]), while periodic conditions are imposed for all variables in the axial zdirection, with a period L z , and in the azimuthal φ-direction. In this first effort to gain a deeper insight into the process of the nonlinear saturation of HMRI, we vary the Λ number at fixed Re and other parameters of the flow (which is equivalent to varying Ha).

Energy and torque
The perturbation energy density, composed of kinetic and magnetic parts, e = (u 2 + Ha 2 P m · b 2 )/2 (here b is the normalized magnetic field perturbation as defined in Eqs. [5][6][7], and torques at the cylinders, governing the evolution of angular momentum l = ru φ (in nondimensional units), are among the important characteristics/diagnostics of the flow. We see that at very small P m, the magnetic energy is much smaller than the kinetic one. The evolution equations for these quantities, taking into account the boundary conditions, are readily derived from Eqs. (5)- (7). As a result, we get for the total energy E = edV , to leading order in Rm [42,65] where D is the dissipation function The first term in Eq. (8) is the Reynolds stress u r u φ multiplied by the shear, rdΩ/dr, and describes energy exchange between the basic flow and perturbations, which is therefore due to the shear and is not affected by the magnetic field. The dissipation function D is always negative definite and consists of two parts. The first part describes the usual viscous dissipation of the kinetic energy, while the second part describes the Joule losses and in general is comparable to the former. Since the magnetic energy is negligible in this case, the Joule losses act as an additional effective dissipation of special kind for the kinetic energy through the kinetic-magnetic exchange term proportional to Λ in Eq. (5) [62]. Note that a net contribution from nonlinear terms has canceled out in Eq. (8) after averaging over the domain due to the boundary conditions. Thus, only Reynolds stress, when positive, can act as a source of perturbation energy, extracting it from the background flow due to the shear. The Maxwell stress, −b r b φ , does not play a role because of the smallness of the induced field in the low-Rm regime. The nonlinear term, not directly tapping into the flow energy and therefore not changing the total perturbation energy, acts only to redistribute energy among different wavenumbers.
In the quasi-steady saturated state, energy injection by the Reynolds stress balances dissipation losses. Since dΩ/dr < 0, this implies that the stress should be positive for the long-term sustenance of this state.
For the total angular momentum, L = ldV , we get where are the total viscous torques exerted, respectively, by the inner and outer cylinders on the flow. These torques also characterize angular momentum transport by the perturbations. Other -advection (i.e., Reynolds stress) and magnetic -torques vanish after volume integration due to the boundary conditions and do not contribute to the total angular momentum evolution. In the equilibrium (laminar) state, the torques at both cylinders have the same absolute value, implying that the inner cylinder tries to increase the angular momentum of the basic flow, while the outer cylinder to decrease. These two torques are in balance, resulting in the stationary rotational profile (4). Below we measure the torque in the perturbed state against the laminar torque, i.e., consider the ratio G i,o /G lam and redefine G i,o /G lam → G i,o . Its advantage is that it does not depend on the normalization of the velocity.

Results
We move to the analysis of the HMRI, from its linear growth to nonlinear saturation. To have an idea about the regimes of the onset and operation of HMRI in the given TC flow configuration, in figure 1 we show the growth rate of the most unstable mode (which is axisymmetric) as a function of Λ and Re, obtained by solving the linear stability problem with the same parameters and boundary conditions. From this figure we see that HMRI operates in a certain interval of Λ provided the Reynolds number is larger than a critical value, Re > Re c = 1244. The instability interval with respect to the interaction parameter broadens with the increase of Re, but converges in the inviscid limit (not shown in this figure). It is also seen in figure 1 that at sufficiently high Re 3000, the maximum growth of HMRI is mainly determined by Λ and attained at its nearly constant value. For the chosen Re = 6000, the instability first appears in the flow at the critical value Λ c = 0.022, reaches a maximum at Λ = 0.15 and then decreases again, disappearing at Λ = 1.38. Having outlined the main linear features of HMRI, we can now move to nonlinear evolution.
We initialize the simulations by imposing random noise perturbations of the velocity and magnetic field on the equilibrium flow. The subsequent evolution of the total energy E and torques G i , G o are shown in figure 2 at different interaction parameters, which are all above the critical value Λ c = 0.022 when HMRI first emerges. (We also checked that, as expected, the runs at Λ < Λ c eventually decay, although a further study with a variable initial amplitude of perturbations is required to see whether there is a possibility of subcritical transition in this case.) As we will see below, the dynamics  in these cases differ qualitatively not only in terms of the temporal evolution, but also in its spatial appearance and spectral characteristics. After an initial transient phase, the most unstable mode eventually emerges and grows. The exponential growth phase lasts until the velocity amplitude becomes large enough for the nonlinear term in Eq.  After the initial exponential growth phase, these quantities settle down into a sustained quasi-steady state, where they oscillate around constant values. and torque increase with increasing the interaction parameter. This behavior of the saturation level with Λ will be investigated below in more detail. Because the volumeaveraged angular momentum is also quasi-stationary in time, both torques have nearly the same absolute values, although the inner torque is more oscillatory at larger Λ due to the appearance of smaller scale turbulent eddies near the inner cylinder (see below). This implies that the HMRI, like the SMRI, transports angular momentum outwards, which is expected, since it also derives energy from shear.

Saturation mechanism
Let us see what is a main mechanism underlying the saturation of HMRI. During its initial exponential growth, the main nonlinear term (advective derivative in Eq. 5) also gradually gains strength and transfers the energy from the most unstable wavenumbers to larger and smaller ones. In the saturated state, the linear energy extraction rate due to the Reynolds stress associated with the unstable wavenumbers matches the nonlinear transfer rate from these wavenumbers to larger ones, where energy is  ultimately dissipated due to viscosity. Together with Λ, the linear growth rate of HMRI is determined by the radial shear of the azimuthal velocity characterized by the Rossby number, Ro = (r/2Ω)dΩ/dr (e.g., [41,46]), which is a function of radius in a TC flow (for a Keplerian flow Ro = −0.75 everywhere). In the exponential growth phase, a slowly increasing feedback of the nonlinear term on the large-scale azimuthal velocity (which is initially a TC profile given by Eq. 4) rearranges the radial distribution of shear (Ro), such as to reduce it by absolute value in the bulk of the flow. This is clearly illustrated in figure 3, which shows the radial profile of Ro, calculated with the total azimuthal velocity (i.e., initial TC flow plus the azimuthal velocity perturbation averaged over azimuthal φ− and axial z-directions) in the saturated state. It considerably deviates from the initial profile corresponding to the TC flow and results in the reduction of the linear growth rate of the HMRI and hence of the energy extraction rate from the mean flow. This is confirmed by figure 4, which shows a parallel evolution of energies in two simulations one starting with the saturated/modified mean azimuthal velocity profile for the run shown in figure 3   Calculating the inclinations of these two linear curves, we infer that the latter velocity profile indeed results in an 3.77 times reduced growth rate of HMRI than the former.
(4) for the same other parameters of the flow, following only early exponential growth phase in the linear regime. Comparing the growth rates, derived from the inclination of ln(E) versus time, we clearly see that the saturated profile indeed results in a reduced exponential growth rate for the perturbation energy in the initial linear regime, which for the chosen value of Λ is smaller than that in the presence of originally imposed TC flow by a factor of 3.77. Saturation sets in when the energy gain by the linear instability, being reduced by the nonlinearity, eventually becomes equal to the nonlinear transfer of the energy. The saturation levels in these two cases also differ, because, as we checked, the corresponding established nonlinear states have different structures -the number of rolls in the axial direction. Thus, there is no unique flow attractor, which would otherwise result in the same saturation amplitudes for these two profiles. Obviously, the deviation of the shear profile from the initial one due to the nonlinear feedback is larger, the higher is Λ, i.e., the more unstable is the flow. So, the saturation, where the energy gain by the linear HMRI is in balance with the nonlinear transfers, occurs through the modification (reduction) of the shear of the mean rotational velocity, which in turn defines the growth rate of HMRI. The modified shear profile in the quasi-steady state then stays practically unchanged during the evolution. This type of saturation mechanism is in fact similar to that for SMRI in TC flow discussed in [33,30,36,31,32]. As for the initially imposed mean magnetic field B 0 , it basically does not change in the saturated state, because, as mentioned above, magnetic field perturbations are usually orders of magnitude smaller the imposed field in the low-Rm regime. Figure 5 shows the time-averaged values of the energy, E (a), the inner torque G i (b) and the dissipation rate D (c) in the quasi-steady saturated state as a function of Λ. Since this state is on average constant in time, the saturated mean values of both torques are essentially the same, G i = G o , at each Λ, as it should be (see also [60]). This figure allows us to clearly see different regimes of the nonlinear saturation of the HMRI in TC flow and transitions between them. When the interaction parameter increases, the instability first appears at Λ = Λ c = 0.022 via classical supercritical Hopf bifurcation [67]. For Λ > Λ c , but still near the instability threshold, the perturbations are weakly nonlinear and the saturated energy is proportional to ∝ Λ − Λ c , as expected for Hopf bifurcation. The saturated torque and dissipation function also exhibit a similar dependence on Λ as the energy. As mentioned above, in this case of |Λ − Λ c | << Λ c the growth rate is small and the saturation time is long, several tens of viscous time.

Properties of the saturated state
A typical spatial structure of the velocity in this regime is shown in figure 6(a), which consists of a well-organized chain of regular Taylor vortices uniformly filling the domain and resembles the corresponding most unstable mode from which it originated. As evident from the isosurface plot of the axial velocity, this weakly nonlinear state is axisymmetric. Five pairs of vortices fit in the domain, implying that the dominant most unstable azimuthal and axial wavenumbers are (m, k z ) = (0, 2πn z /L z ) with n z = ±5, which will be confirmed by spectra calculations presented in the next section. This picture is consistent with the pattern of the unstable modes of the HMRI observed in previous linear stability studies [40] and nonlinear simulations in the similar weakly nonlinear regime [51]. This regime holds until about Λ tr = 0.0232, where an abrupt transition to less organized and irregular (turbulent) regime takes place, which is marked by a sharp inflection point on the curves of the saturated values in figure 5. Beyond this point both E and G i exhibit a different behavior with Λ than the linear one in the weakly nonlinear regime. When Λ is just above Λ tr , they start with a slower increase, followed by a steeper monotonic increase starting from Λ = 0.025. The transition to the strongly nonlinear regime is reflected in the clear change of temporal behavior of the total energy and torques -oscillations emerge in their evolution, which become more and more irregular and of higher amplitude with increasing the effect of nonlinearity as Λ grows. This is illustrated in figure 7, which shows the time evolution of G i before (at Λ = 0.023) and after (at Λ = 0.026) the transition value Λ tr and also in figure 2 at higher Λ. The emerging oscillations in the later case are due to a wider spectrum of higher frequency modes (inertial waves, see below) excited in the stronger nonlinear regime. However, already after Λ ≃ 0.08, the saturated energy and torques do not seem to increase anymore and the character of the time-variation are more or less similar. In the next section, we will see that not only the time evolution properties, but also the power spectra of the saturated state of HMRI drastically and abruptly change when going from the weakly nonlinear to the strongly nonlinear regime. Figures 6(b)-6(d) show the spatial structures of the radial u r and axial u z velocities in the saturated state of HMRI at different Λ after the transition point, which differ qualitatively from those in the weakly nonlinear regime in figure 6(a). Now the eddies of different shapes and sizes have appeared in the domain, being distributed non-uniformly along the axial direction. With increasing Λ, more and more small-scale vortices emerge mostly near the inner cylinder, while larger ones remain near the outer cylinder (figures 6(c) and 6(d)). As a result, the flow takes the form of fully developed turbulence. This is most clearly seen in the maps of the axial velocity. These eddies change shape in a random manner on the turnover (dynamical) time, which is the shorter the smaller is the eddy. As a result, the torque at the inner cylinder, G i , displays faster chaotic oscillations due to the small-scale eddies than the torque at the outer cylinder,G o , due to larger scale eddies, as seen in figures 2(b) and 2(c). Note also how the axisymmetry of the nonlinear state changes with Λ (see also figure 9 below). After the transition point Λ tr , before the emergence of smaller-scale eddies at higher Λ 0.1, the nonlinear state preserves nearaxisymmetry, as evidenced by the isosurface plots of the axial velocity, for example, in the case of Λ = 0.05 in figure 6(b). So, we observe a type of self-sustained quasi two-dimensional (2D, i.e., dependent only on r and z coordinates) MHD turbulence driven/fed by HMRI, analogous to that previously found in hydrodynamically unstable TC flow with much stronger and purely azimuthal background magnetic field [68]. We calculate its spectral characteristics in the next section. However, with further increasing Λ, starting from about Λ = 0.07, the nonlinear state slightly deviates from axisymmetry, as shown in figures 6(c) and 6(d), respectively, for Λ = 0.1 and 0.2 (see also figure 9 Increasing Λ, the saturated state changes from (a) weakly nonlinear regular HMRI wave to the turbulent state with different structures, dominated by (b) larger scale eddies or (c),(d) by many smaller scale eddies near the inner cylinder and larger ones near the outer cylinder. All these states are nearly axisymmetric, although the contribution of non-axisymmetric modes are noticeable in (c) and (d) cases and are mainly attributable to smaller scale eddies. below). One can discern in these plots that the non-axisymmetry is associated mostly with smaller-scale eddies near the inner cylinder, while the larger scale ones near the outer cylinder stay mostly axisymmetric. Noticing that the saturated mean azimuthal velocity profile in figure 3 develops a Rayleigh-unstable region, Ro < −1, near the inner cylinder at this larger Λ, we may conjecture that the emergence of small-scale eddies (spatial defects) in the same radial interval is due to a HD instability of the nonlinearly modified flow profile. The spectral analysis presented below will give a more quantitative characterization of the degree of the non-axisymmetry in this case. In this regard, we would like to caution that studying the nonlinear dynamics of HMRI at any values of the Elsasser number (i.e., for any strength of the imposed magnetic field) with purely axisymmetric simulations might lead to the incomplete picture of the saturation and nonlinear dynamics of HMRI, overlooking the role of these modes. An advantage of the present analysis, including non-axisymmetric modes, is that it enables us to establish at which magnetic field strength (i.e., interaction parameter) the nonlinear saturated state of HMRI deviates from axisymmetric configuration (see below). The transition from weakly nonlinear to strongly nonlinear and ultimately turbulent regime for HMRI is actually analogous to the nonlinear transition found for AMRI with increasing the Hartmann number of the imposed purely azimuthal field at constant Re by Guseva et al. [58]. (Since Reynolds number is fixed, increasing the interaction parameter corresponds to the increase of Hartmann number.) As in our case, in the weakly nonlinear regime, a regular pattern of vortices emerge in the saturated state via Hopf bifurcation, except that for HMRI it is a traveling wave, as will be evident from spectral analysis below, while for AMRI it is a standing wave. With further increase of Hartmann number, so-called defects appear. The spatial structure at this stage is similar to that shown in figure 6(b) with nonuniformly distributed irregular vortices, not far from the transition point. Subsequently, a fully developed turbulent state arises from these defects at even higher Hartmann numbers -a stage which can be identified in our case with that depicted in figures 6(c) and 6(d). It is also expected that a similar hysteresis phenomenon mediated by an edge state dividing weakly nonlinear and strongly nonlinear (turbulent) states, as found by Guseva et al. for AMRI, also exists for HMRI. However, we did not study this possibility here: the system is always followed starting with the same random initial conditions when changing the interaction parameter. In this case, after the exponential growth, the flow always ends up in the chaotic state if Λ is larger than the transition value. Identifying edge states and the hysteresis in the transition process for HMRI will be a subject of future study. Such edge states are known to play an important role in the turbulence transition in HD shear flows (see [69] and references therein).

Spectral characteristics
Further insight into the nature of the quasi-steady states of the HMRI at different interaction parameters, which we have described so far in physical space, can be gained by spectral analysis. It is well known from linear theory that in the case of the background helical magnetic field, z-reflection symmetry is broken [67,40,46] and, consequently, in the given right-handed configuration of the axial and azimuthal fields (ΩB 0φ B 0z > 0), HMRI is characterized by inertial waves that travel downwards (opposite the z-axis) [41,46]. As a result, the subsequent nonlinear state is composed of these downward propagating waves. In a real experiment, both upward and downward propagating waves are present because of the reflection from the endcaps, but the one that is unstable tends to be dominant [21]. Below we analyse the axial and azimuthal spectra of these waves. To this end, we decompose each velocity component in the azimuthal and axial directions (in which they are periodic), where the index i = r, φ, z and m = 0, 1, 2, ... and k z = 2πn z /L z with n z = 0, ±1, ±2, ..., are the azimuthal and axial wavenumbers, respectively. We also define radially integrated spectral energy density, as well as the azimuthal,Ê m , and axial,Ê, spectral energy densities viâ so that their sum over wavenumbers is equal to the total energy, E = m,kz E(m, k z ) = mÊ m = kzÊ . These energy spectra are among the main characteristics/diagnostics in the nonlinear (turbulent) regime in flows. In particular, the azimuthal spectrum characterizes the total contribution of different m in the dynamics and hence can serve as a measure of the non-axisymmetry, while the axial k z spectra provide information on the energy-injection due to axisymmetric HMRI modes (see e.g., [68]). We use For reference, we also show the 2D HD turbulence spectrum k −3 z as well as a flatter k −2 z spectrum of rotating HD turbulence. At small Λ, before the transition value Λ tr = 0.0232, the spectrum is discrete ("spiky" dashed lines), dominated by the most unstable mode and its multiple wavenumbers, while for Λ beyond the transition point, in the strongly nonlinear (turbulent) regime, the spectrum (solid lines) becomes continuous with all axial wavenumbers excited; at Λ 0.05 it exhibits a scaling behavior at intermediate wavenumbers 5 k z 20, which is close to k −3 z .
these spectra below to characterize the dynamics in the quasi-steady state. This state, however, exhibits, as we have seen above, irregular oscillations, making the spectra noisy. To avoid this, we also averageÊ m andÊ in time over the whole duration of the quasi-steady state in the simulations. Figure 8 shows the time-averaged axial energy spectrumÊ in the saturated state at different Λ. It is seen that qualitatively different spectra correspond to weakly and strongly nonlinear saturation regimes, and the rearrangement from the former to the latter type, like the saturated values of energy and torque (figure 5), occurs suddenly near the transition point Λ tr = 0.0232. At small Λ Λ c , the linearly most unstable axial wavenumber of HMRI, k z = 2πn z /L z = 3.14 with the mode number n z = 5, is dominant also in the saturated state and carries most of the power. In physical space, this corresponds to five pair of vortices, as shown in figure 6(a). The modes with wavenumbers, which are multiples of the most unstable axial wavenumber, are excited due to weak nonlinear self-interaction of the dominant mode in the saturation process. These subharmonics (separate spikes in figure 8) have orders of magnitude smaller power compared to the dominant one, whereas at all other k z the spectral power is essentially zero. Remarkably, the energy spectrum changes qualitatively when the interaction parameter approaches a transition value Λ tr , where the effect of nonlinearity is more appreciable. As a result, power in other wavenumbers, lying between the multiple wavenumbers, abruptly increases, i.e., the spectrum is quickly filled and becomes continuous, although the multiple harmonics of the most unstable wavenumber still have larger power than other ones. But already at Λ = 0.025, slightly larger than the transition point, the spectrum takes a smoother shape, continuously well populated at all wavenumbers, as it is characteristic of a turbulent state. With further increasing Λ, this spectrum moves upwards and converges starting from about Λ = 0.05. Now, the maximum comes again at the wavenumber of the most unstable HMRI mode, k z ∼ 2 − 3, which mainly determines the energy injection from the basic flow into turbulent fluctuations. As it is seen in figure 8, at intermediate wavenumbers in the range 5 k z 20, this converged spectrum displays the power-law dependence close to k −3 z typical for 2D HD turbulence [70], as opposed to the Kolmogorov spectrum, k −5/3 , for three-dimensional turbulence or to the k −2 spectrum for rotating HD turbulence [71]. Then, at k z 20 the spectrum decays faster because of viscous dissipation. The appearance of the quasi-2D turbulence in physical space, corresponding to these spectra, we have already seen in figure 6. The energy spectrum with a similar powerlaw dependence were also reported for 2D MHD turbulence in TC flow with a large purely azimuthal magnetic field in the highly resistive (inductionless) limit, Rm ≪ 1, in the above-mentioned work [68]. In that case, however, the turbulence is due to HD instability and forced to be 2D by the imposed azimuthal field. By contrast, in the present case, where the TC flow is hydrodynamically stable, the observed turbulence is of magnetic origin, triggered and energetically supplied by HMRI, and hence quasi-2D from the outset. This implies that the turbulence we observe here cannot be described within standard Kolmogorov phenomenology, which anyway is inapplicable in the case of high resistivity, because the Joule dissipation in this regime is anisotropic and acts at all scales in the flow (see e.g., [62]). In this situation, it is not possible generally to define an inertial range in a classical sense, where only nonlinear term (u · ∇)u in Eq. (5) operates and transfers energy towards large wavenumbers. The presence of the power-law interval in the energy spectrum different from Kolmogorov one that has been found here and in [68] is a confirmation of that. The similarity of the spectra and the associated power-law indices, despite different driving mechanisms of the turbulence in our case and in [68], indicates that these spectral features could be generic to quasi-2D MHD turbulence, which can occur in a magnetized resistive TC flow.
As we have found above, the saturated state remains nearly axisymmetric, despite non-axisymmetric modes emerging with increasing Λ. To quantify the role of the latter modes, in figure 9 we plot the ratio of the sum of azimuthal energiesÊ m of all the non- axisymmetric modes with m ≥ 1 to the energy of the axisymmetric modes with m = 0, measured by the parameter A = m≥1Ê m /Ê 0 , as a function of Λ (a) as well as the typical azimuthal energy spectrum in the strongly nonlinear state (b). From this figure, it is seen that non-axisymmetric modes are essentially absent at small Λ; they start to emerge from about Λ = 0.07, after which the total energy of all non-axisymmetric modes relative to that of axisymmetric ones increases with Λ. However, the former still remains much smaller than the latter, particularly in the range 0.07 ≤ Λ ≤ 0.2, where the saturation level reaches a maximum ( figure 5). This is also confirmed by the azimuthal spectrum of the energy,Ê m , in figure 9(b). Since non-axisymmetric modes are present only at Λ 0.07, we show the azimuthal spectrum at typical values Λ = 0.1 and 0.2, corresponding to the strongly nonlinear (turbulent) regime. It reaches a maximum at m = 0 and rapidly decreases with m. The energy of the first non-axisymmetric modes with m = 1 is already by more than an order of magnitude smaller than that of the axisymmetric ones,Ê 1 /Ê 0 =0.007 and 0.04, respectively, at Λ = 0.1 and 0.2, and the energy of higher non-axisymmetric modes is even lower.

Summary and conclusion
In this paper, we investigated the development of HMRI in an infinite/periodic TC flow domain at very small magnetic Prandtl numbers by following its evolution from the linear growth phase to nonlinear saturation using direct numerical simulations. To focus on the basic dynamics of HMRI arising from the combined effect of differential rotation and helical magnetic field, we simplified the analysis by ignoring the effects of endcaps that would induce Ekman pumping. We analyzed the dynamics with respect to the interaction parameter, or Elsasser number, Λ, which is known to be a central parameter determining the linear evolution of HMRI. As distinct from previous nonlinear analyzes of HMRI, which focused only on axisymmetric modes, an advantage of our study is that by allowing for non-axisymmetric modes in the simulations, it enabled us to establish for which interaction parameters the role of the latter modes becomes important. We confirmed that in the nonlinear regime, just as SMRI, HMRI transports angular momentum outward. We demonstrated that different regimes of nonlinear saturation are realized when the interaction parameter changes. At smaller values of this parameter, just above the stability threshold, the HMRI saturates in the weakly nonlinear regime and the corresponding spatial structure consists of well-organized regular vortices, which are axisymmetric. In this case, the most unstable mode of instability dominates, which is also supported by our spectral analysis. However, with increasing the interaction parameter, at a certain value an abrupt transition to strongly nonlinear (turbulent) state takes place, marked by different, irregularly oscillating behavior of the energy and stress. The saturated values of the energy and stresses as a function of Λ exhibit a sharp inflection point, corresponding to this transition, and then increase with until about Λ = 0.1 and then slowly decrease (figure 5), since the effectiveness of the HMRI itself is decreasing at higher magnetic fields. The spatial structure of the nonlinear state is predominantly axisymmetric and represents a type of self-sustained quasi-2D (in (r, z)-plane) MHD turbulence due to HMRI. Generally, classical 2D MHD turbulence is decaying and requires external forcing for long-term maintenance [72]. So, here we demonstrated the existence of self-sustained quasi-2D MHD turbulence in a TC flow at high resistivity, or small magnetic Prandtl numbers caused and supplied by HMRI at the expense of the flow energy. We also calculated the energy spectrum as a function of axial wavenumber and found that it increases with Λ, but converges from about Λ = 0.05, and at intermediate axial wavenumbers k z exhibits a power-low dependence close to k −3 z , as typical for 2D HD turbulence [70]. It would be interesting to probe in future studies whether this quasi-2D turbulent state persists as the Reynolds number increases or a transition to three-dimensional turbulence occurs.
In both weakly and strongly nonlinear regimes, the saturation mechanism appears to be general, consisting in the modification, or reduction of the radial shear profile of the mean azimuthal velocity in the bulk of the flow, which, in turn, results in the reduction of the exponential growth rate of HMRI to match energy transfer rate due to the nonlinear (advection) term. Such a nonlinear saturation mechanism was already discussed for SMRI in a TC flow [33,29,36,31,32].
The present analysis, although simplified by excluding the effects of endcaps, is a first step towards understanding the experimental manifestation of HMRI, where it is already in the saturated nonlinear regime. The insight gained from this study will be a stepping stone for future numerical studies incorporating endcaps and conducing boundaries, where the situation is complicated by Ekman circulations, penetrating from the endcaps in the bulk of the flow, and for subsequent comparison with planned experiments on HMRI within the DRESDYN project at HZDR. Nevertheless, as demonstrated in [58] on the example of AMRI, the results obtained with periodic boundary conditions along the axial direction are still useful to interpret the experiments. In the future, we plan to do analogous studies of HMRI in a TC flow with endcaps and compare with the present results. To date, such realistic three-dimensional simulations of HMRI, which are necessary for understanding the related experimental results, have not been undertaken yet because of high computational cost. Although previous early simulations explored the nonlinear development of HMRI in TC flow with endcaps, they considered a narrow range of parameters and/or were axisymmetric by design [51,52,53], not capturing different regimes of the nonlinear saturation.
Another interesting venue of research extending the present analysis is to explore the nonlinear dynamics of HMRI and associated angular momentum transport for quasi-Keplerian rotation relevant to protoplanetary disks, whose dense and cold interiors are too resistive for SMRI to operate. The nonlinear development of HMRI has not been explored also for positive shear profiles relevant to the near-equator strip of the solar tachocline. These topics are of a current research interest and have been studied so far only in the linear regime [45,46,73].