Transition to magnetorotational turbulence in Taylor--Couette flow with imposed azimuthal magnetic field

The magnetorotational instability (MRI) is thought to be a powerful source of turbulence and momentum transport in astrophysical accretion discs, but obtaining observational evidence of its operation is challenging. Recently, laboratory experiments of Taylor--Couette flow with externally imposed axial and azimuthal magnetic fields have revealed the kinematic and dynamic properties of the MRI close to the instability onset. While good agreement was found with linear stability analyses, little is known about the transition to turbulence and transport properties of the MRI. We here report on a numerical investigation of the MRI with an imposed azimuthal magnetic field. We show that the laminar Taylor--Couette flow becomes unstable to a wave rotating in the azimuthal direction and standing in the axial direction via a supercritical Hopf bifurcation. Subsequently, the flow features a catastrophic transition to spatio-temporal defects which is mediated by a subcritical subharmonic Hopf bifurcation. Our results are in qualitative agreement with the PROMISE experiment and dramatically extend their realizable parameter range. We find that as the Reynolds number increases defects accumulate and grow into turbulence, yet the momentum transport scales weakly.


I. INTRODUCTION
The magnetorotational instability (MRI) is of great importance in astrophysics. First discovered by Velikhov [1] in 1959, it remained unnoticed until 1991 where Balbus and Hawley [2] realised its application to accretion disc theory. Accretion discs are astrophysical systems that consist of ionised gas and dust orbiting a massive body. Planets and stars are formed from this initially dispersed matter. The physical mechanism of matter contraction is straightforward: a parcel of viscous fluid in the differentially rotating disc loses its angular momentum over time and falls onto the central object. To explain the astrophysically observed rates of accretion, one must assume a turbulent transport of angular momentum in the outward direction [3]. In so-called Keplerian discs the angular velocity profile of gas follows the law which is hydrodynamically stable according to the Rayleigh criterion for rotating fluids [4]: Ionized accretion discs, however, are necessarily magnetized and the MRI acts in rotating flows if the angular velocity decreases with radius, which is true of Keplerian flows (1).
The growth rates of the MRI and the parameter ranges in which it acts were determined in several linear analyses [2,[5][6][7][8], but these cannot provide information about the flow structure and scaling of angular momentum transport after nonlinear saturation. In the last two decades there has been a great deal of numerical work concerned with the nonlinear properties of the MRI. Simulations were usually performed with the shearing sheet approximation, which is a local model of an accretion disc [6]. The main disadvantage of this model is the influence of boundary conditions on the geometry of the observed modes and transport scaling [9].
These theoretical results and numerical simulations inspired physicists to realise the MRI in laboratory experiments. In 2001 Ji et al. [10] and Rüdiger and Zhang [11] independently suggested the possibility of directly observing magnetorotational instabilities in a cylindrical vessel made of two co-axial and independently rotating cylinders containing a liquid metal alloy (see Fig. 1a). The standard form of the MRI (SMRI), which they proposed, emerges when a purely axial magnetic field is imposed, but this has not yet been achieved in ex-periments. The difficulty is that liquid metals have very small magnetic Prandtl numbers (e. g. P m ∼ 10 −6 for gallium alloys), which implies that very large Reynolds numbers are necessary to observe the SMRI (Re 10 7 ). In fact, such high Reynolds numbers have never been achieved even for non-magnetic Taylor-Couette flows. Furthermore, we note that endplates cause hydrodynamic instabilities in quasi-Keplerian Taylor-Couette flows even at low Reynolds numbers [12,13]. This makes it difficult to identify the SMRI unambiguously.
Hollerbach and Rüdiger [7] proposed instead a combination of axial and azimuthal magnetic fields, giving rise to the helical MRI (HMRI) at Re ∼ 10 3 and Hartmann numbers Ha ∼ 10. This was successfully observed [14,15] in the PROMISE facility (Potsdam-ROssendorf Magnetic Instability Experiment). Both the SMRI and HMRI consist of axisymmetric toroidal vortices, which are stationary for the former but travel axially for the latter. Hollerbach et al. [8] realised that although a purely azimuthal magnetic field does not yield any axisymmetric instabilities, non-axisymmetric modes can be destabilized. The resulting azimuthal MRI (AMRI) arises in Taylor-Couette flow at Re ∼ 10 3 and Ha ∼ 10 2 [8], and a recent upgrade of the PROMISE power supply made its experimental observation possible [16].
Despite recent experimental progress in realizing magnetorotational instabilities in the laboratory, little is known about the bifurcation scenario, transition to turbulence and transport properties. In this work we address these points for the AMRI. We perform direct numerical simulations of the coupled induction and Navier-Stokes equations using axially periodic boundary conditions, thereby avoiding undesired end-plate effects. We find that the laminar quasi-Keplerian flow becomes unstable to a wave rotating in the azimuthal direction and standing in the axial direction. Subsequently, we identify a new bifurcation scenario giving rise to spatio-temporal defects via a subcritical subharmonic Hopf bifurcation. As the Reynolds number is further increased, the flow becomes turbulent and outward momentum transport is enhanced, albeit at a weak rate. The results are in good qualitative agreement with the PROMISE observations [16] and substantially extend the parameter range which can be explored experimentally. action of AMRI from [8]. The rotation rate is fixed at µ = 0.26. Green, blue and red lines show the values of Re and Ha numbers we study in this work.

II. GOVERNING EQUATIONS
We consider an incompressible viscous liquid metal that is sheared between two independently rotating cylinders of radii r i (inner) and r o (outer). The angular velocity of the cylinders are Ω i and Ω o , respectively, and an external azimuthal magnetic field (r i /r)B 0 , where r is the radial coordinate, is imposed. The relevant fluid properties are the electrical conductivity σ, the kinematic viscosity ν , the density ρ and the magnetic diffusivity η. The velocity field u is determined by the Navier-Stokes equations (3), whereas the magnetic field B is determined by the induction equation (4), which represents a combination of the laws of Ampere, Faraday, and Ohm. The equations were rendered dimensionless by using the gap between cylinders d = r o − r i for length, d 2 /ν for time, and B 0 for the magnetic field.
In dimensionless form they read together with ∇ · u = ∇ · B = 0. Here p is the pressure, Ha the Hartmann number, and P m the magnetic Prandtl number. The dimensionless parameters of the system are specified in Table I. Following the PROMISE experiment [16], we use a magnetic Prandtl number of P m = 1.4 · 10 −6 (corresponding to the alloy Ga 67 In 20.5 Sn 12.5 ), a radius-ratio of δ = 0.5 and a rotation-ratio of µ = 0.26. This places the velocity profile in the quasi-Keplerian regime, for which the angular velocity decreases radially, whereas the angular momentum increases, i.e. δ 2 < µ < 1.
Periodicity in the axial direction is imposed with basic length L z . The background circular Couette flow U = V (r)e φ is a solution to these equations and boundary conditions, given by The magnetic field is also assumed periodic in the axial direction. In the radial direction the boundary condition depends on the material of the cylinders. Typically two idealized cases are considered in the MRI problem: insulating and conducting cylinders. These lead to slightly different results, as theoretically demonstrated by Chandrasekhar [17]. However, the difference is not dramatic, and here we will consider only the case of insulating boundaries.
See Willis and Barenghi [18] for a derivation for a Fourier decomposition of the magnetic field. if the primary instability is a Hopf bifurcation the resulting states can be either standing or traveling waves in the axial direction [19]. By contrast, a combined helical magnetic field breaks the reflection symmetry and only traveling waves (TW) can be observed [20]. Finally, if the bifurcating solution is non-axisymmetric, as in the AMRI, this will generically be a rotating wave in the azimuthal direction.

III. NUMERICAL METHOD
In the numerical simulations only the deviation u ′ = u − V (r)e φ is computed. Its governing equations read which are supplemented with homogeneous boundary conditions u ′ = 0. Here N stands for the nonlinear term in the Navier-Stokes equations (3), which contains the advective terms and the Lorentz force: Spatial discretisation is accomplished via a Fourier expansion in the azimuthal and axial directions, As the variables are real, their Fourier coefficients satisfy the property A k,m = A * −k,−m , where A * denotes the complex conjugate. Finite differences are employed in the radial direction, on a grid of N points. Derivatives are calculated typically using 9-point stencils.
The equations for u r and u φ are coupled through the Laplacian. They can be separated in a Fourier decomposition, however, by considering where the ± are considered respectively. Using this change of variables the governing equations become where A. Influence-matrix method It is well known that primitive variable formulations are subject to loss of temporal order if inappropriate boundary conditions are enforced. For the velocity field the problem is associated with the lack of a boundary condition for the pressure. For the magnetic field, there appear at first sight to be too few boundary conditions, and further, the components of B are coupled in the boundary condition. It is shown here that the influence-matrix method resolves these issues, the appropriate boundary conditions can be satisfied to machine precision, and temporal order is retained. We show first how the method is applied for time integration of the velocity field.
The governing equations are solved using a 2nd-order predictor-corrector method, based on the implicit Crank-Nicolson method. The influence-matrix technique applied to the velocity is similar to [21]. An analogous approach is applied to the magnetic field.

Method for the velocity field
We write the time-discretised Navier-Stokes equations in the form where q denotes time t q . This form is sixth order in r for u q+1 and second order for p, without the solenoidal condition explicitly imposed. In principle this system should be inverted simultaneously for p and u q+1 with boundary conditions u q+1 = 0 and ∇ · u q+1 = 0. In practice it would be preferable to invert for p first then for u q+1 , but the boundary conditions do not involve p directly. Note that the Y u q term has been included in the right-hand side of the pressure-Poisson equation, so that it corresponds precisely to the divergence of the equation for u q+1 . This also ensures that any non-zero divergence in the initial condition is projected out after a single time-step. We split the system (14) into the 'bulk' solution, with boundary conditionsū = 0 and ∂ rp = 0, and the auxiliary systems with boundary conditions u ′ = 0 and with boundary conditions u ′ , i} or those swapped on r i and r o . The system (16) provides, with the two boundary conditions, two linearly independent functions u ′ j that may be added toū without altering the right-hand side in (15). Similarly the system (17) provides a further six functions. The superposition may be formed in order to satisfy the eight original boundary conditions, u q+1 = 0 and ∇ · u q+1 = 0 on r = {r i , r o }. Substituting (18) into the boundary conditions, they may be written where A = A(g(u ′ )) is an 8×8 matrix. The appropriate coefficients required to satisfy the boundary conditions are thereby recovered from the solution of this small system for a. The error in the boundary conditions g j (u q+1 ) using the influence-matrix technique is at the level of the machine precision.
The auxiliary functions u ′ j (r), the matrix A and its inverse may all be precomputed, and the boundary conditions for u ′ have been chosen so that u ′ ± are purely real, u ′ z is purely imaginary, and A is real. For each timestep, this application of the influence matrix technique requires only evaluation of the deviation from the boundary condition, multiplication by an 8×8 real matrix, and the addition of four functions to each component of u, each either purely real or purely imaginary. Compared to the evaluation of nonlinear terms, the computational overhead is negligible.

Method for the magnetic field
Consider the induction equation (4) time-stepped without ∇ · B = 0 implicitly enforced.
Evolution of ψ = ∇ · B is then governed by the divergence of (4), In addition to the boundary conditions derived in [18], the condition ψ = 0 must be satisfied on the boundary, so that (4) has the appropriate number of boundary conditions and to avoid introduction of divergence into the domain.
To prevent accumulation of divergence from artificial internal sources, i.e. discretisation error, it is commonplace to introduce an artificial pressure Π [22]. The discretised system is then as in (14) where one reads B for u and Π for p. The boundary condition for Π is any choice such that, when one computes Π for a given B q , it is found to be constant for ∇ · B q = 0 [23]. The choice ∂ r Π = 0 is applied here. When comparing with the problem for the velocity, here the difficulty is not the coupling of the boundary condition for B with Π, but between the components of B at an insulating boundary.
Here the system is split as in (15) for the 'bulk' solution, with approximate boundary conditionB = B q on {r i , r o }. This is then corrected precisely via the influence matrix requiring only the simple auxiliary system  Again, the auxiliary functions B ′ j (r), the matrix A and its inverse may be precomputed, and the boundary conditions for B ′ have been chosen so that B ′ ± are purely real, B ′ z is purely imaginary, and A is real. At the end of the timestep, the solution is solenoidal and satisfies the boundary conditions to machine precision.

B. Implementation notes and parallelization
The Taylor

C. Numerical validation
The code was validated against several published linear stability results, as well as threedimensional nonlinear simulations of the coupled induction and Navier-Stokes equations.
We tested the inductionless limit P m = 0 and finite P m, obtaining excellent agreement in all cases.

Linear stability of Couette flow subject to magnetic fields
Linear instabilities were detected in the calculations by monitoring the kinetic energy of the deviation from circular Couette flow after introduction of a small disturbance. In the linear regime we write where λ = σ + iω is a complex number; the imaginary part ω is the oscillation frequency and the real part σ the growth rate of the dominant perturbation. The latter is readily extracted from the relationship log(E) ∼ 2σt.
We first reproduced the classical results of Roberts [24], who considered the inductionless  Finally, we compared results of the axisymmetric HMRI (helical field B = B 0 (e z + γe φ )) obtained with the spectral code of Hollerbach [27]. A typical diagnostic quantity is the torque at the cylinders The laminar flow torque will be used as a scale, so that the dimensionless ratio G/G laminar Because of the symmetries (see §II B), two different Hopf-bifurcation scenarios are possible [20]. In the first one, the z-reflection symmetry is broken and depending on the initial conditions either upward traveling waves (with k > 0 modes) or downward traveling waves (with k < 0 modes) may be observed. In the second scenario, the z-reflection symmetry is preserved and a standing wave emerges. This is a combination of upward and downward traveling waves for which positive and negative k modes are in phase and have exactly the same amplitude. In both scenarios waves rotate in the azimuthal direction.
We found that at Re = 1480 the circular Couette flow (6) becomes unstable at Ha c = 107.
The emerging pattern is a standing wave (SW) with dominant mode (k, m) = (±8, 1), so that 8 pairs of vortices fit in the domain. Figure 2b shows the square of the amplitude of the complex Fourier coefficient A 8,1 for increasing Ha. As expected in a Hopf bifurcation,  axial wavenumber (black curve on the Fig. 3a), so at the same parameter value states with different wavenumber and torque can be realised depending on the initial conditions. Further increasing Ha the instability is gradually damped until it disappears at Ha ≈ 175. Over the whole Hartmann range the torque never exceeds 1% of the laminar flow (see figure 3a), indicating very weak transport of angular momentum.

V. ONSET OF SPATIO-TEMPORAL CHAOS
At Re = 2960 the Hopf bifurcation occurs at larger Ha c = 120 and the emerging SW remains stable until Ha = 160. At this point a catastrophic transition to spatio-temporal chaos is observed: the vortex structure is damaged and the up-down symmetry is broken (Fig. 4a). Between Ha = 130 and 160 there is a hysteresis region in which both SW and the defect structure are locally stable (see Fig. 3a). In this Ha-range, if the initial condition is a SW from another run with slightly different Ha, this remains stable. However, by starting for example from a disturbed Couette profile the flow evolves directly to defects.
This catastrophic transition is suggestive of a subcritical bifurcation. We investigated this hypothesis by computing the unstable branch separating defects and SW. For this purpose we combined time-stepping with a bisection strategy as follows. If the SW is slightly disturbed, then the flow will rapidly converge to the SW because it is locally stable. The same applies to defects. For intermediate initial conditions the flow will take a long time before asymptotically reaching either the SW or the defects. Such initial conditions were generated here by performing a linear combination between two selected flow snapshots of SW and defects. This combination was parametrised with a variable β, for which β = 0 is SW and β = 1 is defects. By refining the bisection procedure in β the resulting initial conditions get successively closer to the manifold (or edge) delimiting the two basin boundaries. The edge is comprised of those initial conditions that tend neither to defects nor to SW, and the attractor in this manifold is referred to as an edge-state [29]. Figure 3b shows that as initial conditions are taken closer to the edge as the result of the bisection procedure, the dynamics approach a periodic oscillation. Thus the edge-state is a relative periodic orbit or modulated wave, which confirms that the SW becomes unstable at a subcritical Hopf bifurcation. The instability consists of a long-wave modulation of the axially periodic SW pattern, whose spiral structure breaks the reflection symmetry (see figure 5). Despite the simple temporal behaviour, the spatial structure is complicated and the periodic edge-state can be seen as a precursor to defects. Because of the computing cost associated with the bisection algorithm, we only computed the edge-state at Ha = 155 and 140, with similar results.

VI. TURBULENT FLOW
When the Reynolds number is further increased the defects grow gradually into turbulence. For Re = 4000 vortices become small at the inner cylinder and remain quite large at the outer cylinder (Fig. 4b). At Re = 9333 this tendency develops into dynamically very fast small vortices at the inner and slow large vortices at the outer cylinder (Fig. 6a). There is no preferred direction in the system; vortices can travel up or down, both at the inner and outer cylinders.
The qualitative difference between standing wave, defects and turbulent flow is apparent in time series of the radial velocity v r taken at the mid-gap between the cylinders (r, φ, z) = (1.5, 0, 0). Figure 7a shows that the radial velocity of the standing wave oscillates around zero at a certain frequency. The edge state adds a long-wave modulation to the main one ( Fig. 7b). When defects appear in the SW flow structure (Fig. 7c), there are several oscillation frequencies in the time series. Turbulent flow has a completely disorganised pulsating velocity field (Fig. 7d).
The transfer-rate of angular momentum between the cylinders is important for accretiondisc modeling. We checked the torque scaling for increasing Re numbers (see Fig. 6b).
The parameters Re and Ha were chosen according to the results of the linear analysis of the AMRI [8]: we used values of Re and Ha across the line of maximum growth rates of disturbance energy for given Re (see Fig. 1b). The dimensionless torque increases according to the scaling law which is surprisingly low compared to hydrodynamic experiments in the Rayleigh unstable regime [30].

VII. DISCUSSION
We showed that the AMRI in Taylor-Couette flow manifests itself as a wave rotating in the azimuthal direction and standing in the axial direction, thereby preserving the reflection symmetry in the latter. In order to compare to experimental observations [16] we computed the angular drift frequency of the wave. This is shown in figure 8  cylinder frequency (dashed line) and slows down as the Hartmann number increases, which is in qualitative agreement with the experimental data. Note however that in the experiment two frequencies are simultaneously measured, corresponding to the up-and down-traveling spiral waves, respectively. Although in the standing wave the two frequencies are identical, in the experiment the asymmetric wiring creates B r and B z components of magnetic field which break the reflectional symmetry. As a result, up and down spirals travel with different frequencies, similar to co-rotating Taylor-Couette flow in which the reflection symmetry is broken by an imposed axial flow [31]. Another difference is that in the experiment the flow becomes unstable at lower Ha, which can be explained by the different boundary conditions in the experiment from our simulations. In the experiment copper cylinders are used, so perfectly conducting walls would be a closer boundary condition for the magnetic field.
More significantly, in the experiments the cylinders are of finite length, so to reproduce their results exactly a no-slip condition on end-plates should be used. We have applied periodic boundary conditions in the axial direction, which more accurately model the accretion disc problem and allow us to compute high Reynolds number flows more efficiently.
As Re increases a catastrophic transition to spatio-temporal chaos occurs directly from the SW. In a range of parameters SW and chaos are both locally stable and can be realised depending on the initial conditions. We have shown that the first step in this transition process is a subcritical Hopf bifurcation giving rise to an unstable relative periodic orbit, which has been computed using an analogue of the edge-tracking algorithm introduced by Skufca et al. [29] in shear flows. This unstable relative periodic orbit consists of a longwave modulation of the axially periodic pattern of the standing wave and destroys the homogeneity of the vortical pattern. It can thus be seen as a temporally simple defect precursor of the ensuing spatio-temporal chaos. Because of the computational cost we could not track further instabilities on the unstable branch, which we speculate result in chaotic flow before the dynamics stabilize at a turning point (Ha = 130 at Re = 2960) and defects can be computed simply by time-stepping.
We believe that such long-wave instabilities are ubiquitous in fluid flows. In linearly stable shear flows, such instabilities of traveling waves were found to be responsible for spatial localisation [32,33]. In fact, in pipe flow the ensuing localised solutions, which are also relative periodic orbits, suffer a bifurcation cascade leading to chaos [34]. One difference is that in shear flows the traveling waves are disconnected from laminar flow, where the standing wave of the HMRI is connected to the circular Couette flow.
Our simulations were performed with a powerful spectral DNS method, which we have developed and validated with published results to excellent agreement. The method allowed us to compute flows up to Re = 2 · 10 4 . As Re increases, defects accumulate and the flow becomes gradually turbulent. Although we found that the AMRI exhibits a weak scaling of angular momentum transport, with G/G laminar ∝ Re 1. 16 , we expect that larger magnetic Prandtl numbers P m (realistic for accretion discs) may result in a stronger scaling.
Astrophysically important issues such as the precise angular momentum transport scalings obtained for different values of P m, and also for different choices of imposed field (e.g. SMRI, HMRI, AMRI) will be the subject of future investigations.

ACKNOWLEDGMENTS
Support from the Deutsche Forschungsgemeinschaft (grant number AV 120/1-1) and the Jülich Supercomputing Centre is acknowledged.