MRI driven wind cycles in local disc models

Magnetorotationally driven wind cycles in local disc models


Jets, from the protostellar to the AGN context, have been extensively studied but their connection to the turbulent dynamics of the underlying accretion disc is poorly understood. Following a similar approach to Lesur et al. (2013), we examine the role of the magnetorotational instability (MRI) in the production and acceleration of outflows from discs. Via a suite of one-dimensional shearing-box simulations of stratified discs we show that magneto-centrifugal winds exhibit cyclic activity with a period of , a few times the orbital period. The cycle seems to be more vigorous for strong vertical field; it is robust to the variation of relevant parameters and independent of numerical details. The convergence of these solutions (in particular the mass loss rate) with vertical box size is also studied. By considering a sequence of magnetohydrostatic equilibria and their stability, the periodic activity may be understood as the succession of the following phases: (a) a dominant MRI channel mode, (b) strong magnetic field generation, (c) consequent wind launching, and ultimately (d) vertical expulsion of the excess magnetic field by the expanding and accelerating gas associated with the wind. We discuss potential connections between this behaviour and observed time-variability in disk-jet systems.

accretion discs – MHD winds – instabilities

1 Introduction

Observations of accretion discs around T-Tauri stars, in X-ray binaries, or on the larger scales of Active Galactic Nuclei (AGN), have revealed the existence of supersonic outflows and collimated jets (Bally et al., 2007; Fender, 2006; Schwartz, 2007). Sequences of geyser-like eruptions occurring at quasi-regular intervals have been observed in protostellar systems (Zinnecker et al., 1998; Hartigan & Morse, 2007), while relativistic jets emitted by AGN also exhibit rich time-variability, characterized by knotty structures revealed in particular by the CHANDRA X-ray imaging of M87 (Harris et al., 2009). From a theoretical perspective, these powerful ejections (steady or not) are thought to be launched from accretion discs and possibly related to the magneto-centrifugal effect as proposed by Blandford & Payne (1982). In this model, a large-scale poloidal magnetic field, anchored in the disc and corotating with it, is believed to accelerate the gas. When both centrifugal and gravitational forces are taken into account, particles are accelerated along the field lines when their inclination from the rotation axis is more than 30\degree. Magnetic field lines are able to extract angular momentum from the disc and therefore might take part in the accretion process. Detailed models of magnetically driven jets (Lovelace et al., 1991; Pelletier & Pudritz, 1992) and simulations of magneto-centrifugal ejection (Romanova et al., 1997; Ouyed & Pudritz, 1999; Krasnopolsky et al., 1999) have been produced. However, they usually assume that the global magnetic configuration is supported by the accretion disc, neglecting the back-reaction of the disc dynamics on the large-scale magnetic field and vice versa. Also, most of them treat the disc as a boundary surface that loads mass on to the magnetic field lines at a given rate, but in fact the conditions for gas to escape from the disc are still not entirely understood.

More recent work has investigated the effect of the vertical structure of the disc on the outflow. In the case of a strong poloidal field that threads the disc with a gas to magnetic pressure ratio at the midplane, it has been shown that significant thermal assistance is required to overcome the potential barrier between the surface of the disc and the first critical point, namely the slow magnetosonic point (Ogilvie, 1997; Ogilvie & Livio, 1998). In a symmetric disc, this barrier is due to the bending of the field lines close to the midplane. Ogilvie & Livio (2001) showed that, in the absence of an additional source of energy, the mass-loss rate decreases very sharply with increasing field strength and is maximized at an inclination of approximately . Ogilvie (2012) demonstrated that quasi-steady outflows in this regime can be simulated in a local 1D shearing-box model as this contains most of the relevant physical elements.

An interesting regime is the case of weaker fields for which at the midplane. The situation is then altered because the disc is subject to the magnetorotational instability (MRI; Velikhov, 1959; Chandrasekhar, 1960; Balbus & Hawley, 1991). The interaction between the MRI dynamics and the wind outflow has been studied in the last few years but still remains debated. Simulations of 3D MRI turbulence in a shearing box, with a very weak net vertical magnetic field ( from to ) have shown that unsteady outflows can be produced, starting about two scale heights above the midplane (Suzuki & Inutsuka, 2009; Fromang et al., 2013; Bai & Stone, 2013a). The occurrence of transient and local outflows appears to be a consequence of the nonlinear development of the MRI in stratified discs, even if the connection between local outflows and the observed large-scale jets remains unclear. Lesur et al. (2013) showed in a 1D model that a wind can be produced even when starting with a purely vertical field, as the field lines are bent during the growth of MRI channel modes. Magnetic pressure gradients enhanced by MRI modes are thought to play an important role in the launching process as they help to overcome the potential barrier, while far from the disc the magneto-centrifugal effect accelerates the gas.

Although steady or quasi-steady outflows seem to appear in the MRI-unstable weak-field regime, their properties depend strongly on the numerical boundary conditions of the simulation. In particular, the mass-loss rate decreases with the height of the box (Moll, 2012; Lesur et al., 2013; Fromang et al., 2013). Simulations that take into account the disk’s vertical structure and its symmetry also fail to obtain steady jets with velocities exceeding the fast magnetosonic speed. The fact that these solutions depend on the boundary conditions applied at the top of the shearing box begs the question of their connection with the external medium, and even casts doubt on their existence in reality. Other important issues, which are not addressed in the present paper, are the sensitivity of MRI-driven outflows to non-ideal effects such as Ohmic resistivity, the Hall effect and ambipolar diffusion (Bai & Stone, 2013b; Gressel et al., 2015) and the problem of the radial transport of poloidal magnetic flux, in which MRI turbulence is thought to play a crucial role (Guilet & Ogilvie, 2013).

In this paper, we aim to extend previous numerical studies showing that the MRI can excite unsteady winds in the disc atmosphere and quasi-periodic outbursts (Lesur et al., 2013; Fromang et al., 2013), by simulating a regime closer to the marginal stability boundary of the MRI () and investigating in more detail the nonlinear interaction between the instability and these outflows. We show that outburst cycles are a robust feature in 1D vertically stratified simulations and depend on the strength of the imposed vertical field. We present an interpretation of the phenomena which treats the cycle as a sequence of quasi-steady MHD equilibria (which we compute), with the engine driving the system from state to state being either MRI channel modes, or the loss of magnetic flux via the ensuing outflow. The cycle can then be broken down into the following successive phases: (a) the emergence of a dominant MRI channel mode, (b) the generation of strong horizontal magnetic field by the channel, (c) suppression of the MRI by the field and the concurrent launching of a wind, and (d) vertical expulsion of the excess magnetic field by the wind and a return to the initial laminar weak field state. This behaviour could potentially be connected to the variability both in accretion disk and jet emission, but is probably limited to regions of the disk where the MRI is not far from marginality, i.e. the magnetic flux or non-ideal MHD effects are locally strong.

The structure of this paper is as follows. In section 2, we describe the equations and justify the use of a one-dimensional local model. In section 3, we explain the choice of our numerical setup and in particular the tall boxes used for this study. In section 4, we present the phenomenology of time-periodic solutions found in the MRI unstable regime and show that they exist in a large range of parameter values. In section 5, we propose a mechanism that might explain the cyclic process. We calculate the MRI stability of magnetohydrostatic disc equilibria, and connect the results with the nonlinear dynamics of periodic outflows. We finally discuss in section 6 the relevance of our model for accretion discs and the possible implications of our cyclic solutions for astrophysical applications.

2 Model and equations

We use the Cartesian local shearing-sheet description of differentially rotating flows (Goldreich & Lynden-Bell, 1965), whereby the radial, azimuthal and vertical directions are denoted respectively , and . The axisymmetric differential rotation is approximated locally by a uniform rotation plus a linear shear flow with for a Keplerian disc. We assume that perturbations are axisymmetric, which means that they do not depend on in the local framework. As the local approximation has a translational symmetry in , it is possible to study magnetohydrodynamic (MHD) solutions that depend only on . This assumption does not reproduce the complex physics occurring in a realistic disc but might capture the essential mechanisms operating in the wind-launching process. We also note that some 3D studies have found that the dynamics reduces naturally to a laminar 1D form (Bai, 2013).

For simplicity we also assume that the gas is isothermal so that its pressure and density are related by with the uniform sound speed. The system of equations governing the evolution of the fluid density , velocity and magnetic field , in the 1D approximation, is (Ogilvie, 2012)


where and are respectively the kinematic viscosity and magnetic diffusivity. It is possible to define as the standard hydrostatic scaleheight of the disc. This quantity will be considered as our unit of length, while is our unit of time. The source term is an artificial mass injection that mimics the mass replenishment that could occur in a realistic disc when curvature effects and radial transport are taken into account. To obtain steady or periodic solutions of equations (1)–(7), this term has to compensate the loss of mass being ejected by MHD winds. The way mass is replenished is explained in more detail in section 3.2 and discussed in section 4.4.1.

The gravity term reduces to in the standard local approximation. However, this approximation is only appropriate when the vertical scales of interest are small compared to the disc radius. For the reasons explained in section 3.3, wind solutions are better studied in the limit with comparable to the radial extent of the disc. It is then necessary to introduce an additional length scale , which can be regarded as the radius of the reference orbit on which the local model is constructed, i.e. the distance of the origin of the shearing box from the central object. For the vertical component of the central gravity starts to decrease and tends to zero further away. Assuming that the disc can be defined by its local thickness or by its dimensionless aspect ratio


the radius can be expressed in our unit system as . The modified vertical gravity taking into account the finite distance from the central object is then


where . Similar corrections to the vertical gravity have been used by previus authors when modelling accretion discs (Matsuzaki et al., 1997) or galaxies (Kuijken & Gilmore, 1989). For , we retrieve the standard gravity of the local thin disc approximation. The typical for an astrophyscial protoplanetary disc, inferred from recent observations of T-Tauri stars, is between 0.03 and 0.2 (Andrews et al., 2009; Gräfe et al., 2013). Note that this gravity law is not fully consistent with the local approximation of Keplerian flows as the radial gravity force also changes with altitude. However, in order to simplify the problem and still use the one-dimensional model, we will assume that this variation of the radial gravity is unimportant for the dynamical flow we are studying.

3 Numerical setup and parameters

3.1 Numerical Code

We used the shearing-box approach to simulate locally the one-dimensional disc atmosphere. As the flow is compressible and possibly gives birth to a variety of shocks, we used the PLUTO code (Mignone et al., 2007), a finite-volume method with a Godunov scheme that integrates equations (1)–(7) in their conservative form. The fluxes are computed with the HLLD Riemann solver but we checked that solutions are independent of the choice of solver. The time-integration is achieved by a Runge–Kutta method of third order. We emphasize that 1D simulations can be time-consuming, especially in regimes or locations where the typical Alfvén speed is large compared to .

3.2 Boundary conditions and mass replenishment

As in Ogilvie (2012) and Lesur et al. (2013), we impose symmetry with respect to the disc midplane. Simulations are then restricted to the upper half of the disc, as the lower half is obtained by taking , , , , and . This symmetry implies that the poloidal magnetic field is vertical in the midplane, i.e . With such symmetry, we guarantee that field lines bend in the disc midplane, which makes the disc structure more stable.

For the upper boundary , we use standard outflow boundary conditions for the velocity field but impose hydrostatic balance in the ghost cells for the density, in a manner similar to Simon et al. (2011). In this way, we reduce significantly the excitation of artificial waves near the boundary. As explained in more detail in section 4.4.2, various boundary conditions for the magnetic field have been used. A first possible way (Type I boundary) is to fix (and therefore the inclination of the poloidal field) and impose zero gradient for as employed by Lesur et al. (2013). Another is to make a linear extrapolation of and into the ghost cells (Type II boundary). Unless it is explicitly mentioned, all simulations have been run with the first type of boundary condition.

In order to reproduce the mass replenishment due to the radial redistribution of material, we inject mass near the midplane at each numerical time step. The source term in the mass conservation equation (1) is similar to the one prescribed by Lesur et al. (2013),


where is the mass injection rate (into the upper half of the disc) and is a free parameter that corresponds to the altitude below which most mass is replenished. There are two different ways to replenish the mass. The first one is to maintain a constant disc surface density in time. Then, by computing the amount of mass that leaves the box, it is straightforward how much must be injected each time step. Perhaps a more physically relevant method is to inject mass at a constant rate . This approach is probably more realistic but suffers from the fact that the mass contained in the computational domain can vary in time (periodically if the solution of interest is a cycle). We use preferentially the first method and will see in section 4.4.1 how results are affected by choosing the other method. Note that mass is injected without any velocity, although it leaves the box with a non-zero velocity, which implies that the total momentum in the box is not conserved. In fact the loss of horizontal momentum from the box (which can also occur through magnetic stresses acting at the upper boundary) drives a mean horizontal flow, which can be interpreted as an accretion flow together with a small departure from Keplerian rotation.

3.3 Box size and resolution

As pointed out by Fromang et al. (2013) and Lesur et al. (2013), wind properties and mass-loss rates obtained from shearing-box simulations depend on the box size and the upper boundary conditions. This might be a consequence of the fact that critical points, such as the Alfvén point or the fast magnetosonic point, often lie outside the computational domain [definitions of these points are given by Ogilvie (2012)]. Crossing the Alfvén point for is possible by having a tall box with . However, crossing the fast magnetosonic point is more delicate. Previous numerical studies in 1D or 2D models have failed to compute realistic steady solutions that cross this critical point (Lesur et al., 2013), unless the temperature of the atmosphere is strongly enhanced (Ferreira & Casse, 2004) or the vertical symmetry of the disc is broken (Vlahakis et al., 2000). In our case, we did not find any solutions that pass through the fast point. Therefore the boundary conditions and the box size may still influence the numerical solutions, although we found surprisingly that several wind properties seem to converge when a large vertical height is reached (see section 4.4.2).

Another drawback inherent to previous shearing-box simulations is that increases linearly with in the upper disc atmosphere, corresponding to a gravitational potential well of limitless depth; this might lead to solutions that are artificial or dependent on boundary conditions. In particular it is possible to show that any steady solutions with as cannot exist in our setup (see appendix A). This suggests that the behaviour of solutions at high altitude might be better assessed by using the modified gravity of equation (9) provided that .

We then chose a suitable box size that satisfies the different conditions discussed before (except the crossing of the fast point) for the typical range of parameters we used. Note that the main limitation of the shearing box model when considering such large is that magnetic field lines in the upper atmosphere are anchored to radii far from . In reality, these radii will rotate at different periods from those implied by the uniform shear of the local approximation, suggesting that our model is unlikely to be consistent with the global disc equilibrium. However, given the robustness of the phenomena we uncover, we doubt that this inconsistency impacts on our results unduly. A second problem is that the MHD approximation may fail in such tall boxes as the mean free path in the upper atmosphere may be too long. In actual fact this is rarely the case since outflows supply enough mass in this region: the typical density measured at the upper boundary never drops below of the midplane density.

Figure 1: Space-time diagrams of a periodic ideal MHD wind solution computed for , and . Top, middle and bottom colourmaps show respectively the mass flux , radial () and azimuthal () magnetic field. The black/dashed line at corresponds to the separation between the refined and the coarse grid patches. It is also used conveniently to dissociate the ‘lower’ atmosphere (, including the midplane) from the ‘upper’ atmosphere ().

Our numerical grid is a combination of two uniform patches with in the domain and for . The reason for using a finer grid in the midplane region is that the disc scaleheight can reach when the material is subject to strong magnetic compression. In addition, small-scale structures like pressure waves arise more naturally near the midplane. For comparison, we also used a uniform grid with resolution and checked that the solutions do not differ significantly.

3.4 Physical parameters

In all simulations (except the one with constant shown in Fig. 5 & 6), we fixed the surface density to be equivalent to the column mass of a hydrostatic disc with midplane density , which sets our unit of mass. In the one-dimensional approximation, the vertical magnetic field is independent of and . A suitable dimensionless parameter for this problem is then


We focus in this paper on the dynamics in the range , which translates to an average plasma , defined as:


which remains close to the midplane if the disc is not significantly compressed. We expect the MRI to be active in this regime, giving rise to significant outflows.

We assumed ideal MHD (), except in section 4.4.4 where the role of ohmic diffusion is studied. We introduced a very tiny kinematic viscosity in the midplane, . Without any substantial explicit viscosity near the upper boundary, strong numerical instabilities appear, in particular for , for which outflows can be 80 times faster than the sound speed. In order to avoid these instabilities, we forced to be inversely proportional to density at large . Actually this prescription might be physically relevant, at least in the isothermal case. Indeed, as long as the plasma remains collisional, can be estimated as the product of the thermal velocity and the mean free path, which scales as (Maxwell, 1866). Whether this prescription is realistic or not, we checked that viscosity plays a minor role in the results. Viscous forces in the upper atmosphere remain very small compared to inertial forces, while changing the value of at the base does not strongly modify the solutions, as long as numerical instabilities are damped.

4 Time-periodic wind solutions

4.1 Initial conditions and capturing a cycle

The set of initial conditions that converge towards a periodic solution is not necessarily easy to guess. These solutions do not appear spontaneously in small boxes and their study requires the use of large . However, starting from an MRI-unstable hydrostatic equilibrium in very tall boxes with small is numerically unachievable because of excessively high Alfvén speeds. We therefore began with a small box of size and . The initial condition consists of a hydrostatic disc embedded in a vertical field (originally we chose ). Small perturbations were introduced and were amplified by the MRI instability. The nonlinear state evolved rapidly into a steady outflow, as described by Lesur et al. (2013). We then progressively continued this steady solution by increasing the box size and found that it bifurcates to a periodic outflow at . As shown later in section 4.4, cyclic wind solutions are fairly robust to any change in the numerical or physical parameters and we believe they have a physical reason to exist. Unlike the steady outflow obtained for , the properties of these cyclic solutions seem to converge as the box size is increased (see section 4.4.2).

Another method to obtain periodic outflows is to start from the perturbed hydrostatic equilibrium in a very tall box by taking . In that case, the density drops moderately at large altitude so that the Alfvén speed (and therefore the computational time) do not diverge too rapidly with . By means of parameter continuation, is finally decreased to obtain more realistic solutions.

4.2 Phenomenology of wind cycles

We present in this paragraph an example of a cyclic solution obtained for (), and . Figure 1 shows the spatial variation and temporal evolution of the magnetic components and (normalized in the same way as eq. (11)) and vertical mass flux during 70 . The period of the cycle is roughly equal to , i.e. two orbital periods. The dashed line at denotes the separation between the coarse grid and the fine grid. Actually this virtual line remains very close to the sonic line during the whole simulation. This is true for a uniform grid as well, and thus is not an artefact of splitting the box into two domains. It is convenient to use this particular altitude as a separation between the ‘lower atmosphere’, dominated by the disc dynamics, and the ‘upper atmosphere’ dominated by the wind flow. Another important altitude, located below , is where the toroidal current and magnetic pressure gradient change sign (going from positive to negative).

The first colourmap (Figure 1, top) shows that the cycle is made of two phases. One is quiescent with a very weak outflow, , occurring between and near the midplane. The second consists of a violent ejection of gas starting around in the lower atmosphere, with a mass flux measured at the upper boundary a few later. The second colourmap shows that, during the cycle, the amplitude of the radial field varies quite strongly in the lower atmosphere ( varies between 0.3 and 1.4) whereas it remains confined to smaller values in the upper atmosphere. The toroidal field varies at the same frequency with a small positive time lag, comparable to .

Figure 2: Evolution of maximal inclination , maximal normalized toroidal field , midplane density , density at , effective disc scaleheight and Alfvén point altitude during two periods of the cycle of Fig. 1. The left labels correspond to the blue/solid curves whereas the right labels are for the red/dashed curves.

The outflow is initiated near the midplane (about one scale height above it) when is maximum and starts to decrease. The inclination of the poloidal field at this time and for is around 60\degree with respect to the vertical axis, allowing gas to be accelerated along field lines by the magneto-centrifugal effect. This inclination is smaller in the upper atmosphere (beyond the Alfvén point), but remains on average larger than 30\degree. A strong toroidal field is produced in the lower atmosphere, well before the sonic point and Alfvén point are reached. Angular momentum is then extracted from the disc by the toroidal field, allowing the flow to accrete (we checked that near the midplane). In the upper atmosphere and beyond the Alfvén point, its amplitude starts decreasing slowly with altitude, indicating that angular momentum is on the contrary transferred to the accelerating gas. Most of the phenomenology described here is reminiscent of the acceleration mechanism proposed by Blandford & Payne (1982). Note that the model is symmetric under reflection in , so that another solution must exist in which the field lines bend the opposite way and angular momentum is added to the disc, resulting in a radial outflow.

To analyse the cycle in more detail, we plotted in Fig. 2 the evolution of several disc and outflow properties during two periods of the cycle. In particular, we show two dimensionless quantities that will be widely used in the next sections. They are the maximum poloidal inclination in


and the maximum normalized toroidal magnetic field along the same axis,


At , the magnetic variables and are close to their minimum values. The disc is in a relatively expanded state as its midplane density is minimum while its effective scaleheight , defined by , is maximum (close to 0.5 ). The Alfvén point is located far above the disc, in the upper atmosphere, indicating that the wind is rather weak. Between and , the radial and toroidal components of the magnetic field off the midplane are amplified. The magnetic pressure resulting from this amplification rapidly compresses the disc material; its midplane density is multiplied by a factor 3 while its effective scaleheight is divided by a similar factor. At reaches a maximum while the poloidal field lines are strongly inclined, with an angle of 70\degree. The midplane density and remain constant for a while, while the density near keeps decreasing, suggesting that the gas in the disc is being compressed by magnetic pressure. However a fraction of the material in the disc atmosphere, located at the altitude where the current changes its sign, is expelled outward. Fig. 2 shows that the Alfvén point drops suddenly and approaches the location of the lower atmosphere, indicating that a wind is developing in the lower atmosphere. Note that the vertical velocity (as inferred from Fig. 1) is rapidly oscillating during this phase and before the emergence of the jet. Magnetosonic waves propagating upward seem to be excited and are probably associated with the rapid compression of the disc as its height is oscillating with the same frequency.

Figure 3: Left: Projection of the cyclic dynamics in a plane . Two periods are shown here, except for for which only one period has been obtained. Solutions have been computed for a fixed , and but with different ( varies from 2.5 to 16). The red/dashed curve corresponds to the largest and the purple one to the smallest . Around , the solution bifurcates to a steady solution (fixed point) which is reduced to a point in this representation.
Figure 4: Left: Evolution of mean (blue) and maximum (red) mass-loss rate as a function of ( varies from 2.5 to ). The circle markers represent periodic solutions whereas diamond markers stand for steady solutions. Right: Cycle frequency as a function of . The top and bottom figures correspond respectively to () and ().

Between and , the wind accelerates rapidly and turns into a powerful outflow, manifested by the sudden increase of the density at . The jet finally ends when the mass in the lower atmosphere has been completely depleted. During the ejection phase, the amplitude of the magnetic field decreases and the poloidal field lines tend to their initial configuration with a small angle of . As the magnetic compression on the disc decreases, its typical height increases and comes back to its initial value. At , the magnetic field is amplified again, provoking a second outflow and the cycle repeats for ever.

Our simulation shows clearly that the disc and its periodic outflows are driven by the recurrent dynamics of the magnetic field. Off-midplane magnetic pressure is responsible for disc compression and assists the launching of a wind near the midplane, whereas the magneto-centrifugal effect accelerates the wind further away. However, the reasons for which the magnetic field evolves into such different configurations remain unclear at this stage and will be investigated in section 5.

4.3 Dependence of solutions on the vertical field

We examined how the cycle depends on the strength of the vertical field, which is an important quantity in this problem. Because the cycle of Fig. 1 is stable in the parameter space studied, it is straightforward to perform its continuation in . As is slightly decreased (or increased), it takes generally a certain time before the flow converges to a new periodic solution. Figure 3 shows the cycle projections for six different , assuming and . The results indicate that the amplitude of the cycle decreases as decreases. The solution bifurcates to a steady solution (fixed point) for . The average amplitude of the toroidal and poloidal field around which periodic solutions are orbiting increases as is decreased. This dependence on is discussed in section 5.4. We found also that the cycle seems to disappear for , although we failed to explore larger owing to the increasing computational time in this regime. The state obtained in this highly magnetized disc is non-steady and rather chaotic, which is in contrast with the steady flow obtained by (Ogilvie, 2012) for .

Figure 4 (top left) shows that for , the average mass loss rate does not vary significantly with . In case of , we found that a new branch of cyclic solutions emerges after the bifurcation to a steady solution. Their mass-loss rates are similar to the ones obtained for larger , , We found that the cycle persists for small (high ), but the mass loss rate is significantly weaker . Note that the cyclic solutions obtained for such small are characterized by a very slow convergence towards a periodic solution. The transient state is more chaotic and shows episodic bursts of ejection, each of them having different mass-loss rates. Fig. 4 (top right) reveals also that the cycle period increases significantly when is decreased. It is approximatively multiplied by two when is divided by four. We checked whether these results hold for the critical case . Fig. 4 (bottom) shows that for the maximum and average mass-loss rate do not vary significantly with , as long as the solution has not bifurcated. The mass-loss rate reaches a maximum for while it slowly decreases for larger . As soon as the cycle bifurcates to a steady wind (for ), the mass-loss rate decreases strongly by an order of magnitude. Similarly to , the period increases as is decreased.

4.4 Robustness of the cyclic behaviour

Mass replenishment

In the simulations described in 4.2, we kept a constant surface density by artificially injecting mass near the disc midplane. However, this procedure has one major drawback, as depends on time and the mass injected near the midplane depends on the properties of the outflow at the upper boundary. If a violent ejection occurs, then a strong injection occurs as soon as the mass has left the box. Replenishing mass in this way by imposing a fixed could lead to artificial cycles whose period might increase with the box size. To make sure that the cycle we observed is not artificially maintained by the replenishment method, we first checked that the period does not change significantly as is increased or decreased. We also altered the way mass is replenished by assuming a constant mass replenishment with given by equation (10), but with a constant whose value is guessed from Fig. 1 by averaging the mass-loss rate over one period of the cycle. By using this new replenishment method, we obtained a cycle with very similar properties, shown in Fig. 5 (yellow/dashed curve).

Figure 5: Left: Projection of the cyclic dynamics for in a plane where and are respectively the maximum poloidal inclination and normalized toroidal field over the domain. The red/dashed solution corresponds to the reference cycle presented in Fig. 1 . The yellow/dashed curve corresponds to the case of a constant replenishment rate ( varying periodically), the solid/green curve to a modified boundary condition (type II) while the solid/cyan curve has been computed for a lower (altitude below which mass is injected). Right: time variation of the mass-loss rate for the different cases.
Figure 6: Evolution of the mass-loss rate as a function of time for different mass inflow . The initial condition is the cycle for , and . In the bottom figure, the purple curve shows the evolution of .

One may argue that our replenishment is still artificial. Indeed in a real astrophysical disc, there is no reason why the mass inflow rate would take such a specific value. It can vary with time and radius, sometimes independently of the outflow dynamics. We then ran three simulations with different mass inflow rates , and , starting from the cycle computed for and where initially. Figure 6 shows that for both and , the initial state relaxes toward a new periodic solution with different surface density (averaged in time). In fact, changing is somehow equivalent to changing , as is inversely proportional to . One notable result is that the system self-organizes and adjusts so that . For a more extreme case , bursts of ejections are still obtained in the first ten orbits but ultimately they decay into a weak outflow. The magnetization inferred from increases but saturates around unity. This result is reminiscent of the behaviour observed by Lesur et al. (2013) where a transition is observed from a weakly magnetized state with strong outflow to a more strongly magnetized state with very weak outflow. However, in our case, the disc is disrupted after although the mass has not been entirely depleted. In summary, a quasi-periodic state can be sustained for a sufficiently long time, independently of the history of the mass inflow, provided that this inflow is not too small or reduced to zero.

We also checked that changing the typical height below which mass is injected does not significantly affect the solution provided that is not much larger than . In Fig. 5 we compare two different values, (cyan/solid curve) and (red/dashed curve), and show that differences are small between the two periodic solutions. In the case where , the cycle bifurcates to a steady solution. The reason might be that in such a situation there is no need to extract the gas from the disc to launch a wind, as the mass reservoir in the lower atmosphere is directly replenished. Then the dynamics of the midplane region, which we believe to be responsible for the cyclic dynamics, probably becomes less important. However, having a that differs significantly from the effective scaleheight (typically in our simulations) is clearly not realistic as the mass should be mainly replenished in the bulk of the disc.

Boundary conditions and convergence

To study the effects of boundary conditions on periodic solutions, we first analyse their dependence on box size. As explained in section 3.2, the location of the Alfvén point might play an important role in MHD wind simulations. If this point leaves the box then the poloidal field inclination is no longer determined by the Alfvén point crossing condition but ought to be specified through the upper boundary condition as it is determined by global considerations beyond the scale of the box (Ogilvie, 2012). In our setup, is this critical point crucial for obtaining converged solutions? Also, are the cyclic solutions converged with respect to when the aspect ratio of the disc tends to 0?

Fig. 7 shows some properties of the reference cycle of Fig. 1, such as its average mass-loss rate and the maximum normalized toroidal field , as a function of . We analysed three different sets of parameters. The first case is , for which the Alfvén point lies 80% of the time inside the computational domain and leaves the box at regular intervals. The second is , for which the Alfvén point is inside the domain 100% of the time. Finally the last case corresponds to the limit with an Alfvén point that is inside the box most of the time. For the two first configurations, the average mass-loss rate and , as well as (not shown here), seem to converge to asymptotic values for large . In particular the mass-loss rate appears to hover around 0.026 for the first case and 0.035 for the second one. The convergence appears slightly better when the Alfvén point stays fully inside the domain. Note that the solutions depend strongly on the upper boundary location for , in agreement with previous studies. In the case , for which the potential diverges as , we found that the mass-loss rate starts to converge (but were unable to follow numerically the solution at because of a too strong vertical flow), while keeps increasing linearly at large . This result suggests that periodic solutions are converged when the box size is much larger than and when the Alfvén point remains fully in the computational domain. However, we admit that the convergence is rather slow (one needs to reach high altitudes ) and this might be due to the limitation of the shearing box model discussed in section 3.3 and the fact that the fast magnetosonic point is not reached.

Figure 7: Top: average mass-loss rate (left) and maximum of (right) over one period as a function of the box size for three different sets of parameters. Diamond markers indicate that the solution has bifurcated to a steady state. Bottom: Alfvén point altitude as a function of time in these three configurations.

We checked also that our cycles are not strongly affected by changing the nature of the upper boundary conditions. Figure 5 shows a comparison between the cycle of Fig. 1 computed with Type I boundary conditions (red/dashed curve), for which the poloidal magnetic field inclination is forced to be 0, and the same solution computed with Type II boundary conditions (solid/green curve), for which the horizontal magnetic field is extrapolated away from the domain. This result suggests that the cyclic mechanism is robust to the conditions imposed by the external medium in the upper atmosphere.


To understand how the gravity affects the dynamics of these cycles, we performed several simulations by changing the aspect ratio . Simulations are initialized from the reference cycle of Fig.1 obtained for . We scanned a large range of , going from to . Figure 8 shows that a recurrent dynamics is obtained for all aspect ratios, suggesting that the cycle process is independent of the shape of the gravitational potential at large . In particular for (which correpsonds to and marginally for ) the cycle is only weakly dependent on . However, we note that the amplitude of the magnetic cycle increases as decreases. We obtained also a stronger vertical velocity for and . This is consistent with the fact that a more powerful wind has to be produced to oppose a stronger gravity. A shock propagating upward is also observed when tends to 0, as a result of the strong disc compression. Note that there are two limits for which the cycle deviates from the original one. First the limit . The fact that the solution does not converge in that limit is consistent with the results of Sect. 4.4.2 which shows that convergence requires at large (the case of a uniformly increasing gravity is clearly problematic). Second, the limit which is inconsistent with the thin-disc approximation and probably irrelevant for astrophysical discs. In that case, the vertical gravity decreases sharply over a scale comparable to , which might change the structure of the disc outflow.

Figure 8: Left: Projection of the cyclic dynamics in a plane . Solutions have been computed for a fixed , and but with different . The red/dashed curve corresponds to the extreme and unrealistic case . Right: time variation of the mass-loss rate for different (decreasing from top to bottom).

Ohmic dissipation

The simulations described so far have been done in the ideal limit where the effect of microscopic ohmic diffusion is not taken into account. When non-ideal physics is considered, it is possible to define the magnetic Reynolds number Rm as the typical ratio between the ideal MHD term and the resistive term in the induction equation:


We checked that the cycle dynamics is maintained for a large range of Rm . Protoplanetary discs are generally poorly-ionised but can have magnetic Reynolds number Rm larger than 10, especially in the outer regions. AGN and binary discs have a wide range of resistivity and quite large Rm near their centres (Balbus & Henri, 2008). Therefore we believe that the existence of large-scale periodic outflow solutions in astrophysical discs might not be limited by such dissipative effects, provided that a net magnetic poloidal flux can be maintained for a long time in the disc.

5 Investigation of the cycle mechanism

5.1 An MRI-driven cycle?

We showed in the previous section that the evolution of mass-loss rate and disc structure during the cycles is intimately connected to that of the magnetic field. However, the time variability of and remains unexplained. In order to unveil the nature of the cyclic process, it seems necessary to understand the dynamics of the magnetic field and its interplay with the disc atmosphere.

As mentioned earlier, the MRI is one obvious candidate to explain the growth of the horizontal magnetic field and the launching of a wind. Some clues, detailed further in section 5.3, suggest that the MRI is active during the first phase of the cycle. However its nonlinear evolution and connection to periodic behaviour are not obvious as it evolves in a complicated background. In this section, we explore the interplay between the MRI, the stratification, and the wind in order to understand the cyclic magnetic activity. As part of this, it is first necessary to probe in detail the dynamics of the MRI on a background that is different from the pure hydrostatic equilibrium, taking into account the strong horizontal magnetic fields associated with the wind flow.

5.2 MRI stability of magneto-hydrostatic equilibria

As the cycle period is significantly longer than the typical MRI growth time (), the dynamics of the midplane region and lower atmosphere can be modelled to a first approximation as a succession of quasi-magnetohydrostatic equilibria. The action of the MRI is one key ingredient that might control the disc migration from state to state. We computed a particular class of quasi-steady equilibria, as solutions to suitable boundary value problems, and their stability. Through this study, we aim to identify the conditions for which the MRI is excited and to characterise its behaviour when a background poloidal inclination and a toroidal field are present.

Families of 1D nonlinear MHD equilibria

Nonlinear equilibrium solutions are the local manifestations of the global equilibria calculated by Ogilvie (1997) and connect to the local winds calculated by Ogilvie (2012). They also generalise the weakly nonlinear states calculated by Liverts et al. (2012) at the onset of the MRI. But being steady, they must be distinguished from the exponentially growing channel flow solutions of Latter et al. (2010). The wind equilibria are not trivial to obtain and their existence is still matter of debate. To simplify the problem, we first make the assumption that . This approximation is acceptable for our purposes because the vertical velocity is small compared to the sound speed in the region where the MRI is active (i.e. the disc midplane and the lower atmosphere). We also do not include resistivity and restrict our study to symmetric equilibria, i.e with and . Although equilibria also exist in the asymmetric configuration, we found that they take the shape of magnetically ‘levitating’ discs in cases of and might be strongly unstable to buoyancy-type instabilities.

The steady nonlinear equations and their solutions are briefly presented in Appendix B. The system to solve comprises five first-order differential equations and depends on the same number of governing parameters. Two of them are already set by the symmetry. The three remaining parameters can be chosen as the surface density , the maximum of the normalized toroidal field and the maximum poloidal field inclination The two last parameters are directly associated with the quantities we used to describe the magnetic state of cycles. is always equal to the asymptotic value taken by and quantifies the toroidal magnetic pressure in the upper atmosphere of the disc. Note that although we do not consider any wind in the equilibrium solutions, the existence of a toroidal field and its related torque is inevitably associated with the physical conditions at the outer boundary. In a realistic scenario, they can only be maintained by a significant outflow coming from the disc.

Instability in 1D and role of the toroidal field

Figure 9: Stability of magnetohydrostatic equilibria as a function of and for (top) and (bottom). The colourmaps show the MRI growth rates of the mode. A blue pixel () indicates a stable equilibrium. The white contour lines corresponds to different iso-values of the Alfvénic angular frequency defined in equation (16), for which is obtained numerically. The solid line seems to lie on the marginal stability border of the mode in the case . The yellow dotted line is the marginal stability border calculated analytically in the limit and . Note that growth rates were difficult to obtain in the regime of large and small .

The states computed in the previous subsection constitute exact steady nonlinear solutions to the governing equations. However, the dynamical system, if started at (or near) one of these states, may evolve away because of instability, in particular the MRI but also potentially the magnetic buoyancy instability (Parker, 1966). In the classical analysis, the MRI occurs when magnetic tension is relatively weak, , while the buoyancy instability (in isothermal gas) emerges when the equilibrium is strongly deformed by the poloidal field so that where is the norm of the horizontal field (Hughes et al., 2012). In this subsection, however, we investigate only the possibility of MRI. Previous studies in the vertically stratified box have focused exclusively on the MRI stability of a purely hydrostatic equilibrium with a uniform vertical field. (Gammie & Balbus, 1994; Latter et al., 2010). Our results here generalise this to include basic states exhibiting strong and non-uniform toroidal fields. Some of these results have been computed earlier by Ogilvie (1997, 1998), who examined the stability of discs with bending poloidal magnetic fields. Axisymmetric and non-axisymmetric magnetic instabilities have also been studied in the disc context by Terquem & Papaloizou (1996). The stability problem for such equilibria can be solved by numerical methods which are briefly outlined in appendix C.

The three governing parameters are , and , which together set the background equilibrium state. First, for a given inclination , we constructed a two-dimensional stability map in the plane. For each point in this map, we computed the corresponding equilibria and solved the linear eigenvalue problem. When the disc is unstable, it is then possible to report the maximum growth rate of the instability at that point. Figure 9 shows two of these maps for different inclinations and . Growth rates are computed for the largest-scale unstable mode which corresponds to in our linearised symmetric subspace ( if the anti-symmetric solutions are taken into account). The instability is clearly identified as the MRI, since for growth rates are identical to those obtained by Latter et al. (2010) in the hydrostatic case. When (top figure) and is small, the instability reaches its maximum growth for , while it is suppressed beyond a critical . As increases, the critical below which the largest MRI mode becomes unstable decreases. In other words, for a fixed , there is a critical above which the MRI is suppressed. The same result is obtained when the poloidal background field is inclined with respect to the rotation axis. In the case of , we found however that the MRI stabilizes at a larger critical for small . For both inclinations, we note that the MRI growth rate becomes small, even negligible, when the toroidal field .

We also computed the same map for larger inclinations, up to (not shown here). When the inclination becomes large enough, another instability seems to arise for and small with a larger growth rate than the MRI. These modes are more difficult to interpret but could be the magnetic buoyancy modes as the background density distribution is strongly deformed by the poloidal field (see the density profile in Fig. 16 for ).

MRI criterion

A strong magnetic pressure induced by the toroidal field seems to quench the MRI. This result can be interpreted as follows: as increases, the disc is compressed and its typical height is reduced (with defined in section 4.2). This forces the modes to be of smaller vertical size, and if this characteristic length is less than the Alfvén length the mode is stabilised. In the classical analysis this is equivalent to . By analogy with the unstratified case, let us define the Alfvénic angular frequency as the product of the Alfvén speed measured in the midplane times the typical wavenumber of the unstable mode (Latter et al., 2010). Then


We then expect marginal stability when (in our units). We compute () numerically for each equilibrium and plot its iso-contours on top of the first colourmap of Fig. 9. For , it seems clear that the curve lies close to the stability border of the mode, indicating that is a suitable quantity to infer the stability of the MRI. We checked that (dotted line in top of Fig. 9) or equivalently lies also on the marginal stability of the mode. For larger inclinations (), the curve reproduces the stability border of the mode only for lower . The discrepancy at stronger vertical fields might issue from the fact that in this regime, equilibria display an extended height (see Fig. 16 of appendix C), allowing the MRI mode to be of larger scale and therefore unstable for larger . Another reason might be that the strong positive vertical shear (a result of Ferraro’s isorotation law), combined with the vertical motions of perturbations, can increase locally the shear rate seen by two fluid particles.

We show in Appendix C that an analytical expression for the marginal stability border can be obtained in the limit of large toroidal field . Our calculation shows that this border, drawn in Figure 9 (dotted yellow curve), is given by for the mode in the symmetric subspace. It is straightforward to check that this curve lies along the MRI stability border at large for both inclinations. The theoretical curve fits the numerical stability border for in the case of moderate inclination and for in the case of zero inclination.

Figure 10: Cycles (plain curves) and marginal stability borders of the largest MRI symmetric mode (dashed curves) for different . Red, green and blue colors corresponds respectively to , and . Cycles are computed for and . Note that for , the branch of equilibrium that we follow cannot be continued at larger inclination.

5.3 Nonlinear evolution of the MRI and connection with the cycles

Figure 11: Nonlinear evolution of an MRI eigenmode growing on top of a magnetohydrostatic equilibrium, for and . Top: projection of the dynamics in a plane (left) and in a plane (right) where is the midplane density and is the mass-loss rate. Simulation starts at from the perturbed magnetohydrostatic equilibrium. The green curve is the trajectory in phase space that corresponds to the linear amplification of the MRI eigenmode. The red dashed curve is the periodic solution toward which the dynamics is converging. Bottom: Spacetime diagram showing the evolution of , and . In the last colourmap, the vertical extent is reduced to for visibility reasons. The white line corresponds to the disc scaleheight .

The aim of this subsection is to link the stability analysis of the equilibria studied in the last section with the cycles found in section 4. We showed that for a given cycle, a strong toroidal field can build up in the upper atmosphere. The results of section 5.2 suggest therefore that the MRI is stabilized by such toroidal field at some stage of the cycle. Fig. 10 indicates that for a fixed , the projection of the reference cycle into a plane periodically intersects the MRI stability border. This means that it alternates from an MRI-unstable to an MRI-stable configuration. For a lower close to the bifurcation point, the cycle becomes tangent to the stability curve while right after the bifurcation, for , it is completely embedded in the unstable region. This figure indicates that a connection might exist between the cyclic dynamics and the change in stability of MRI channel modes.

Although the suppression of the MRI during the second phase of the cycles seems to play a central role in their dynamics, it does not explain precisely why the disc returns to its weakly magnetized state. To understand more deeply the link between the MRI stabilization and the decay of the magnetic field, we examined the nonlinear evolution of an MRI mode using PLUTO. We allow a wind to be produced: . The initial condition consists of an equilibrium in the three-dimensional parameter space perturbed in the direction given by its unstable MRI mode whose theoretical growth rate is . We used a box of size , with a uniform grid of 4000 points and (this large ratio is unrealistic but allows us to simulate the linear phase of the MRI in a reasonable amount of time).

The results of the nonlinear simulation are shown in Fig. 11. First, as expected, between and (in units of ), the MRI mode is linearly amplified as and grow exponentially with rate . The projection of the dynamics in a two dimensional plane (top panels) shows that the toroidal component and the mass flux are also growing exponentially during this phase. The magnetic pressure is then enhanced in the lower atmosphere and the disc is compressed, exactly as described for the cycle of section 4. At , the poloidal inclination and reach a maximum. This is followed by a very violent mass ejection and the emergence of shocks produced in the bulk of the disc. At , the magnetic field amplitude decreases rapidly while the dynamics becomes more gentle and less chaotic. Finally the magnetic field is amplified again and the flow seems to converge toward a stable cycle indicated by a dashed red line. The cycle phenomenology looks like the one described in Fig. 8 for .

Figure 12: Magnetic energy budget for (top) and (bottom) as a function of time. Each curve is associated with the contribution of a term in the induction equation. Solid lines represent positive terms (energy sources) whereas dotted lines represent negative terms (energy sinks).
Figure 13: Total magnetic energy budget as a function of time. The plain/purple curve is the vertical Poynting flux (outward-directed, so that it is positive when leaving the domain ) and the green/dashed curve is the integral of in the domain where is the Lorentz force.

The periodic evolution of and the role of the MRI in generating such cyclic dynamics can be understood by examining the energy budget of the poloidal and toroidal fields. It is derived from the induction equation and reads:


The right-hand side terms in the radial field budget are associated respectively with the linear induction, compressible and advective terms. The contributions to the toroidal field are similar except that the linear induction has an additional term , owing to the shear, also known as the -effect. Fig. 12 shows these different energy contributions integrated in the lower atmosphere from to . During the first phase of the cycle, between and , is amplified and mainly gains its energy from the induction term whereas gains its energy from the -effect. This budget is similar to this situation found when the initial MRI eigenmode is linearly amplified (before ) and is reminiscent of previous budgets of MRI unstable flow (Riols et al., 2015). Therefore, this result suggests that the MRI is active during this phase.

On the contrary, between and , when the horizontal components of the magnetic field are maximized, the inductive term associated with the MRI drops sharply and the dominating term becomes . The suppression of the MRI during this phase seems to be directly related to the effect described in section 5.2.3. The magnetic compression reduces the disc scaleheight and forces the MRI modes to be of smaller scale. This results in a stabilizing magnetic tension, in the same manner that the MRI is suppressed when the fluid is strongly magnetized. Quantitatively, we have when the disc reaches its maximum compression and the diagram of Fig. 10 indicates that for and , the MRI is actually suppressed in the range of inclinations considered here. Note however that we cannot extrapolate directly the marginal stability curves to this nonlinear calculation as they have been computed for .

Finally - and this was our main question - why does the disc return to its initial weakly magnetized state (since we showed that the highly compressed disc is MRI-stable)? Actually, the second term in the energy equation is always negative (the advective term being also negative but negligible). Therefore if the MRI is quenched, there is no other source of energy and the magnetic field has to decay. The existence of this negative inductive term, associated with the wind acceleration, clearly breaks the stability of the system and brings it away from its highly compressed state. Note that if the radial field directly loses energy through this term, the toroidal field is affected indirectly through the -effect. In other words, the wind generated by the MRI and the magnetic pressure gradients in the lower atmosphere have a negative feedback on . However, one may wonder how this fraction of magnetic energy is lost physically. Is it transferred to the outflow through Lorentz force, or simply expelled outward by the wind through the expansion of the gas (which differs actually from a pure advective process)? To answer this question, one can re-arrange the induction equation and write the total magnetic energy equation in the form:


Here is the Poynting flux which represents the amount of magnetic energy released by the disc through a given vertical boundary. is the work done by Lorentz force and quantifies the exchange between kinetic and magnetic modes. Figure 13 shows that during the second phase of the cycle, when the MRI is suppressed and the magnetic field starts decreasing, (between and for example), the Poynting flux through the boundary is the dominating term, indicating that magnetic energy is mainly expelled outward. At the same time, the term including the Lorentz force work and the -effect is smaller. Although the can be negative, there is no significant transfer from magnetic to kinetic energy on average. Note that the Lorentz force work is, on the contrary, enhanced during the first phase of the cycle, when the field is amplified by the MRI. By combining the results of Fig. 13 and Fig. 12, it becomes clear that the acceleration and the resulting expansion of the gas in the lower atmosphere (owing to the wind) releases the surplus of magnetic energy accumulated in the disc by the MRI and allows a cycle to exist. The same process occurs in a steady wind, but in that case the magnetic energy is continually produced by the MRI and released by the wind without accumulating.

A similar simulation has been achieved at much lower . An equilibrium has been selected in this regime and we perturbed it along its unstable eignmode. We used it as an initial condition in a box of . Again, a cycle was found in the nonlinear phase, with a period slightly longer than that obtained for . Surprisingly, we found that the average mass-loss rate over one period is , larger than the one obtained for higher . We also found a similar cycle at with a completely different code (developed by G. I. Ogilvie, which uses high-order finite differences) where . These two solutions seem to be different from the branch of solutions described in Fig. 4, which undergoes a bifurcation to a fixed point at . This result shows that powerful cyclic outflow solutions can be found in a large range of and can involve small-scale modes when the vertical field becomes weak.

5.4 Vertical field and condition of existence

We showed numerically that the reference cycle studied in section 4.2 exists above a certain critical . At this critical value, it bifurcates to a steady wind that clearly lies in the MRI-unstable region. Figure 9 shows that the growth rate associated with the largest-scale MRI mode decreases sharply near this value of , independently of the inclination of the poloidal field. In addition, Fig. 4 seems to indicate that the maximum mass-loss rate is obtained for close to 0.3–0.4, which corresponds to the maximum growth rate of the MRI mode according to Fig. 9. These two results together suggest that the mechanism involved in the periodic behaviour and summarized by Fig. 14 is sensitive to the growth rate of the MRI mode that drives the cyclic dynamics. In particular, if the growth rate is too small, the wind has time to develop and release the surplus of magnetic energy produced by the MRI. Therefore the system does not develop a strong toroidal field, nor a powerful jet, but still ejects continuously a small amount of mass.

Although the cycle disappears for in the case of and for a small range of for , other periodic solutions seem to exist at much smaller but are supported by smaller-scale MRI modes (see end of section 5.3). Indeed, the only way to maintain a significant growth rate in that regime is to excite smaller-scale structures. It is then possible to imagine a family of cyclic solutions for which the dominant MRI mode cascades to small scales when is reduced. The minimum that allows a cyclic dynamics will be probably fixed by the dissipative scale. Also, one could expect even more complex 1D structures than explored so far, involving the excitation or interaction of multiple MRI modes, as suggested by the long transient chaotic dynamics obtained at small before converging to cyclic solutions. However, it is important to emphasize that the cycles or more chaotic 1D structures found in the weak field regime (high ) might be highly idealized and probably never seen in the fully 3D turbulent discs. This is discussed in more detail in Sect. 6.2

6 Discussion and astrophysical implications

6.1 Summary of the results

By using one-dimensional isothermal simulations in the local (shearing box) model of an accretion disc, we have shown the existence of a rich set of solutions: cycles, steady states, chaotic magnetically dominated states, transient bursts of ejections. We examined in particular periodic wind solutions that appear when the disc contains a vertical magnetic field with midplane plasma beta . The periods of these cyclic solutions are around and seem to decrease with increasing field strength . They are characterized by a sudden and violent ejection, during which a significant fraction of the mass and angular momentum of the local region of the disc is carried away, followed by a ‘quiescent’ stage. They appear to be robust with respect to the variation of several physical parameters (disc aspect ratio, resistivity, vertical field strength, mass inflow) and computational details (grid, boundary conditions, mass replenishment method). The magnetorotational instability (MRI) appears to play a crucial role in the dynamics of these solutions.

Figure 14 sketches a possible mechanism driving the cyclic dynamics found in numerical simulations and based on the different processes discussed in section 5.3. Starting from a weakly magnetized and compressed disc, the MRI amplifies horizontal magnetic components, resulting in a strong increase of magnetic pressure in the lower atmosphere. An outflow is triggered by magnetic pressure gradients and accelerated by the magnetocentrifugal effect higher in the atmosphere. At the same time, magnetic pressure vertically compresses the disc, making the MRI ineffective. The amplification of the horizontal magnetic field stops and the magnetic flux accumulated in the lower atmosphere is expelled outward by the expansion associated with the wind. The first phase of the cycle is then characterized by an MRI amplification whereas the second phase is dominated by the wind dynamics. If the role of the MRI on the wind-launching process has already been highlighted by many authors, we show here that the wind itself, by acting on the large-scale magnetic field, can indirectly affect the conditions of excitation of the MRI. In particular we showed that the amplification of a toroidal field (mainly by the -effect near the midplane) seems to quench the MRI. This feedback makes possible the existence of regular outbursts in the disc atmosphere.

Figure 14: A sketch of the possible cycle mechanism

6.2 Connections with other MRI-driven wind simulations

In this work, we have been able to study disc wind solutions in a particular geometrical configuration (axisymmetric, no radial dependence) and a magnetic regime where the disc flow is relatively large-scale and cannot be qualified as ‘turbulent’. The disc is marginally unstable to the MRI so that it oscillates between MRI-unstable and stable configurations. One may ask whether these solutions are relevant or not for astrophysical discs.

A first point that needs to be discussed is the connection between our results and the 3D local simulations with very high , a regime that has been explored by several authors and for which the disc is supposed to be fully turbulent. Can we expect similar periodic structures in that regime with significant mass-loss rates? The simulations of Fromang et al. (2013); Suzuki & Inutsuka (2009); Suzuki et al. (2010) suggest that a quasi-periodic outflow can exist when the turbulence is fully developed. Similarly to our results, the variation of the magnetic fields is correlated with the evolution of the vertical mass flux, suggesting that horizontal fields are involved in the cyclic lauching mechanism. However their solutions display much weaker oscillations and mass-loss rates. The orientation of the poloidal field lines also changes periodically ( goes from positive to negative values), which is not observed in our case3. Their simulations also indicate a more chaotic behaviour, which is expected since many active modes exist further from criticality. A first step will be probably to explore 2D axisymmetric simulations in the local or global models with larger to understand whether cyclic solutions exist in this setup and whether they bifurcate to a more complicated dynamics. In particular it would be interesting to see if the cycles develop instabilities similar to those discussed by Lubow & Spruit (1995), recently obtained numerically by Moll (2012) in 2D and Lesur et al. (2013) in 3D. In addition, they might be affected or even destroyed by the non-axisymmetric small scale MRI turbulence. Although the cyclic solutions we computed might be unstable or subdominant in the 3D case with small , it is however not excluded that they can form the backbone of a sustained chaotic or turbulent wind dynamics in that regime.

It is also worth comparing our 1D shearing box calculations with the case of 3D strongly magnetised discs ( comparable to ours), recently simulated by Bai & Stone (2013a) (=) and Salvesen et al. (2016) () in the shearing box. Their simulations show some activity related to the MRI but do not exhibit any cyclic wind motions. It is however possible that their box size (typically for one half disc) is not tall enough to capture such periodic wind solutions. The strongly magnetized simulations are particularly relevant to the so-called “magnetically arrested disks” (MAD, e.g. Tchekhovskoy et al., 2011; McKinney et al., 2012) which are thought to exist around black holes. These discs transport a large amount of magnetic flux that accumulates at small radii. Simulations of such discs show strong episodic outflows that release the excess of magnetic flux at the vicinity of the black hole (Tchekhovskoy et al., 2011). However, it is still unclear how the MRI operates in these objects. It might be quenched as a result of strong magnetic fields. Thus, the cyclic mechanism studied in this work might only be effective during the early phase of the disc, when the magnetic field is building up.

Another issue we addressed in this paper is the dependence of solutions to boundary conditions and in particular the height of the box. Several authors have argued that this effect might result exclusively from the use of the shearing-box model. Actually we showed that one can circumvent this problem by using a modified gravity (which represents the true variation of vertical gravity with height at a fixed cylindrical radius) and ensuring that the height of the box is larger than , which is the radial distance of the centre of the box from the central object. However it remains unclear in what way a sub-fast magnetosonic outflow depends on the boundary conditions and to what extent the external medium can influence the properties of the outflow. Our preliminary work tends to show that they are independent of the nature of the boundary condition (fixed inclination, extrapolation) if the flow remains super-Alfvénic, but more work is needed to understand the role played by the boundaries.

Finally, it is worth discussing the relevance of shearing-box simulations as they clearly do not reproduce all aspects of a real disc. First, they do not deal with the problem of the evolution of the vertical magnetic flux. How is the large-scale poloidal field maintained during the disc lifetime? Some studies suggest that the advection of field lines by the accreting disc can be sufficient to accumulate a strong poloidal magnetic field in their center (Guilet & Ogilvie, 2013), although these studies do not resolve the MHD turbulence inside the disc, using instead prescriptions for turbulent diffusion. Moreover, to simulate periodic solutions, we artificially injected mass near the midplane. However, this mass cannot be computed in a consistent way from the local accretion rate because of the horizontal invariance of our solutions in the shearing box. Global simulations are needed to simulate consistently the accretion-ejection problem and understand how quasi-periodic structures can be obtained in such cases.

6.3 Connections with observations

The idealized setup of the numerical simulations presented in this work (1D local model, isothermal gas) makes it hard to extrapolate directly to astrophysical discs. However, our model might incorporate enough physical ingredients to allow us to understand some apsects of the interplay between the wind and the MRI dynamics. First, there are some observational evidence that real accretion discs (X-ray binaries or AGN) may possess a strong or moderate poloidal field (Fender et al., 2004; Miller et al., 2006; Martí-Vidal et al., 2015), which suggests that the periodic behaviour described in this paper is permitted. Second we think that the cycles and the mechanism described in Fig. 14 might play a role in accretion discs wind dynamics and have a connection with the variability observed in some protostellar or protoplanetary discs, if the conditions of excitation of the MRI are satisfied in such discs. They could also possibly be involved in the quasi-periodic oscillations (QPOs) measured in the light curves of some discs around white dwarfs, neutron stars or black holes. These phenomena generally have a short period that can be similar to the orbital period . For comparison the wind cycles studied here have a period . The jet variability measured in our simulations also seems strong enough to be detected in light curves. However, this preliminary work does not allow any firm conclusion on the possible link between a cyclic MHD mechanism and the observed oscillations. A more exhaustive study, comparing instrumental spectral signatures and synthetic data from simulations will have to be achieved to understand if such connections exist.

Appendix A Existence of steady solutions in the linearised gravity

We show in this appendix that steady wind solutions with the standard linearised gravity do not exist. Steady states are obtained by taking in equations (1)–(7). Assuming that the viscous term is negligible, we can rewrite equation (4) as




is the combined pressure of the gas and the horizontal magnetic field. (The pressure associated with is constant.) For simplicity we assume that (mass is replenished only at ), which implies by mass conservation that . The argument can be generalized to allow for a situation in which all the mass is injected below a certain non-zero height. Integrating equation (20) with respect to from 0 to an arbitrary positive altitude gives


where the subscript refers to values at the midplane. This equation relates the change in the vertical momentum flux to the weight of the column and the difference in the combined pressure. Since and are positive, while is negative, we deduce that


(In fact, for symmetrical solutions, , so .) This equation implies that the vertical velocity is bounded above, which in turn means that is bounded below and cannot tend to at large . In the standard shearing-box model with , it is then straightforward to see that the weight integral diverges to as . This makes it impossible to balance equation (22) in the limit of large . Therefore steady solutions cannot exist in this case.

In the more realistic modified gravity given by equation (9), which has for small but for large , the weight integral is finite in the case that both and tend to non-zero constant values as . Therefore steady solutions can exist, in principle, in this case, although they would have an infinite column density.

Appendix B Nonlinear steady states and numerical methods for their stability

The dynamics of MRI channel modes in a magnetized background can be studied by computing the stability of 1D magneto-hydrostatic equilibria that satisfy the equations (2)–(6) with and . The set of nonlinear equilibria satisfies the following equations:


For a given , as we are studying solutions that are symmetric about the mid-plane, with , , being even in while and are odd, the number of free parameters that describe the whole family is reduced to three. As explained in section 5.2, we can fix the surface density (which is well-defined as long as decays sufficiently rapidly with altitude), the maximum of the toroidal field and the maximum inclination .

Acceptable solutions of the steady equations require that the radial momentum tends to zero at infinity, which implies that reaches a constant value at . It is straightforward to show that the quantity , defined as the maximum of , is directly related to this asymptotic value. The same argument does not hold for since the poloidal field inclination at is not equal to . The main reason for choosing as a parameter instead of is that for a fixed inclination , there is a critical below which physical equilibria (with decreasing density profile) are not permitted. In order to scan the whole parameter space, it is then better to use the maximum inclination (13), located most of the time near the midplane at .
All these solutions have a uniform accretion flow in , which is driven by the magnetic stress acting above the midplane and is directly related to the strength of the toroidal field:


To compute one of these equilibria, we integrate the system of differential equations using a Runge–Kutta integrator of 4th order in space, starting from an estimated midplane density and azimuthal velocity. These two guessed quantities are then refined with a Newton algorithm to match the specified inclination and toroidal field .

Figure 15: Magnetohydrostatic equilibria for . Density (left) and toroidal field (right) as a function of are shown for different . The solid/red curve corresponds to the pure hydrostatic equilibrium.
Figure 16: Magnetohydrostatic equilibria for two different inclinations and . Quantities , , and are shown as a function of (from left to right and top to bottom). The toroidal field at infinity is fixed so that .

Figure 15 shows solutions for and three different . For this inclination both and are identically zero. Increasing clearly compresses the disc by reducing its characteristic vertical scale and by increasing the midplane density. Figure 16 shows a second set of equilibria computed for inclinations of and . Increasing gives rise to a magnetic pressure component localized around which has been induced by the poloidal field. For large (dashed curves), this results in a drastic deformation of the density profile, since a local extremum starts developing at large (this is particularly visible for ). For small (solid curves), and oscillate near the midplane with a wavenumber that scales as . Note that departure from the Keplerian flow increases linearly at large distance as a consequence of field lines co-rotating with the disc and . However the total momentum integrated over remains finite.

To study the stability of these equilibria, we added small perturbations on top of them and linearized the system. Eigenvalues and eigenmodes of the problem are obtained by using a Chebyshev collocation method on a Gauss–Lobatto grid. This results in a matrix eigenvalue problem that can be solved using a QZ method (Golub & Van Loan, 1996; Boyd, 2001). Numerical convergence is guaranteed by comparing eigenvalues for different grid resolutions and eliminating the spurious ones. Eigenvalues that are not symmetric with respect to the imaginary axis are also eliminated. In the symmetric subspace defined in section 2, only half of the disc needs to be considered and boundary conditions for perturbations are imposed at the midplane. Thus , , and and are forced to be 0 in order to be consistent with the background equilibrium symmetry. We checked that the eigenfunctions have their momentum and horizontal magnetic field components that tend to 0 as goes to infinity when such boundary conditions are used.

To obtain the marginal stabiliy border of Fig. 10, we used a different but complementary approach. These curves are obtained by computing the marginally stable equilibria directly by a shooting method. The ordinary differential equations governing the vertical structure of a symmetric equilibrium and a symmetric marginally stable linear eigenmode are integrated from an upper boundary (at several scaleheights above the midplane) to the midplane. At the upper boundary, the horizontal components of the magnetic field are set to specified values, and the density is set to a very small specified value. The value of , and the remaining unknown quantities at the upper boundary, are adjusted by Newton-Raphson iteration so that the symmetry conditions at the midplane are satisfied.

Appendix C MRI stability in a disc compressed by a strong toroidal field

An analytical model can be developed for the equilibrium and stability of a disc that is compressed by a strong toroidal magnetic field associated with a wind torque. This corresponds to the limit of the problem