Strongly magnetized accretion disks: structure and accretion from global magnetohydrodynamic simulations

TL;DR

This study uses global magnetohydrodynamic simulations to explore how net vertical magnetic fields influence the structure and accretion processes in geometrically thin disks, revealing magnetic pressure's significant role and the presence of winds.

Contribution

It provides the first comprehensive analysis of the impact of net vertical magnetic flux on thin accretion disk structure using global simulations, highlighting magnetic effects beyond local models.

Findings

Magnetic pressure significantly influences vertical disk structure.

Accretion occurs mainly in surface layers driven by large-scale magnetic stress.

Winds are present but not the main driver of disk evolution.

Abstract

We use global magnetohydrodynamic simulations to study the influence of net vertical magnetic fields on the structure of geometrically thin ($H/r \approx 0.05$) accretion disks in the Newtonian limit. We consider initial mid-plane gas to magnetic pressure ratios $\beta_0 = 1000,\, 300$ and $100$, spanning the transition between weakly and strongly magnetized accretion regimes. We find that magnetic pressure is important for the disks' vertical structure in all three cases, with accretion occurring at $z/R\approx 0.2$ in the two most strongly magnetized models. The disk midplane shows outflow rather than accretion. Accretion through the surface layers is driven mainly by stress due to coherent large scale magnetic field rather than by turbulent stress. Equivalent viscosity parameters measured from our simulations show similar dependencies on initial $\beta_0$ to those seen in shearing box simulations, though the disk midplane is not magnetic pressure dominated even for the strongest magnetic field case. Winds are present but are not the dominant driver of disk evolution. Over the (limited) duration of our simulations, we find evidence that the net flux attains a quasi-steady state at levels that can stably maintain a strongly magnetized disk. We suggest that geometrically thin accretion disks in observed systems may commonly exist in a magnetically "elevated" state, characterized by non-zero but modest vertical magnetic fluxes, with potentially important implications for disk phenomenology in X-ray binaries (XRBs) and active galactic nuclei (AGN).

Equations

$$\rho (R,z)=\rho (R,z=0)\exp \left[\frac{GM}{c_s^2}\left(\frac{1}{\sqrt{R^2+z^2}}-\frac{1}{R}\right)\right],$$

$$\rho (R,z=0)=\rho (R_0,z=0){\left(\frac{R}{R_0}\right)}^p.$$

$$T(R,z)=T(R_0){\left(\frac{R}{R_0}\right)}^q,$$

$$v_{\varphi }(R,z)=v_K{\left[(p+q){\left(\frac{c_s}{v_{\varphi ,k}}\right)}^2+1+q-\frac{qR}{\sqrt{R^2+z^2}}\right]}^{1/2},$$

$$\frac{dE}{dt}=-\frac{E-c_{\nu }\rho T_0}{t_{\mathrm{cool}}},$$

$$\rho_{\mathrm{fl}}=\left\{\begin{array}[]{ll}\rho_{\mathrm{fl,0}}\left(\frac{R}{R_{0}}\right)^{p}\left(\frac{1}{z^{2}}\right)\,\,\mathrm{if}\,\,R>r_{\mathrm{min}}\\ \rho_{\mathrm{fl,0}}\left(\frac{r_{\mathrm{min}}}{R_{0}}\right)^{p}\left(\frac{1}{z^{2}}\right)\,\,\mathrm{if}\,\,R<r_{\mathrm{min}}\,\,\mathrm{and}\,\,r>3r_{\mathrm{min}}\\ \rho_{\mathrm{fl,0}}\left(\frac{R}{R_{0}}\right)^{p}\left(\frac{1}{z^{2}}\right)\left(5-2\frac{r-r_{\mathrm{min}}}{r_{\mathrm{min}}}\right)\left(4\frac{r_{\mathrm{min}}-R}{r_{\mathrm{min}}}+1\right)\\ \mathrm{if}\,\,r<r_{\mathrm{min}}\,\,\mathrm{and}\,\,r<3r_{\mathrm{min}}\end{array}\right.$$

$$A_{\varphi }=\frac{B_0}{R_0^m}\frac{R^{m+1}}{m+2}+\frac{B_0{r}_{\min }^{m+2}}{R_0^m}\left(\frac{1}{2}-\frac{1}{m+2}\right),$$

$$⟨a(R,z){⟩}_{\varphi }=\frac{1}{2\pi }{\int }_0^{2\pi }a(R,z)d\varphi .$$

$$⟨a{⟩}_t=\frac{1}{(t_f-t_i)}{\int }_{t_i}^{t_f}adt$$

$$T_{\mathrm{r}\varphi ,\mathrm{Rey}}$$

$$T_{\mathrm{r}\varphi ,\mathrm{Max}}$$

$$T_{\theta \varphi ,\mathrm{Rey}}$$

$$T_{\theta \varphi ,\mathrm{Max}}$$

$$⟨T_{r\varphi ,\mathrm{Max}}^{\mathrm{Coh}}{⟩}_{\varphi }=-⟨B_{\mathrm{r}}{⟩}_{\varphi }⟨B_{\varphi }{⟩}_{\varphi }$$

$$⟨T_{r\varphi ,\mathrm{Max}}^{\mathrm{Turb}}{⟩}_{\varphi }=⟨-B_{\mathrm{r}}{B}_{\varphi }{⟩}_{\varphi }+⟨B_{\mathrm{r}}{⟩}_{\varphi }⟨B_{\varphi }{⟩}_{\varphi }.$$

$$\alpha =\frac{⟨T_{\mathrm{r}\varphi }{⟩}_v}{⟨P_{\mathrm{gas}}{⟩}_v}.$$

$${\alpha }_{\mathrm{SS}}=\frac{⟨T_{r\varphi ,\mathrm{Max}}^{\mathrm{Turb}}{⟩}_v}{⟨P_{\mathrm{gas}}{⟩}_v}.$$

$${\dot{M}}_{\mathrm{acc}}(r)={\int }_0^{\pi }{\int }_0^{2\pi }\rho v_r{r}^2\sin \theta d\theta d\varphi .$$

$${\dot{M}}_{\mathrm{wind}}(r)={\int }_0^{70^{\circ }}{\int }_0^{2\pi }2\rho v_r{r}^2\sin \theta d\theta d\varphi .$$

$$j_{\mathrm{r}}(r)={\int }_{70^{\circ }}^{110^{\circ }}{\int }_0^{2\pi }\left(\rho v_{\mathrm{r}}{v}_{\varphi }R-R{B}_{\mathrm{r}}{B}_{\varphi }\right)r^2\sin \theta d\theta d\varphi$$

$$j_{\theta }(r)=2{\int }_{0.1}^R{\int }_0^{2\pi }\left(\rho v_{\theta }{v}_{\varphi }R-R{B}_{\theta }{B}_{\varphi }\right)rdr\sin \theta d\varphi .$$

Full text

††pagerange: Strongly magnetized accretion disks: structure and accretion from global magnetohydrodynamic simulations–LABEL:lastpage

1 Introduction

The $\alpha$-viscosity prescription for accretion disks (Shakura & Sunyaev, 1973) — which sets the vertically integrated radial stress proportional to pressure — assumes that angular momentum transport is driven by turbulence and related to the magnetic field. The physical basis for this prescription in the magnetorotational instability (MRI: Balbus & Hawley, 1991, 1998) seems well-supported by simulations. However, geometrically thin disk models based on the $\alpha$-prescription suffer numerous failures in matching the phenomenology of disks in accreting systems such as cataclysmic variables (CVs), XRBs and AGN. For example, $\alpha$-models were found to be thermally (Shakura & Sunyaev, 1976; Piran, 1978) and viscously (Lightman & Eardley, 1974) unstable in the inner, radiation pressure-dominated regions of luminous disks, yet observations of black hole XRBs indicate stable accretion all the way to the innermost radius. In CVs and XRBs, spectral signatures indicate disks that are hotter and geometrically thicker than predicted by the standard theory (Begelman & Pringle, 2007, and references therein).

The standard thin disk model also fails to support the expected accretion flows around AGN. The outer regions of thin AGN disks are expected to be gravitationally unstable, fragmenting into clumps that trigger star formation. This curtails accretion and hence fails to explain not only the fueling of AGN but also the growth of supermassive black holes (Kolykhalov & Syunyaev, 1980; Shlosman & Begelman, 1987; Goodman, 2003). More directly, observations of AGN disk sizes, spectra and variability timescales also point to hotter and geometrically thicker disks in the inner regions (Dexter & Begelman, 2019).

These difficulties can be ameliorated and perhaps resolved if accretion disks are geometrically thicker and faster-accreting than predicted by the standard Shakura & Sunyaev (1973) theory. The development of strong toroidal ($B_{\phi}$) magnetic fields, amplified by the dynamo process associated with MRI, provides a plausible way to obtain such thickening by exerting magnetic pressure in the vertical direction. A strongly magnetized disk, with a magnetic pressure larger than the gas or radiation pressure, would have a lower density than a standard thin disk with the same accretion rate, thus avoiding gravitational instability (Pariev et al., 2003; Begelman & Pringle, 2007; Gaburov et al., 2012). Lower densities could increase the color correction (Blaes et al., 2006), leading to higher temperatures, while geometric thickening could increase the degree of radiative reprocessing (also contributing to higher temperatures) and accretion inflow speeds. Finally, the decoupling of radiative transport, dissipation and vertical pressure support would eliminate the thermal-viscous instability (Begelman & Pringle, 2007; Sa̧dowski, 2016).

MRI-driven turbulence tends to amplify the $\phi$-component of the magnetic field more than the other components, but when the net flux of the disk is negligible, the saturated pressure of $B_{\phi}$ remains small compared to the background pressure. However, this changes if the MRI is seeded by a net poloidal (vertical) magnetic flux of sufficient magnitude. Hawley et al. (1995) performed three-dimensional local (shearing box) MHD simulations of the growth and saturation of MRI seeded by poloidal magnetic field of different strengths. They found that angular momentum transport was dominated by magnetic (Maxwell) stresses rather than Reynolds stress, and characterized by an effective $\alpha$ that increases with increasing initial magnetic field strength. Recent stratified shearing box simulations with higher poloidal fluxes extend the trend measured by Hawley et al. (1995) to ratios of gas to magnetic pressure as low as $\beta_{0}=100$ (Bai & Stone, 2013) and $\beta_{0}=10$ (Salvesen et al., 2016). These studies find that values of $\beta_{0}\lesssim 1000$ lead to a strongly magnetized disk midplane and regions of large Maxwell stress extending to several scale heights on either side.

Quasi-periodic alternations in the disk’s large scale magnetic field (often referred to as a “dynamo") were observed in vertically stratified local shearing box simulations (Brandenburg et al., 1995), and are characteristic of the MRI. In such a dynamo process, large scale toroidal magnetic fields are generated which show a “butterfly diagram” with even parity across the disk midplane and flips in orientation every $\sim 30$ local orbital periods. Although the dynamo develops even with negligible poloidal flux, the characteristics of the field reversals are governed by the initial magnetic field strength. A stronger net poloidal magnetic field delays the reversal of $B_{\phi}$, leading to field structures that resemble long-lived global patterns rather than turbulence (Bai & Stone, 2013; Salvesen et al., 2016).

Early global simulations of unstratified (cylindrical) disks that included a locally non-zero vertical field yielded moderately large $\alpha$ values $\sim 0.1$ (Armitage, 1998), qualitatively in accord with local results. Most stratified simulations, however, have focused on the zero net flux regime (Mishra et al., 2016; Sa̧dowski, 2016; Hogg & Reynolds, 2018). Simulations with a significant net field require a smaller time-step to follow regions with low plasma $\beta$, and also require careful treatment of the inner boundary to avoid numerical artefacts. With a careful choice of the inner boundary condition, Zhu & Stone (2018) were able to evolve a fully global numerical MHD model for two different thermal scale heights ($H/r=0.1$ and $H/r=0.05$). Zhu & Stone (2018) used seed magnetic fields with $\beta_{0}=10^{4}$ and $1000$, with the stronger field case barely (based on local expectations) reaching the boundary of the strongly magnetized regime.

In this paper, we use global simulations to study disks whose net vertical flux is chosen (based on local results) to span the transition from weak to strong magnetic field strengths. Our goal is to test whether the local simulation results carry over to the global regime, and hence whether strongly magnetized disks are potentially relevant to the open problems discussed above. We ran a weak magnetic field case simulation ($\beta_{0}=10^{3}$) for 23 orbits (at a fiducial radius $R=1$). Intermediate ($\beta_{0}=300$) and strong ($\beta_{0}=100$) magnetic field cases are run for 23 and 50 orbits, respectively. We note that Zhu & Stone (2018) attempted a model similar to our strong field simulation, but observed extremely rapid mass loss in that case. We find that this problem can be alleviated with a change to the assumed radial temperature structure of the disk, allowing physical results to be obtained in the strong field case.

The paper is organized as follows. In Section  2, we describe the numerical setup of the disk which includes the hydrodynamical disk profile and magnetic field configurations. In Section  3 we present detailed findings from our numerical simulations. The results are categorized into vertical structure (Section 3.1), radial structure (Section 3.2), wind properties (Section 3.3), magnetic flux evolution (Section 3.4) and comparison with local (Section 3.5) shearing box simulations and other global (Section 3.6) simulations. In Section 4, we conclude and discuss the astrophysical applications and relevance of our numerical model.

2 Numerical Setup

We model the accretion disk using the Athena++ (White et al., 2016) code in the ideal MHD limit. Athena++ uses a higher-order Godunov scheme and employs constrained transport (CT) to maintain divergence-free evolution of the magnetic field. We use the HLLD solver implemented in Athena++. Our initial conditions and methods are generally similar to those used by Zhu & Stone (2018), with the main difference being that we model a disk with a radially constant disk aspect ratio ($H/r$, where $H$ is the thermal scale height of the disk). The density profile of the disk in hydrostatic equilibrium in cylindrical coordinates is given by

[TABLE]

where the initial midplane density profile is

[TABLE]

The initial temperature is taken to be constant on cylinders,

[TABLE]

and is maintained constant throughout the simulation using the same methods as in Zhu & Stone (2018). The initial velocity profile in the disk has only an azimuthal component given by

[TABLE]

where $R$ and $z$ are cylindrical radius and height. We choose $p=-1.5$ and $q=-1$, in contrast to the values $p=-2.25$ and $q=-0.5$ chosen by Zhu & Stone (2018). We are interested in geometrically thin disks and choose an initial disk aspect ratio of $H/r=0.05$. Note that we wrote the disk profile in cylindrical coordinates for simplicity but the simulations are performed in spherical polar coordinates.

We maintain effectively isothermal conditions by rapidly relaxing the temperature to a pre-specified value at each cylindrical radius. The cooling rate is computed by

[TABLE]

where $E$ is internal energy per unit volume and $c_{\nu}$ is heat capacity per unit mass (Zhu et al., 2015). The rapid cooling time is assumed to be the numerical time step of the simulation. We remove any excess internal energy in the simulation and maintain the same temperature as its initial value at each cylindrical radius.

Strong magnetic field leads to small time steps in the inner polar regions of the disk. In order to maintain a reasonable time step we used the radial power-law density profile of floor values prescribed in Eq. 10 of Zhu & Stone (2018):

[TABLE]

where, $\rho_{\mathrm{fl,0}}=10^{-9}$ is a constant. The simulation domain extends from $r_{\mathrm{min}}=0.1$ to $r_{\mathrm{max}}=100$ in code units. The purpose of keeping the outer boundary far away is to avoid any effects of the outer boundary conditions on the simulation region of interest. The polar and azimuthal domains extend from $\theta_{\mathrm{min}}=0$ to $\theta_{\mathrm{max}}=\pi$ and $\phi_{\mathrm{min}}=0$ to $\phi_{\mathrm{max}}=2\pi$. The radial grid is evenly spaced in logarithmic space. We use the static mesh refinement (SMR) feature of Athena++ to improve the resolution of the dense regions of the disk. We add three levels of static mesh refinement over a base resolution of $256\times 32\times 128$, giving an effective resolution of $2048\times 256\times 1024$.

The static mesh refinement levels (shown in Fig. 1) are applied as follows: all refinement levels extend to the minimum radius $r=0.1$, which is the inner boundary of the simulation domain. The outer radial limits for the first, second and third refinement level are at $r=40,\,30,\,20$, respectively. The refinement domains in $\theta$ are applied symmetrically (Fig. 1) about the equator but Athena++ requires slightly different minimum and maximum values of $\theta$ for best mesh-block configuration of the simulation domain. This leads to the following numbers for $\theta$. The first refinement level extends from $\theta=0.8$ to $1.8$, the second refinement level extends from $\theta=1.2$ to $1.77$ and the third refinement level covers the narrow domain of $\theta=1.4$ to $1.6$ (all in radians). We use an outflowing boundary condition radially for gas and for magnetic field. The azimuthal boundary condition is periodic. In order to avoid significant flux loss from the polar region we use the polar boundary condition described in Zhu & Stone (2018). The polar boundary condition in a full $2\pi$ azimuthal domain acts in a way that assigns cell- and face-centered values into the ghost cell in $\theta$ using the opposite active cell in $\phi$.

The use of an effectively radially isothermal equation of state (and the neglect of physics such as radiative transport) means that the simulations are scale-free. We report all results in code units. We take $G=M=1$ and report densities and velocities in terms of the midplane density and Keplerian velocity at $R=1$.

2.1 Magnetic Field

As illustrated in Fig. 2, we can identify three regimes of disk magnetization, all of which are potentially of astrophysical interest. Weakly magnetized disks are those whose poloidal field, even if non-zero, is weak enough that the properties of the MRI are essentially unaltered from the strictly zero net flux limit. Strongly magnetized disks, the subject of this study, are those where the poloidal field, via its ability to source larger toroidal fields, is large enough to affect disk structure qualitatively in the vertical direction, but not so large as to quench MRI and dominate the radial force balance. Finally, magnetically arrested disks (MAD; Narayan et al., 2003; McKinney et al., 2012) occur when the poloidal field on its own is strong enough to dominate the disk dynamics.

The results of local simulations (Salvesen et al., 2016) imply that the transition between weak and strong magnetization occurs for initial values of (poloidal) $\beta=P_{\rm gas}/P_{\rm b}$ in the range between $10^{3}$ and $10^{2}$. We therefore seed the initial disk with a purely vertical magnetic field with a uniform ratio of gas pressure to magnetic pressure as a function of radius at the disk midplane. To ensure an initially divergenceless magnetic field, we initialize the magnetic field via the magnetic vector potential defined in Zhu & Stone (2018),

[TABLE]

where $m=(p+q)/2$ with $p=-1.5$ and $q=-1$. We adopt three values of initial midplane $\beta=P_{\mathrm{gas}}/P_{\mathrm{b}}$, corresponding to weak, intermediate and strong initial magnetic field strength cases with initial $\beta_{0}=1000\,,300$ and $100$, respectively. We add a sinusoidal component to the initial vertical magnetic field to suppress the growth of channel flows in the strongly magnetized cases. A model similar to our $\beta_{0}=1000$ case (although using a flaring disk profile) has been studied by Zhu & Stone (2018) and serves as a point of comparison.

3 Results

We first define the post-processing steps that we employ. An azimuthally averaged quantity $a$ at fixed radius is defined by

[TABLE]

Similarly, a time averaged quantity is defined by

[TABLE]

where the time averaging is performed between $t_{i}$ and $t_{f}$. We ran the $\beta_{0}=100$ case for $50$ orbits at $R=1$, while the total runtime of the $\beta_{0}=1000$ and $\beta_{0}=300$ simulations is 23 orbits at the same radius. We typically present time averages between $t_{i}=19$ and $t_{f}=23$. All the vertical profiles reported are azimuthally averaged at $R=1$ (similar to Zhu & Stone, 2018). The radial and polar Reynolds and Maxwell stresses are defined by,

[TABLE]

and,

[TABLE]

respectively. The Maxwell stress defined above is normalized with a factor of $4\pi$ to be consistent with the Athena++ definition of the magnetic field strength. Here $B_{r}$, $B_{\theta}$ and $B_{\phi}$ correspond to components of the magnetic field in spherical polar coordinates, and $\delta v_{\phi}$ corresponds to velocity fluctuations defined by $\delta v_{\phi}=(v_{\phi}-\langle v_{\phi}\rangle_{\phi})$, where $\langle v_{\phi}\rangle_{\phi}$ is azimuthally averaged toroidal velocity.

The total radial Maxwell stress $T_{r\phi\,,\mathrm{Max}}$ can be decomposed into coherent and turbulent components. The coherent component is defined by,

[TABLE]

and the turbulent component is,

[TABLE]

The viscosity parameter $\alpha$ is defined as the ratio of the volume averaged total radial stress (including turbulent and coherent components) to the similarly averaged gas pressure

[TABLE]

We distinguish this from the “Shakura-Sunyaev" $\alpha_{\rm SS}$,

[TABLE]

which is a similar ratio based only on the turbulent part of the stress, i.e. (Eq. 17). We will see later that the coherent stress dominates the accretion in magnetically dominated disks, leading to $\alpha\gg\alpha_{\rm SS}$ (see Sect 3.1.2).

3.1 Vertical Structure

The basic structure of the simulated disks is shown in Fig. 3 for the most strongly magnetized disk with $\beta_{0}=100$. This time slice was taken at $t=50$ orbits at $R=1$, at which time the disk and its magnetic field structure had achieved a quasi-steady state (inflow equilibrium) inside $R\approx 3$. Three regions of the flow can be distinguished. Near the disk midplane the magnetic pressure is relatively weak, because of the toroidal field reversal that must occur there when the flow is threaded by net vertical flux. The gas near $z=0$ is outflowing, accompanied by a pushing or pinching of the poloidal magnetic field outwards. Above the midplane, strong accretion occurs in magnetic pressure supported regions, which lie at $z/r\approx 0.2-0.4$ in our models. The poloidal magnetic field is pinched radially inward in this zone of the flow. Finally, the polar or atmosphere regions above the main body of the disk exhibit outflow. Due to the low density in this region (approximately six orders of magnitude lower compared to the density in accreting regions) accretion greatly dominates outflow, overall.

The primary finding of surface or elevated accretion is consistent with previous studies based on global numerical simulations (Suzuki & Inutsuka, 2014; Zhu & Stone, 2018; Jiang et al., 2019). However, it is physically different from the structures that had been suggested on the basis of local simulations (Bai & Stone, 2013; Salvesen et al., 2016) and analytic models (Begelman et al., 2015). Although similar values of $\beta_{0}$ trigger the onset of magnetic pressure support for the vertical structure, the toroidal field reversal at $z\approx 0$ seen globally means that our magnetically elevated disks can remain gas pressure dominated in the disk midplane. This physical configuration for the toroidal field can occur in local simulations, but is not guaranteed to be present.

3.1.1 Profiles of the disk variables

Fig. 4 shows the vertical structure of the simulated disks. The top panel shows the vertical profile of the density, which has a maximum at the disk midplane and a steeply decreasing profile up to $z/R~{}\approx 0.1$. The density drops more gradually (compared to a Gaussian profile) at higher altitudes of the disk. The density is highest in the elevated region for the strong field case and lowest in the weak field case. The results clearly indicate that the extent to which a disk will be magnetically elevated depends on the initial magnetic field strength.

The middle panel of Fig. 4 shows the vertical profile of the radial mass flux density computed by azimuthally and time averaging (over orbits 19–23 for each magnetic field strength case) the product of the density and the radial component of velocity. The inflowing region appears at higher altitudes ($z/r\approx 0.2$, see Fig. 3), while the midplane shows outflow. The geometrically thin disk ($H/r=0.05$) case of Zhu & Stone (2018) showed similar outflowing regions close to the midplane. The strong field case shows the largest accretion mass flux and the weak field case the smallest accretion flux close to $z/r\approx 0.2$. Unlike the strong field case and intermediate field case, the weak field model shows dominant accretion at $z/R\approx 0.1$ with a secondary smaller inflow region at $z/R\approx 0.3$.

The lowest panel of Fig. 4 shows the vertical profile of azimuthally and time averaged total $\beta_{t}$, obtained by taking the ratio of the azimuthal averages of gas and magnetic pressure at $R=1$. The disk is initialized with an initial $\beta_{0}=1000,\,300$ and $100$, and saturates with a final disk midplane total $\beta_{t}=100,\,10$ and $2$ for the weak, intermediate and strong case, respectively. At higher altitudes in the disk, where accretion is preferentially occurring, magnetic pressure dominates in all three cases. In the accreting regions ($H/r\approx 0.2$), $\beta_{t}$ for the weak field case is approximately 4. The intermediate case achieves a final $\beta_{t}\approx 0.2$ at this height, and the strong field case attains a final $\beta_{t}\approx 0.08$. We find that the vertical profiles of total $\beta_{t}$ in our global simulations differ from those in the local shearing box simulations of Salvesen et al. (2016). The local shearing box simulations can achieve a much more strongly magnetized disk midplane because there is not necessarily an equatorial current sheet, in contrast to the global simulations (see Section 3.5).

The azimuthally and time averaged vertical profile of $B_{\phi}$ (upper panel Fig. 5) shows the dependence of the time-evolved toroidal magnetic field on the initial poloidal magnetic field strength. The weak magnetic field case ends up with a weak toroidal field and the strong magnetic field case leads to a stronger resultant toroidal magnetic field strength. Close to the disk midplane, $B_{\phi}$ is zero due to the formation of a current sheet there. This also gives intuitive and qualitative insight into how the toroidal field is being generated in these numerical models. Once the MRI triggers, an initial poloidal field develops a radial and toroidal magnetic field. MRI amplifies the toroidal field (Latter et al., 2010; Lesur et al., 2013) and leads to a magnetic pressure dominated accreting region. In our chosen spherical polar coordinate system, a sheared radial and poloidal component of the magnetic field will produce a positive and a negative $B_{\phi}$ above and below the disk midplane, respectively. This is consistent with previous global simulations reported in Zhu & Stone (2018), but differs in significant respects from older models in which large-scale magnetic torques drove predominantly midplane accretion (Pudritz, 1985). The lower panel of Fig. 5 shows $B^{2}_{\phi}$ to illustrate the gradual decrease of toroidal magnetic pressure with altitude. A dip in the vertical profile of $B^{2}_{\phi}$ close to the disk midplane is due to the current sheet.

3.1.2 Stress profiles in the disk

The time averaged viscosity parameter $\alpha$ is shown in the upper panel of Fig. 6 for each of our models. We emphasize that this value of $\alpha$ includes the total stress (sum of radial Reynolds and Maxwell stresses), in contrast to the standard $\alpha_{\rm SS}$ (Eqn. 17) which only includes the turbulent component. We find the expected trend that the weak initial magnetic field case gives a lower value of the viscosity parameter $\alpha$ and the strong magnetic field case gives larger values of $\alpha$. In the disk midplane, the difference in measured $\alpha$ from each of the simulations is small but in the elevated regions, where most of the accretion is occurring, the difference between inferred $\alpha$ between strong and weak field cases is much larger. The weak magnetic field case gives $\alpha\approx 1$ at $z\approx 0.2$, the intermediate initial magnetic field strength case gives $\alpha\approx 5$ and the strong magnetic field case gives $\alpha\approx 10$. Note that these $\alpha$ parameters are normalized to the gas pressure, which is much lower than the magnetic pressure in the elevated layers for the intermediate and strongly magnetized cases. If the viscosity parameter is instead defined normalized to the magnetic pressure, all three cases have similar values in the accreting layers.

In the middle panel of the Fig. 6 we show the time and azimuthally averaged total radial stress ($T_{r\phi}=T^{\mathrm{Rey}}_{r\phi}+T^{\mathrm{Max}}_{r\phi}$), which is dominated by Maxwell stress. The strong magnetic field case has almost a factor of eight larger value of total radial stress at $z/R\approx 0.2$ than the weak field case. The intermediate field strength case lies in between the two extreme cases of initial field strength. The polar stress ($T_{\theta\phi}$) is almost a factor of four smaller than the radial stress in the accreting regions. The vertical profile of the radial stress is always positive, implying a net outward angular momentum flux for Maxwell stress dominated flows. Above the disk midplane the azimuthal component of the magnetic field is positive and the radial component is negative (streamlines in Fig. 3) leading to positive radial stress above the disk midplane. Similarly, the vertical net flux seeding the magnetic field in our model leads to a positive radial component of the magnetic field (Fig. 3) and a negative azimuthal component below the midplane, also causing positive radial stress. This profile of radial stress ensures that there should be net accretion within the disk regions.

The lowest panel in Fig. 6 shows the vertical profile of $T_{\theta\phi}$. A negative/positive value of $T_{\theta\phi}$ above/below the disk midplane implies an outward angular momentum carried away by winds. $T_{\theta\phi}$ behaves differently in the regions close to the disk midplane compared to higher altitudes. The green curve in this figure shows large amplitude of negative and positive $T_{\theta\phi}$ above and below the disk region, respectively. The amplitude of $T_{\theta\phi}$ for higher altitudes in intermediate and weak field cases is much smaller compared to the strong field case. This indicates that the strong field case can launch winds more efficiently compared to the weak field case. Close to the disk midplane, $T_{\theta\phi}$ is pointing towards the disk midplane (specifically in the intermediate field case shown with solid black curve) and this could lead to angular momentum deposition on the disk midplane, causing outflows. These midplane outflows are not surprising because such backflows have also been predicted earlier in accretion disk solutions reported in Kluzniak & Kita (2000).

The total radial and polar stress is dominated by Maxwell stress. The total Maxwell stress can be decomposed into two components, one coherent (Eq. 14) and the other turbulent (Eq. 15). The streamlines of the magnetic fields in Fig. 3 show a coherent field configuration in elevated accreting regions. In order to quantify it further, we plot the total Maxwell stress, its turbulent component and coherent component in Fig. 7. Comparing the middle panel of Fig. 6 with the topmost panel of Fig. 7 it is evident that the total radial stress is dominated by the Maxwell stress. Further decomposition of the Maxwell stress shows that it is dominated by the coherent component (shown in the middle panel of Fig. 7) everywhere except close to the midplane. Close to the disk midplane the coherent component vanishes. On the other hand, the turbulent component shown in the lowest panel dominates in the disk midplane but is about a factor of four smaller compared to the coherent component in the accreting regions. The turbulent $\alpha_{\mathrm{SS}}$ can be easily calculated given the total $\alpha$ in Fig. 6 and the turbulent Maxwell stress component. We find that $\alpha$ due to the coherent component dominates over turbulent $\alpha_{\mathrm{SS}}$. Our strong field and intermediate field cases are thus in qualitative agreement with Zhu & Stone (2018). The weak field case is quantitatively in agreement with Zhu & Stone (2018), where the coherent component of the stress is a factor of two larger than the turbulent stress in accreting regions.

3.2 Radial Structure

In Fig. 8, we show off-equatorial ($z/R=0.2$) slices of the mass density plotted at $t=20$ orbital periods (measured at $R=1$). The upper panel (initial $\beta_{0}=1000$) shows a weak spiral azimuthal density profile. The middle panel and lower panels show the density profile for initial $\beta_{0}=300$ and $\beta_{0}=100$, respectively. As we increase the initial magnetic field strength, the disk shows dramatically stronger spiral structures with an increasingly inhomogeneous density profile.

In Fig. 9 we plot the toroidal component of the magnetic field at higher altitude ($z/R=0.2$) where most of the accretion is occurring. We only show the toroidal field slice for the strong and intermediate field cases. The weak field case shows a corresponding density profile but does not show significant inhomogenities. The color-bar shows the range of the toroidal magnetic field using (blue, red) for (negative, positive) $B_{\phi}$.

The strong ($\beta_{0}=100$) and intermediate ($\beta_{0}=300$) magnetic field cases show spiral structures that can be approximately fitted with logarithmic spirals of the functional form $r=ae^{-b\phi}$, where $a$ and $b$ are constant coefficients for a chosen setup. In Fig. 9, the black curves show the fits to the spiral arms in the azimuthal component of the magnetic field, with $b=0.3$ in the upper panel ($\beta_{0}=100$) and $b=0.088$ in the lower panel ($\beta_{0}=300$). This trend is consistent with $b\propto\beta_{0}^{-1}$. The fit is created by taking the logarithm of $B_{\phi}$ to eliminate the negative component and distinguish or isolate positive $B_{\phi}$ regions in Fig. 9. We trace the the ridge of the field strength and fit a logarithmic spiral through it. A detailed study of the formation of these spirals will be reported in an upcoming paper.

The midplane radial profiles of the azimuthally averaged density, mass flux, poloidal field ($B_{\theta}$) and toroidal field ($B_{\phi}$) are shown in Figs. 10 and 11, respectively. The weak and intermediate field cases have higher density at the disk midplane compared to the strong field case. Since the strong field case accretes at higher rate, the low midplane density simply reflects the fact that matter is draining more quickly.

The lower panel of Fig. 10 shows the product of the density and radial component of velocity at the disk midplane. We notice both outflowing ($R>0.95$) as well as inflowing ($R<0.95$) structures. This is consistent with the vertical profile of the same variable (Fig. 4). Accretion is mainly occurring in the surface of the disk (see Fig. 3), with negligible accretion along the disk midplane (except very close to the inner boundary). The vertically integrated mass flux in Fig. 12 also shows a dominant accretion across all the radii and a significantly large inflow inside ($R<0.5$), which is very close to outflowing inner boundaries.

In Fig. 11, we find an enhanced magnetic field strength in the inner regions. Note that we overplot the initial profile of the magnetic field strength $B_{\theta}\propto R^{-5/4}$ (and $B_{\phi}=0$) as thin curves. The poloidal component of the magnetic field has flipped its orientation at $R<1.0$ for each magnetic field case. The strong field case also reverses its sign at larger radii. The bottom panel shows that the toroidal field strength dominates over the poloidal component. In both poloidal and toroidal magnetic field radial profiles, we find sign reversals. This has also been seen in previous global simulations of weakly magnetized disks by Zhu & Stone (2018). As we discussed earlier in connection with the vertical stress profiles, the disk midplane is dominated by turbulent stress, leading to large fluctuations in magnetic field orientation.

3.3 Accretion, winds and outflows

We find a dominant accretion flow in all three models we simulated. However, in addition to inflow, the disk also exhibits winds (in polar regions) and outflows (disk midplane). If we closely look into the flow above or below the disk midplane, it has three distinct regions as shown in the top panel of Fig. 14. The disk midplane shows outflows followed by elevated accreting regions at $z/R=0.2$. Above these accreting regions, there are winds which are distinct from the midplane outflows. This suggests that the vertical structure of the disk is much more complicated than just turbulent inflowing gas in the disk and winds along the polar regions. In this subsection, we focus on each of these individual flows. First we discuss accretion, secondly winds and finally midplane outflows.

In order to compute the mass accretion rate, we used the vertically integrated mass flux at each spherical radius,

[TABLE]

We evolved the strong field case for $50$ orbital periods at $R=1$, which gives a steady state accretion flow up to $R\approx 3$. This can be seen in Fig. 12, where we show the radial profile of mass accretion rate for the strong field case. The three different curves in Fig. 12 show an almost horizontal radial profile of the mass flux for three time windows separated by 10 orbital periods (measured at $R=1$) each. We also note that, despite large outflow velocity magnitudes at high latitudes shown in Fig. 3, the flow is dominated by accretion. A large inflow rate in the $R<0.5$ regions is influenced by inner outflowing boundary conditions. This region is also influenced by density floor values applied to maintain a reasonable time step ($10^{-7}$) and to avoid magnetic flux loss near the poles.

It is also important to estimate the mass and angular momentum carried away by winds from the accretion disk. In order to compute the mass outflow rate from the wind regions, we use a polar domain spanning from $\theta=0$ to $\theta=70^{\circ}$.

[TABLE]

The factor of 2 in Eq. 19 accounts for the wind mass loss rate from the lower hemisphere. The lower panel of Fig. 12 shows that in our strong field case, mass loss due to the wind varies between $2-6\%$ of the total accreted mass at different time windows separated by 10 orbital periods at $R=1$. This also suggests that winds are not steady but time varying. A time averaged mass loss rate due to winds (Fig. 13) shows that about $6\%$ of the accreted mass is carried away by winds.

In computing the mass loss rate due to winds, we encounter low resolution regions close to the boundary of the third and second refinement levels, leading to small fluctuations in the radial profile of the wind mass loss rate. Despite this caveat due to low resolution in the disk wind regions, we find an order of magnitude larger wind driven mass loss rate compared to Zhu & Stone (2018) for their initial $\beta=1000$ case. This suggests that strongly magnetized accretion disks can launch winds more efficiently compared to weakly magnetized accretion disks.

We further quantify this in Fig. 13 by plotting the mass accretion rate and mass loss rate due to winds as a function of initial $\beta_{0}$. The red triangle data points show mass accretion rate ($\dot{M}_{\mathrm{acc}}$) and blue data points show mass loss rate ($\dot{M}_{\mathrm{wind}}$) computed at $r=1$. Note the different scales for mass accretion and loss. Both mass accretion and mass loss rates decrease as we decrease the magnetic field strength.

We also estimate the angular momentum budget in winds and the disk region. We use the following diagnostics, to compute the radial angular momentum flux

[TABLE]

and the polar angular momentum flux

[TABLE]

Positive $j_{\mathrm{r}}$ corresponds to angular momentum flux pointing outwards whereas negative $j_{\theta}$ computed at fixed $\theta=70^{\circ}$ corresponds to flux pointing away from the disk. Comparing radial and polar angular momentum fluxes, we find that approximately $7\%$ of the angular momentum is carried away by the winds. Unlike the wind region, the disk region is dominated by accretion with outflows at the disk midplane. We also compute the angular momentum flux from a domain very close to the disk midplane by using 8 cells in $\theta$. Similar to the schematic diagram, we find that the angular momentum flux is pointing inwards at the disk midplane with fluctuations across zero. Comparing the middle panel of Fig. 14 (which shows the total outward pointing angular momentum flux between $\theta=70^{\circ}$ and $\theta=110^{\circ}$) and Fig. 15, we find positive radial flux in the former and negative in the latter. The angular momentum flux at the midplane is pointing inwards, carrying approximately 15% of the total outward pointing angular momentum flux shown in Fig. 14. This is also qualitatively consistent with Zhu & Stone (2018), where the authors found that approximately $10\%$ of the angular momentum is carried away by outflows at the disk midplane. At small radii, the midplane angular momentum flux points outwards, suggesting that the disk midplane is weakly accreting there. The black curve in Fig. 15 fluctuates across zero values at the disk midplane, causing more complex inflowing and outflowing structures at the midplane.

3.4 Quasi-steady state of Magnetic Flux

Given the large effect of net magnetic flux on disk structure, it is important to determine how the flux evolves with time as a function of disk radius (Begelman & Armitage, 2014). If the accretion flow advects magnetic field effectively, a large magnetic flux will accumulate close to the accreting object and eventually may lead to a magnetically arrested disk (Narayan et al., 2003; McKinney et al., 2012). Alternatively, rapid diffusion of the field will lead to weakly magnetized inner regions that are more prone to thermal instability (in X-ray binaries) and disk fragmentation in the AGN context. Here we address the competition between advection and diffusion that we observe in our global simulations. Our simulations include the poles and hence we do not lose magnetic flux from polar boundaries. This makes these simulations ideal to study magnetic flux evolution in strongly magnetized disks.

In Fig. 16 we show how the poloidal magnetic flux evolves in an enclosed cylinder with $0.1<R<1.0$. The green dashed curve shows the strong field case which we have evolved for $50$ orbits. The magnetic flux in this case varies about $8\%$ from $t=20$ to $t=50$ orbits. It increases from $t=15-20$ orbits and then fluctuates from $t=30-50$ orbits. The intermediate and weak field cases are evolved for only $23$ orbits but show magnetic flux variation of about $5\%$ over $18$ orbits. This behavior has also been reported in Zhu & Stone (2018), where the authors found a steady or quasi-steady state balance of advection and diffusion of the magnetic field. Fig. 3 shows that the field is advecting inward in accreting regions and being carried outward in the disk midplane.

Overall, our results indicate that the flux is neither accumulating at the center and approaching a MAD state, nor escaping from the inner disk through diffusion. Within the limited context of our specific setup, this suggests that magnetically elevated accretion may be long-lived once it is established.

3.5 Comparison with Local Shearing-box Simulations

As this is one of the first sets of global simulations to probe the strongly magnetized regime, it behooves us to compare our results with previous studies of strongly magnetized accretion, carried out within the local, shearing box approximation. In Fig. 17, we compute the correlation of the viscosity parameter $\alpha$ with the initial $\beta_{0}$. If we fit a line to the three data points from our simulations in the $\log\beta_{0}-\log\alpha$ plane, we find a slope of $-0.65$, which is similar to the value (-0.53) reported for local shearing box simulations (Salvesen et al., 2016). If we compute an analogous viscosity parameter taking magnetic pressure in the denominator of Eq. 16, the slope in the $\log\beta_{0}-\log\alpha$ plot is approximately zero, indicating that $\alpha$ in strongly magnetized disks remains unchanged if computed using magnetic pressure only.

As discussed earlier, we find that $\alpha$ values themselves (at a given $\beta_{0}$) are in approximate agreement between local and global simulations. Furthermore, the power law slope of the $\alpha$–$\beta_{0}$ fit is in broad agreement with the results of Salvesen et al. (2016), though it is slightly steeper.

Figure 17 also displays the stress normalized by the magnetic pressure, i.e., $\alpha_{\rm mag}$ (see Salvesen et al., 2016). As with the standard $\alpha$, we find excellent agreement in the value of $\alpha_{\rm mag}$, which is $\sim 0.2$–0.4 in both Salvesen et al. (2016) and this work. We also find essentially no (or very weak) dependence of $\alpha_{\rm mag}$ on $\beta_{0}$, consistent with local studies (Salvesen et al., 2016).

Upon examining the vertical profiles of relevant quantities (Fig. 4), we find general qualitative agreement between local and global simulations. In particular, gas density, pressure, and magnetic energy fall off quite rapidly away from the mid-plane for all $\beta_{0}$ values. As $\beta_{0}$ decreases, the disk column at $R=1$ tends towards being magnetically dominated, $\langle\beta_{t}\rangle_{\phi}<1$, with only a very small region of subthermal magnetic field at the midplane in our global simulations (compared with no such region in local simulations).

One exception to this agreement is in the vertical profile of the Maxwell stress, as shown in Fig. 4, which depicts both the coherent large scale stress and the turbulent stress in all of our simulations. As shown, the laminar stresses dominate nearly everywhere, apart from the midplane where the turbulent stresses dominate. This result can be contrasted with similar profiles in Bai & Stone (2013) (see their Fig. 6), which shows that for $\beta_{0}=100$ and $\beta_{0}=1000$, the turbulent stress never dominates.

Presumably, this difference in the stress profile arises from the presence of a current sheet at the midplane in our simulations. Due to the symmetry of the shearing box, as discussed in detail in, e.g., Bai & Stone (2013), such current sheets need not be present in local domains. However, in global setups, these current sheets act to dissipate large scale magnetic fields, which can lead to an increase in turbulent stresses at the mid-plane over laminar stresses.

In addition to the vertical stress profile, our global simulations differ from local simulations when considering radial and azimuthal structure, particularly in the magnetic field. While both types of simulations exhibit very strong inhomogeneities within the domain (Salvesen et al., 2016), local simulations are by definition unable to generate the coherent large-scale spiral structures that we see here.

In summary, there are both important similarities between local and global simulations, including the empirical relation between $\alpha$ and $\beta_{0}$, together with a number of results, largely related to the behavior of large-scale magnetic fields, that differ between the two approaches. These differences are not surprising given the obvious limitations and inherent symmetries of the the shearing box. They indicate that when probing the strongly magnetized limit of accretion, local shearing boxes should be treated with caution.

3.6 Comparison with previous global simulations

The closest comparison we can make to prior global simulations is with the work of Zhu & Stone (2018), who considered disks with $\beta=10^{4}$ and $1000$ for two different disk scale heights of $H/r=0.1$ and $0.05$. We also performed a $\beta=1000$ case to compare with results reported in Zhu & Stone (2018). The simulations reported here use the same code, essentially identical numerical methods, and a similar initial setup. The main distinction (apart from our consideration of lower initial $\beta$ values, leading to more strongly magnetized accretion flows) is that we adopt a disk profile with a radially constant $H/r$. We also model the full $2\pi$ in azimuth, whereas Zhu & Stone (2018) used a restricted azimuthal domain for their thin disk case. In the parameter regime where the two studies overlap, we find generally good agreement as to the predicted structure of the accretion flow. In particular, Zhu & Stone (2018) and Zhu et al. (2019) also found that the elevated accreting regions have stress dominated by coherent stress, whereas turbulent stress dominates near the midplane. Winds are sub-dominant as a driver of accretion in both sets of simulations.

Reducing the initial $\beta$ below the value of $10^{3}$ considered by Zhu & Stone (2018) results in stronger disk magnetization, and some qualitatively different results. The radial structure of strongly magnetized flows shows spiral structures which were not prominent in the weakly magnetized case. The vertical structure also shows much higher density at higher altitudes and enhanced elevated accretion flow on the disk’s surface layers. We also find a much higher viscosity parameter $\alpha$ due to our strong seed field. On the other hand, the result that accretion in the elevated regions is driven primarily by a coherent, rather than turbulent, Maxwell stress, persists across the weak-to-strong magnetization transition. More strongly magnetized disks show a larger mass loss rate due to winds, although we find that the mass and angular momentum loss due to the wind remain small even in the most highly magnetized case.

Zhu & Stone (2018) briefly mention that a numerical simulation with $H/r=0.1$ and $\beta_{0}=100$ was attempted but due to rapid mass loss the simulation gave very different outcomes to their weak field cases. We carried out our strong field run with $H/r=0.05$ for $50$ orbits without substantial mass loss (only $0.3\%$ of total mass is lost during the entire run for the strong field case). It is not clear whether the important difference here is due to the different values of the aspect ratio, or to detailed differences in the numerical setup between the two sets of simulations.

4 Summary and Discussion

We have presented global numerical simulations of idealized accretion disk flows threaded by weak but dynamically important net poloidal magnetic fields. Our runs span a range in the initial ratio of midplane gas to poloidal magnetic pressure between $10^{3}$ and $10^{2}$, which corresponds to the locally predicted transition between weakly and strongly magnetized disks. Using Athena++ with three levels of static mesh refinement, we were able to attain a resolution of approximately 26 cells per disk scale height in a fully global model with $H/r=0.05$. We find that:

• •

The simulated disk becomes magnetic pressure supported at higher altitudes with accretion mainly occurring in the elevated regions at $z/R\approx 0.2$. Accretion in these surface layers is mainly driven by the coherent component of Maxwell stress rather than the turbulent component. The disk midplane, on the other hand, exhibits radially outward flow and is dominated by turbulent Maxwell stress. Winds are present but are not the dominant driver of accretion.

• •

The viscosity parameter $\alpha$ (incorporating both laminar and turbulent stress) is approximately 10 in the elevated accreting regions of the strong magnetic field case. In the disk midplane, we find an $\alpha\approx 0.5$ for the strong field case. The intermediate and weak magnetic field cases give an $\alpha$ value close to unity in the elevated regions and $\alpha<1$ close to disk midplane. A comparison of the inferred viscosity parameter $\alpha$ from our numerical simulations and previously performed local shearing box simulations (Salvesen et al., 2016) gives similar results.

• •

In addition to elevated accreting regions, the disk shows spiral structures, characterized by large density inhomogeneities, that can be modeled approximately as logarithmic spirals. The spiral structures are more open as the initial magnetic field strength increases.

• •

We find a quasi-steady state of the magnetic flux in all three simulations reported in this paper, suggesting that advection and diffusion of flux achieve a rough balance.

Our results, taken together with earlier studies, suggest that the qualitative structure of accretion flows is strongly influenced by the net magnetic flux threading the disk. The net flux may be almost as important as the accretion rate, or the efficiency of cooling, in determining observable properties of the flow. All of the three possible regimes of magnetization — effectively zero net flux, the magnetically elevated disks considered in this paper, and magnetically arrested disks — can plausibly be realized in astrophysically relevant situations. Magnetically elevated disks may be good candidates for systems where the basic geometry of accretion appears to be that of a geometrically thin disk, but where in detail the predictions of standard thin disk models fail. We have not constructed full physics models of magnetically elevated disks (which would have to include the physics of radiation pressure and radiative transport), but it is highly plausible that such disks would be thicker, with faster inflow speeds and enhanced thermal and viscous stability, compared to standard models. They would also have decreased self-gravity and highly inhomogeneous density structures. Qualitatively, these differences go in the right direction to address longstanding problems with standard accretion disk models.

We have presented evidence that the net poloidal fields necessary for sustaining magnetically elevated disks may attain a stable equilibrium (i.e., not run away toward either a zero net flux or MAD state). Nonetheless, how the required net flux is established at different radii in the disk remains an open question. The most straightforward source for this flux would be the material feeding the flow. However, it is not clear under what conditions a thin disk starting with a weaker net flux (say, $\beta_{0}>10^{4}$), would be able to overcome magnetic field diffusion according to the criterion proposed in Lubow et al. (1994). Alternatively, stochastic processes could lead to the emergence of patches of net flux in the form of large fluctuations (Beckwith et al., 2011; Begelman & Armitage, 2014). Longer duration numerical simulations, and simplified models for the evolution of the net flux, are needed to explore such issues.

Acknowledgements

BM thanks to Zhaohuan Zhu for disscussions at early stages of the project. We acknowledge computational support from the PROMETHEUS supercomputer in the PL-Grid infrastructure in Poland, and acknowledge support from NASA Astrophysics Theory Program grants NNX16AI40G and NNX17AK55G. PJA acknowledges support from NASA Astrophysics Theory Program grant 80NSSC18K0640.