Global 3D radiation magnetohydrodynamic simulations for FU Ori’s accretion disc and observational signatures of magnetic fields

ABSTRACT

FU Ori is the prototype of FU Orionis systems that are outbursting protoplanetary discs. Magnetic fields in FU Ori’s accretion discs have previously been detected using spectropolarimetry observations for Zeeman effects. We carry out global radiation ideal MHD simulations to study FU Ori’s inner accretion disc. We find that (1) when the disc is threaded by vertical magnetic fields, most accretion occurs in the magnetically dominated atmosphere at z ∼ R, similar to the ‘surface accretion’ mechanism in previous locally isothermal MHD simulations. (2) A moderate disc wind is launched in the vertical field simulations with a terminal speed of ∼300–500 km s⁻¹ and a mass-loss rate of 1–10 per cent the disc accretion rate, which is consistent with observations. Disc wind fails to be launched in simulations with net toroidal magnetic fields. (3) The disc photosphere at the unit optical depth can be either in the wind launching region or the accreting surface region. Magnetic fields have drastically different directions and magnitudes between these two regions. Our fiducial model agrees with previous optical Zeeman observations regarding both the field directions and magnitudes. On the other hand, simulations indicate that future Zeeman observations at near-IR wavelengths or towards other FU Orionis systems may reveal very different magnetic field structures. (4) Due to energy loss by the disc wind, the disc photosphere temperature is lower than that predicted by the thin disc theory, and the previously inferred disc accretion rate may be lower than the real accretion rate by a factor of ∼2–3.

accretion, accretion discs, astroparticle physics, dynamo, instabilities, MHD, turbulence

1 INTRODUCTION

Accretion discs have been observed in a wide range of astrophysical systems, ranging from around low mass stars (Hartmann, Herczeg & Calvet 2016) to around compact objects and supermassive black holes (Begelman, Blandford & Rees 1984). The accretion process not only helps to build the central object, but the released radiation energy allows us to identify and study the central object (e.g. X-ray binaries). The high-resolution M87 image by the Event Horizon Telescope (Event Horizon Telescope Collaboration et al. 2019) is an excellent example that we can constrain the properties of black holes by studying their surrounding accretion discs.

The leading theory to explain the accretion process involves magnetic fields, especially for sufficiently ionized discs.¹ Magnetic fields can drive turbulence through the magnetorotational instability (MRI; Balbus & Hawley 1991, 1998) or/and launch disc winds through the magnetocentrifugal effect in non-relativistic discs (Blandford & Payne 1982). The strengths of both MRI turbulence and disc winds depend on the field strength. Normally turbulence and wind are more prominent in systems having stronger magnetic fields (Hawley, Gammie & Balbus 1995).

Despite the importance of magnetic fields, the observational evidences for magnetic fields in accretion discs remain to be scarce. The collimated jets/outflows provide some indirect evidences of magnetic fields since the confinement of jets may require the presence of magnetic fields (Pudritz et al. 2007; Frank et al. 2014). Another indirect evidence is from magnetic field measurements from meteorites. Paleomagnetic measurements by Fu et al. (2014) suggest that Semarkona meteorites were magnetized to 0.54 G in the solar nebulae.

The most direct evidence of magnetic fields in accretion discs comes from Zeeman splitting of atomic or molecular lines. Current Zeeman measurements of molecular lines using ALMA (Vlemmings et al. 2019) have only placed upper limits on the field strength (<30 mG). So far, the only direct measurement of magnetic fields in accretion discs is the detection of Zeeman splitting of atomic lines coming from the inner disc of FU Ori (Donati et al. 2005).

FU Ori is the prototype of FU Orionis systems: a small but remarkable class of variable young stellar objects that undergo outbursts in optical light of 5 mag or more (Herbig 1977). While the outburst has a fast rise time (≲ 1–10 yr), the decay time-scale ranges from decades to centuries (Audard et al. 2014; Connelley & Reipurth 2018). Although more FU Orionis outbursts have been discovered recently thanks to large-scale all-sky surveys (e.g. Semkov et al. 2010; Kraus et al. 2016; Kóspál et al. 2017; Hillenbrand et al. 2018), the occurrence rate of these objects among young stars is still illusive (Scholz, Froebrich & Wood 2013; Hillenbrand & Findeisen 2015) with rates ranging from less than 1 outburst per young star to more than tens of outbursts per young star.

Such intense outbursts are due to the sudden increase of the protostellar disc’s accretion rate from |$\sim 10^{-8}\rm M_{\odot }\, yr^{-1}$| (Class I–II rates) to |$\sim 10^{-4}\rm M_{\odot }\, yr^{-1}$| (Hartmann & Kenyon 1996). The strong accretion is accompanied by the strong disc wind (Calvet, Hartmann & Kenyon 1993; Milliner et al. 2019). Although the outburst triggering mechanism is not clear,² the inner discs (≲1 au) during the outbursts are hot enough (∼6000 K, Zhu et al. 2007) to be sufficiently ionized and MRI should operate in these discs. Since these inner discs with |$\sim 100\rm L_{\odot }$| are much brighter than the central stars and all the light we see are from these accretion discs, FU Orionis systems are ideal places to study accretion physics.

Taking advantage of many atomic lines available in these systems, Donati et al. (2005) have used the high-resolution spectropolarimeter to detect signals of Zeeman splitting in FU Ori. By splitting the circular polarization signal into symmetric and antisymmetric components, they constrain the magnetic fields in both the azimuthal and radial directions. Assuming that the disc’s rotational axis is 60° inclined with respect to our line of sight, their best-fitting model suggests that the vertical component of the fields is ∼ 1 kG at 0.05 au and points towards the observer, while the azimuthal component (about half as strong) points in a direction opposite to the orbital rotation.

In spite of these stringent observational constraints, theoretical work still lacks behind and its connection with observations has not been established. To study FU Ori using theoretical numerical simulations, high enough numerical resolution is necessary for capturing MRI, while a large simulation domain is needed to study the disc wind. Only recently, with the newly developed Athena++ code that has both mesh-refinement and the special polar boundary condition, we can simulate the whole 4|$\pi$| sphere around the central object with enough resolution to capture MRI (Zhu & Stone 2018). Besides magnetic fields, radiative transfer is also crucial for understanding FU Ori’s inner accretion disc. For example, thermal instability was previously suggested to explain FU Ori’s outburst (Bell & Lin 1994). Although local shearing box MHD simulations with radiative transfer (Hirose et al. 2014) do not support the thermal instability theory for FU Ori outbursts (Hirose 2015), the disc’s thermal structure is still important for both the accretion physics (Zhu et al. 2009b) and the boundary layer physics (Kley & Lin 1999). Furthermore, radiative transfer is important for making connections with observations (e.g. understanding the physical condition at the disc’s photosphere).

Thus, in this work, we include radiative transfer in the global MHD disc simulations to study the accretion structure of FU Ori’s inner disc. We will also compare our simulations with previous Zeeman magnetic field observations and disc wind observations. In Section 2, the theoretical framework for energy transport in accretion discs is presented. We will describe our numerical method in Section 3. The results are presented in Section 4. After connecting with observations and a short discussion in Section 5, the paper is concluded in Section 6.

2 THEORETICAL FRAMEWORK

Angular momentum transport and energy transport are two important aspects of the accretion disc theory. Angular momentum transport is essential for the mass build-up of the central object, while energy transport is crucial for revealing disc properties using observations. Previously in Zhu & Stone (2018), we have done detailed analyses on angular momentum transport for discs threaded by net vertical magnetic fields. In this work, we will focus on energy transport in accretion discs.

The fluid equations with both magnetic and radiation fields are

$$\begin{eqnarray*} \frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho {\boldsymbol v}\right)&=&0\nonumber \\ \frac{\partial \rho {\boldsymbol v}}{\partial t}+\nabla \cdot \left(\rho {\boldsymbol v}{\boldsymbol v}-{\boldsymbol B}{\boldsymbol B}+{\sf P^*}+{\sf \sigma }\right)&=&-{\boldsymbol S_{r}}({\boldsymbol P})+{\boldsymbol F}\nonumber \\ \frac{\partial E}{\partial t}+\nabla \cdot \left[\left(E+P^{*}\right){\boldsymbol v}-{\boldsymbol B}\left({\boldsymbol B}\cdot {\boldsymbol v}\right)+{\sf \sigma }\cdot {{\boldsymbol v}}\right]&=&-cS_{r}(E)+{\boldsymbol F}\cdot {\boldsymbol v}\nonumber \\ \frac{\partial {\boldsymbol B}}{\partial t}-\nabla \times \left({\boldsymbol v}\times {\boldsymbol B}\right)&=&0, \end{eqnarray*}$$

(1)

where E = E_g + ρv²/2 + B²/2 is the total gas energy density, E_g = P/(γ − 1) is the internal energy, |${\sf P^*}\equiv (P+B^2/2){\sf I}$| is the pressure tensor (with |${\sf I}$| the unit tensor), and F is the external force (e.g. gravity). We also include the dissipation tensor |${\sf \sigma }$| in the equations. Although dissipation is not explicitly added in the simulations, shock dissipation is implicitly included in the Riemann solver, and dissipation terms are important for the energy analysis. The radiation equations are

$$\begin{eqnarray*} \frac{\partial E_{r}}{\partial t}+\nabla \cdot {\boldsymbol F_{r}}=cS_{r}(E) \end{eqnarray*}$$

(2)

$$\begin{eqnarray*} \frac{1}{c^2}\frac{\partial {\boldsymbol F_r}}{\partial t}+\nabla \cdot {\sf P_{r}}={\boldsymbol S_{r}({\boldsymbol P})}, \end{eqnarray*}$$

(3)

where the radiation flux F_r and the radiation energy density E_r are Eulerian variables, and they are related to the co-moving flux F_r,0 through |${\boldsymbol F_{r,0}}={\boldsymbol F_{r}}-({\boldsymbol v}E_{r}+{\boldsymbol v}\cdot {\sf P_r})$|⁠. The radiation pressure tensor |${\sf P_r}$| is related to the energy density though the variable Eddington tensor |${\sf P_r}={\sf f}E_{r}$|⁠. The source terms cS_r(E) and S_r(P) are given in Jiang, Stone & Davis (2013).

To study the energy budget, it is also helpful to write the equation for the gas’ internal energy density. The kinetic and magnetic energy equation is

$$\begin{eqnarray*} &&{\frac{\partial }{\partial t}\left(\frac{\rho v^2}{2}+\frac{B^2}{2}\right)}\nonumber \\ &&{\quad+\,\nabla \cdot \left[{\boldsymbol v} \left(\frac{\rho v^2}{2}\right)-{\boldsymbol B}\left({\boldsymbol B}\cdot {\boldsymbol v}\right)+({\sf P^*}+{\sf \sigma })\cdot {\boldsymbol v}\right]}\nonumber \\ &&{\quad-\left(P-\frac{B^2}{2}\right)\nabla \cdot {\bf v}+ \left({\boldsymbol v}\cdot \nabla \right)\frac{B^2}{2}-\left({\sf \sigma \cdot \nabla }\right)\cdot {\boldsymbol v}}\nonumber \\ &&{\quad=-{\boldsymbol v}\cdot {\boldsymbol S_{r}({\boldsymbol P})}+{\boldsymbol F}\cdot {\boldsymbol v}, } \end{eqnarray*}$$

(4)

so that the internal energy density is

$$\begin{eqnarray*} \frac{\partial E_g}{\partial t}+\nabla \cdot \left(E_g{\boldsymbol v}\right)+P\nabla \cdot {{\boldsymbol v}}+\left({\sf \sigma }\cdot \nabla \right)\cdot {\boldsymbol v}=-cS_r(E)+{\boldsymbol v}\cdot {\boldsymbol S_{r}(P)}, \nonumber\\ \end{eqnarray*}$$

(5)

which suggests that the change of the internal energy is due to the Pdv work, the dissipation, and radiative transport.

We can use either the equation for the total energy (equation 1) or the equation for the internal energy (equation 5) to derive the disc luminosity. Here, we rewrite the total energy equation as

$$\begin{eqnarray*} \frac{\partial E}{\partial t}+\nabla \cdot {\boldsymbol A}=-Q_{\rm cool}+ {\boldsymbol F}\cdot {\boldsymbol v}, \end{eqnarray*}$$

(6)

where A = (E + P*)v − B(B · v), and Q_cool is the radiative cooling rate. A can also be rewritten as

$$\begin{eqnarray*} {\boldsymbol A}=\left(\frac{\gamma }{\gamma -1}P+\frac{1}{2}\rho v^2\right){\boldsymbol v}+{\boldsymbol B}\times ({\boldsymbol v}\times {\boldsymbol B}) \end{eqnarray*}$$

(7)

using vector identities.

We will first review the thin disc theory under the cylindrical coordinate system and then we will write similar equations under the spherical-polar coordinate system that has been adopted in our simulations. The perturbed equation for the angular momentum under the cylindrical coordinate system can be written as

$$\begin{eqnarray*} \frac{\partial \langle \rho \delta v_{\phi }\rangle }{\partial t}&=&-\frac{1}{R^2}\frac{\partial (R^2 \langle T_{R\phi }\rangle)}{\partial R}-\frac{\langle \rho v_{R}\rangle }{R}\frac{\partial Rv_{K}}{\partial R}\nonumber \\ &&-\,\frac{\partial \langle T_{\phi z}\rangle }{\partial z}- \langle \rho v_{z}\rangle \frac{\partial v_{K}}{\partial z}, \end{eqnarray*}$$

(8)

where

$$\begin{eqnarray*} T_{R\phi }&\equiv& \rho v_{R}\delta v_{\phi }-B_{R}B_{\phi },\nonumber \\ T_{\phi z}&\equiv& \rho v_{z}\delta v_{\phi }-B_{z}B_{\phi }, \end{eqnarray*}$$

(9)

and 〈〉 denotes that the quantity has been averaged in the azimuthal (ϕ) direction. Assuming a steady state, we have

$$\begin{eqnarray*} \frac{\dot{M}}{2\pi }\frac{\partial Rv_{K}}{\partial R}=\frac{\partial (R^2 \langle T_{R\phi }\rangle)}{\partial R}+R^2\frac{\partial \langle T_{\phi z}\rangle }{\partial z}+R^2\langle \rho v_{z}\rangle \frac{\partial v_{K}}{\partial z}, \end{eqnarray*}$$

(10)

where |$\dot{M}\equiv -2\pi R\langle \rho v_{R}\rangle$|⁠. Thus, the accretion is driven by the T_Rϕ stress within the disc or the T_ϕz stress at the disc surface. If we assume that |$\dot{M}$| is a constant along R, we have

$$\begin{eqnarray*} \langle T_{R\phi }\rangle =\frac{\dot{M}v_{K}}{2\pi R}-\frac{C}{R^2}-\frac{1}{R^2}\int R^2\left(\frac{\partial \langle T_{\phi z}\rangle }{\partial z}+ \langle \rho v_{z}\rangle \frac{\partial v_{K}}{\partial z}\right){\rm d}R. \nonumber\\ \end{eqnarray*}$$

(11)

The energy equation (equation 6) under the cylindrical coordinate system is

$$\begin{eqnarray*} \frac{\partial \langle E\rangle }{\partial t} = -\frac{1}{R}\frac{\partial (R \langle A_{R}\rangle)}{\partial R}-\frac{\partial \langle A_{z}\rangle }{\partial z}-\langle Q_{\rm cool}\rangle +\langle {\boldsymbol F}\cdot {\boldsymbol V}\rangle , \end{eqnarray*}$$

(12)

where the leading terms in A_R (after removing the second-order terms) are

$$\begin{eqnarray*} A_{R}=\frac{\gamma }{\gamma -1}Pv_{R}+\frac{1}{2}\rho v_{R} v_{K}^2+ v_{K}T_{R\phi }, \end{eqnarray*}$$

(13)

and the leading terms in A_z are

$$\begin{eqnarray*} A_{z}=\frac{\gamma }{\gamma -1}Pv_{z}+\frac{1}{2}\rho v_{z} v_{K}^2+ v_{K}T_{\phi z}. \end{eqnarray*}$$

(14)

If we ignore the pressure term in A_R, assume v_z ∼ 0 in A_z, and assume a steady state, we have

$$\begin{eqnarray*} \langle Q_{cool}\rangle =-\frac{1}{R}\frac{\partial \left(\!\left\langle\! -\frac{1}{4\pi }\dot{M} v_{K}^2+ Rv_{K}T_{R\phi } \!\right\rangle\!\right)}{\partial R} -\frac{\partial \langle v_{K}T_{\phi z} \rangle }{\partial z} +\langle {\boldsymbol F}\cdot {\boldsymbol V}\rangle. \nonumber\\ \end{eqnarray*}$$

(15)

If we plug in T_Rϕ, ignore the T_ϕz term, replace F with the gravitational force, and only consider the disc mid-plane, we have

$$\begin{eqnarray*} 2\pi \langle Q_{cool}\rangle =-\frac{1}{2}\frac{\dot{M}v_{K}^2}{R^2}+\frac{\dot{M}v_{K}^2}{R^2}-\frac{3}{2}\frac{Cv_{K}}{R^3}+ \frac{\dot{M}v_{K}^2}{R^2}, \end{eqnarray*}$$

(16)

where the first term on the right is due to the radial derivative of the Keplerian kinetic energy flux, the second and third terms on the right are due to the radial derivative of the R − ϕ stress, and the last term on the right is the release of the gravitational potential energy. With the traditional zero stress inner boundary condition (⁠|$C=\dot{M}R_{in}v_{K,in}$|⁠), the cooling rate is

$$\begin{eqnarray*} \langle Q_{\rm cool}\rangle =\frac{3\dot{M}v_{K}^2}{4\pi R^2}\left(1-\left(\frac{R_{\rm in}}{R}\right)^{1/2}\right). \end{eqnarray*}$$

(17)

After the vertical integration, this cooling rate becomes what we normally use in the thin disc approximation,

$$\begin{eqnarray*} \sigma T_{\rm eff}^4=\frac{3G\dot{M}M}{8\pi R^3}\left(1-\left(\frac{R_{\rm in}}{R}\right)^{1/2}\right). \end{eqnarray*}$$

(18)

If we integrate over the whole disc starting from R_in, the total cooling rate is half the release rate of the gravitational potential energy (⁠|$GM\dot{M}/2R_{\rm in}$|⁠). On the other hand, far away from the central star (R ≫ R_in), the cooling rate (⁠|$3\dot{M}v_{K}^2/4\pi R^2$|⁠) is actually higher than the energy release rate from the gravitational contraction (⁠|$\dot{M}v_{K}^2/2\pi R^2$|⁠). The additional |$\dot{M}v_{K}^2/4\pi R^2$| energy release is due to the energy transport in the radial direction. We note that the same equation can also be derived using the internal energy equation but with an additional step to derive the dissipation term.

On the other hand, our simulated discs are very thick, and the disc photosphere flares roughly following the radial direction in the spherical grids. Thus, we want to derive similar equations for the spherical-polar coordinate system so that we can study energy transport in our simulations. The perturbed angular momentum equation under the spherical-polar coordinate system is

$$\begin{eqnarray*} \frac{\partial \langle \rho \delta v_{\phi }\rangle }{\partial t}&=&-\frac{1}{r^3}\frac{\partial (r^3 \langle T_{r\phi }\rangle)}{\partial r}-\frac{\langle \rho v_{r}\rangle }{r}\frac{\partial rv_{K}}{\partial r}\nonumber \\ &&-\,\frac{1}{r {\rm sin}^2\theta }\frac{\partial ({\rm sin}^2\theta \langle T_{\theta \phi }\rangle)}{\partial \theta }- \frac{\langle \rho v_{\theta }\rangle }{r {\rm sin}\theta }\frac{\partial ({\rm sin}\theta v_{K})}{\partial \theta }, \end{eqnarray*}$$

(19)

where

$$\begin{eqnarray*} T_{r\phi }&\equiv& \rho v_{r}\delta v_{\phi }-B_{r}B_{\phi },\nonumber \\ T_{\theta \phi }&\equiv& \rho v_{\theta }\delta v_{\phi }-B_{\theta }B_{\phi }. \end{eqnarray*}$$

(20)

Assuming a steady state, we have

$$\begin{eqnarray*} \dot{\widetilde{M}}\frac{\partial rv_{K}}{\partial r}&=&\frac{\partial (r^3 \langle T_{r\phi }\rangle)}{\partial r}+\frac{r^2}{{\rm sin}^2\theta }\frac{\partial ({\rm sin}^2\theta \langle T_{\theta \phi }\rangle)}{\partial \theta }\nonumber \\ &&+\,\frac{r^2\langle \rho v_{\theta }\rangle }{ {\rm sin}\theta }\frac{\partial ({\rm sin}\theta v_{K})},{\partial \theta } \end{eqnarray*}$$

(21)

where |$\dot{\widetilde{M}}=-r^2\langle \rho v_{r}\rangle$|⁠. Note that this |$\dot{\widetilde{M}}$| definition is different from the |$\dot{M}$| definition in the cylindrical coordinate system. If we assume that |$\dot{\widetilde{M}}$| is a constant along r, we can integrate the equation to derive

$$\begin{eqnarray*} \langle T_{r\phi } \rangle &=& \frac{\dot{\widetilde{M}}v_{K}}{r^2}-\frac{C}{r^3}-\frac{1}{r^3}\int \frac{r^2}{{\rm sin}^2\theta }\\ &&\times \left(\frac{\partial ({\rm sin}^2\theta \langle T_{\theta \phi }\rangle)}{\partial \theta }+{\rm sin}\theta \langle \rho v_{\theta }\rangle \frac{\partial ({\rm sin}\theta v_{K})}{\partial \theta }\right) {\rm d}r \end{eqnarray*}$$

(22),(23)

The energy equation (equation 6) under the spherical-polar coordinate system is

$$\begin{eqnarray*} \frac{\partial \langle E\rangle }{\partial t} \!=\! -\frac{1}{r^2}\frac{\partial (r^2 \langle A_{r}\rangle)}{\partial r}-\frac{1}{r {\rm sin}\theta }\frac{\partial ({\rm sin}\theta \langle A_{\theta }\rangle)}{\partial \theta }-\langle Q_{cool}\rangle \!+\!\langle {\boldsymbol F}\cdot {\boldsymbol V}\rangle . \nonumber\\ \end{eqnarray*}$$

(24)

The leading terms in A_r are

$$\begin{eqnarray*} A_{r}=\left(\frac{\gamma }{\gamma -1}P+\frac{1}{2}\rho v_{K}^2+\rho v_{K}\delta v_{\phi }\right)v_{r}-v_{K}B_{r}B_{\phi } \end{eqnarray*}$$

(25)

or

$$\begin{eqnarray*} A_{r}=\frac{\gamma }{\gamma -1}Pv_{r}+\frac{1}{2}\rho v_{r} v_{K}^2+ v_{K}T_{r\phi } \end{eqnarray*}$$

(26)

The leading terms in A_θ are

$$\begin{eqnarray*} A_{\theta }=\frac{\gamma }{\gamma -1}Pv_{\theta }+\frac{1}{2}\rho v_{\theta } v_{K}^2+ v_{K}T_{\theta \phi }. \end{eqnarray*}$$

(27)

In Section 4.2, we will measure the energy transport due to the A_r and A_θ terms from our simulations. On the other hand, in this section, we will continue the derivation by making several assumptions. If we ignore the pressure term in A_r, assume v_θ ∼ 0 in A_θ, and assume a steady state, we have

$$\begin{eqnarray*} \langle Q_{cool}\rangle &=&-\frac{1}{r^2}\frac{\partial \left(\left\langle -\frac{1}{2}\dot{\widetilde{M}} v_{K}^2+ r^2 v_{K}T_{r\phi } \right\rangle\right)}{\partial r}\nonumber \\ &&-\,\frac{1}{r {\rm sin}\theta }\frac{\partial ({\rm sin}\theta \langle v_{K}T_{\theta \phi } \rangle)}{\partial \theta }+\langle {\boldsymbol F}\cdot {\boldsymbol V}\rangle . \end{eqnarray*}$$

(28)

If we plug in T_rϕ from equation (23) and ignore T_θϕ terms, we have

$$\begin{eqnarray*} \langle Q_{cool}\rangle =-\frac{1}{2}\frac{\dot{\widetilde{M}}v_{K}^2}{r^3}+\frac{\dot{\widetilde{M}}v_{K}^2}{r^3}-\frac{3 }{2}\frac{C v_{K}}{ r^4}+\frac{\dot{\widetilde{M}}v_{K}^2}{r^3}. \end{eqnarray*}$$

(29)

Thus, if we can ignore the θ direction energy advection/stress and the boundary C term, the cooling rate equals the release rate of the gravitational potential energy (the last term on the right) plus the radially advected energy (the first two terms on the right). Unfortunately, as will be shown in Section 4.2, the energy advection in the θ direction and the T_θϕ stress cannot be ignored. Accordingly, the cooling rate is modified significantly.

3 METHOD

We solve the magnetohydrodynamic (MHD) equations in the ideal MHD limit using Athena++ (Stone et al. 2020, in press). Athena++ is a newly developed grid-based code using a higher order Godunov scheme for MHD and the constrained transport (CT) to conserve the divergence-free property for magnetic fields. Compared with its predecessor Athena (Gardiner & Stone 2005, 2008; Stone et al. 2008), Athena++ is highly optimized for speed and uses a flexible grid structure that enables mesh refinement, allowing global numerical simulations spanning a large radial range. Furthermore, the geometric source terms in curvilinear coordinates (e.g. in cylindrical and spherical-polar coordinates) are specifically implemented to converse the angular momentum to machine precision. In this work, we adopt the second-order piecewise-linear method for the spatial reconstruction, the second-order Van-Leer method for the time integration, and the HLLC Riemann solver to calculate the flux.

The time-dependent radiative transfer equation has been solved explicitly and coupled with the MHD fluid equations using the radiation module of Jiang, Stone & Davis (2014a). The general radiative transfer equation for the static fluid is

$$\begin{eqnarray*} \frac{1}{c}\frac{\partial I_{\nu }}{\partial t}+ {\boldsymbol n}\cdot \nabla I_{\nu }=-(\sigma _{\nu ,a}+\sigma _{\nu ,s})I_{\nu }+j_{\nu }+\sigma _{\nu ,s}^{\rm eff}J_{\nu } \end{eqnarray*}$$

(30)

where I_ν(x, t, n) is the intensity at the position x, time t and along the direction of n. J_ν = (4|$\pi$|⁠)⁻¹∫I_νdΩ and j_ν/σ_ν,a = B_ν, while σ_ν,a and σ_ν,s are the absorption and scattering opacity at the frequency of ν. However, for a fluid that is moving at v, additional correction terms of the order of (v/c) and (v/c)² need to be added (Jiang et al. 2014a). Jiang, Stone & Davis (2019a) has adopted a mixed frame approach to solve the radiative transfer equation for moving fluid consistently. After integrating the radiative transfer equation over frequency, the equation becomes

$$\begin{eqnarray*} \frac{1}{c}\frac{\partial I}{\partial t}+ {\boldsymbol n}\cdot \nabla I= S(I,{\boldsymbol n}). \end{eqnarray*}$$

(31)

After carrying out the transport step in the lab frame, the source terms on the right-hand side are added. But instead of adding the source term S(I, n) with all the (v/c) and (v/c)² corrections to the intensity, the lab frame specific intensity I(n) at angle n is first transformed to the comoving frame intensity I₀(n₀) via Lorentz transformation. Then the source terms in the comoving frame (S₀(I₀, n₀)) are added to I₀(n₀),

$$\begin{eqnarray*} S_0(I_0,{\boldsymbol n_0})&=&\sigma _{a,R}\left(\frac{a_r T^4}{4\pi }-I_0\right)+\sigma _s \left(J_0-I_0\right)\nonumber \\ &&+\,\left(\sigma _{a,P}-\sigma _{a,R}\right)\left(\frac{a_rT^4}{4\pi }-J_0\right), \end{eqnarray*}$$

(32)

where σ_a,R = κ_a,R × ρ, and σ_a,P = κ_a,P × ρ. κ_a,R and κ_a,P are the Rosseland mean and Planck mean opacities. After this step to update I₀(n₀), I₀(n₀) are transformed back to the lab frame via Lorentz transformation. Then, the radiation momentum and energy source terms that are used in the fluid equations are calculated by the differences between the angular quadratures of I(n) in the lab frame before and after adding the source terms.

For our particular FU Ori problem, we find that using the higher order PPM scheme (Colella & Woodward 1984) for the transport step is crucial for deriving the correct radiation fields in the extremely optically thick regime (see Section 3.2). Thus, the PPM scheme has been used in all our simulations for solving the radiative transfer equation. Since the characteristic speed in the transport step is the speed of light, solving this equation explicitly requires very small numerical time-steps. Thus, we adopt the reduced speed of light approach as in Zhang et al. (2018). We reduce the speed of light by a factor of 1000, which still achieves a good time–scale separation between the radiation transport and fluid dynamics. More discussions and tests on the reduced speed of light approach are in Section 3.2. We solve the radiative transfer equation along 80 rays in different directions. Integration of the specific intensity over angles yields various radiation quantities and source terms for the fluid equations.

The opacity that is adopted in the radiative transfer equation is generated in Zhu et al. (2007) and Zhu et al. (2009a). With this opacity, Zhu et al. (2007) find an excellent agreement between the synthetic spectral energy distributions and observations for FU Ori. This gives us great confidence to adopt it in this work for FU Ori MHD simulations. Both Rosseland mean and Planck mean opacities are shown in Fig. 1. The dust opacity that is below ∼1500 K is derived based on the prescription in D’Alessio, Calvet & Hartmann (2001). The molecular, atomic, and ionized gas opacities have been calculated using the Opacity Distribution Function (ODF) method (Castelli & Kurucz 2004; Sbordone et al. 2004; Kurucz 2005) that is a statistical approach to handling line blanketing when millions of lines are present in a small wavelength range (Kurucz, Peytremann & Avrett 1974 ). More details on these opacities can be found in Zhu et al. (2009a) and Keith & Wardle (2014). On the other hand, we adopt a simple equation of state with a constant γ = 5/3 and μ = 1 to avoid any complications due to the change of γ and μ with the temperature.

Figure 1.

The Rosseland mean (solid black curves) and Planck mean (red dashed curves) opacities adopted in the simulations. Different curves represent opacities under different pressures (10⁻³ to |$10^5\, \rm dyn\ cm^{-2}$|⁠). Curves with overall lower values correspond to lower pressures.

Our grid set-up is similar to Zhu & Stone (2018), where the whole 4|$\pi$| sphere is covered by the spherical-polar (r, θ, ϕ) grids with the special polar boundary condition in the θ direction (details in the appendix of Zhu & Stone 2018). The grid is uniformly spaced in ln(r), θ, ϕ with 128 × 64 ×64 grid cells in the domain of [ln(0.25), ln(100)]× [0, |$\pi$|] ×[0, 2|$\pi$|] at the root level. Two levels of mesh refinement have been adopted at the disc mid-plane with θ = [|$\pi$|/4, 3|$\pi$|/8] and [5|$\pi$|/8, 3|$\pi$|/4] for the first level and θ = [3|$\pi$|/8, 5|$\pi$|/8] for the second level. The outflow boundary conditions for flow variables, magnetic fields, and radiation fields have been adopted at both the inner and outer radial boundaries. Additionally, v_r at the radial boundaries is set to prevent the inflow to the simulation domain.

The disc’s initial density, temperature, and velocity profiles are also similar to Zhu & Stone (2018) but with the mid-plane density slope of p = −2.125, the temperature slope of q = −3/4, and H/R = 0.2 at R = 1. This structure is consistent with the structure of a viscously heated α disc. The initial disc scale height is thus resolved by 16 grids with two levels of mesh refinement. The density floor is also similar to Zhu & Stone (2018) except that an additional factor of r_min/r was multiplied to equation (10) of Zhu & Stone (2018) to further decrease the floor value at the disc atmosphere.

Simulations with both net vertical and net toroidal magnetic fields have been carried out. The net vertical field set-up is similar to that in Zhu & Stone (2018) with a constant plasma β at the disc mid-plane initially. In the net toroidal field simulations, magnetic fields are only present within 2 disc scale heights above and below the mid-plane initially, and the plasma β is a constant anywhere within this region.

3.1 FU Ori parameters and simulation runs

Our simulations adopt the disc parameters that are consistent with FU Ori observations. The detailed disc atmospheric modelling (Zhu et al. 2007, 2008) suggests that FU Ori’s inner accretion disc extends from 5 R_⊙ to ∼ 1 au with an accretion rate of |$2.4\times 10^{-4}\rm M_{\odot }\, yr^{-1}$|⁠. The mass of the central star is 0.3 M_⊙. The rotational axis of the disc is 55° inclined with respect to our line of sight. Although these derived parameters are subject to change due to the recent Gaia distance measurement and ALMA disc inclination measurement for FU Ori (see Section 5.3), we will use these numbers as a guidance for our simulation parameters.

The length unit (R = 1) in the simulation is chosen as 0.1 au so that the whole domain extends from 5 R_⊙ to 10 au. The density unit is chosen as 10⁻⁸ g cm⁻³ and the initial mid-plane density is 10⁻⁷ g cm⁻³ at 0.1 au. The time unit is chosen as 1/Ω at 0.1 au around a 0.3 M_⊙ star. In this paper, we use T₀ to represent the orbital period (2|$\pi$|/Ω) at 0.1 au around a 0.3 M_⊙ star, which is 21 d. With these units, the initial disc surface density is Σ₀(r) = 7.5 × 10⁴(R/0.1 au)⁻¹g cm⁻².

Three main simulations have been carried out (1) the disc that is initially threaded by net vertical magnetic fields with the strength of β₀ = 1000 at the disc mid-plane, labelled as V1000, (2) the disc that is threaded by vertical fields with β₀ = 10⁴, labelled as V1e4, and (3) the disc that is initially threaded by net toroidal fields with the strength of β₀ = 100, labelled as T100. We run these simulations to T ∼ 60 T₀, which is equivalent to ∼3 yr. This time is equivalent to 500 innermost orbits in the simulation, and the disc at R = 1 has reached to a steady state as shown below.

3.2 Code tests

Although the radiative transfer scheme has been tested extensively (e.g. Jiang et al. 2014a, 2019a), we still need to test if the scheme is applicable to our particular FU Ori disc set-up. Thus, we set up a 1D plane-parallel atmosphere with a density profile of

$$\begin{eqnarray*} \rho =\rho _{0} {\rm e}^{-z^2/2H^2}, \end{eqnarray*}$$

(33)

to represent the disc’s vertical density structure at R = 1 in our 3D FU Ori simulations. H is chosen as 0.02 au, and ρ₀ is chosen as 10⁻⁸ g cm⁻³. All other parameters are the same as our 3D FU Ori simulations. To maintain this density structure, we do not update the density and velocity during the simulation, and only allow the disc temperature to change. Only two rays have been used in this set-up so that we can use two-stream approximation to calculate the analytical solution.

To represent the viscous heating in the accretion disc, we manually include a heating source term with the heating rate that is proportional to the disc local density (ρ) as

$$\begin{eqnarray*} \frac{{\rm d}E}{{\rm d}t}=C\times \rho . \end{eqnarray*}$$

(34)

We have done three tests, two of which are steady-state tests with a constant C and one of which is the increasing heat test where C suddenly increases at some time.

In the steady state tests, two different values of C (0.0002316 and 0.02316 in the code unit) have been used to test if the disc can reach to the correct temperature in both low and high-temperature regimes. The lower heating rate only heats the disc to T ∼ 10³ K, when the opacity is dominated by the dust and molecular opacities (the upper panels in Fig. 2). The higher heating rate heats the disc to T ∼ 10⁴ K, when the opacity is dominated by the free–free and bound-free opacities (the lower panels in Fig. 2).

Figure 2.

Plane-parallel atmosphere tests for atmospheres having two different heating rates (the simulation with the lower heating rate is shown in the upper panels). Density, temperature, and Rosseland mean opacity at t = 1000T₀ are shown from the left to right panels. The black crosses and curves are results from low resolution simulations while the red curves are from the simulations with 10 times higher resolution. The blue dotted curves in the middle panels show the analytical temperature profiles.

These steady-state tests show that we can accurately simulate the disc thermal structure for some cases, but also reveal the limitation of our set-up. The black crosses in Fig. 2 are results from simulations with 160 grids from −0.1 to 0.1 au (the same resolution as our 3D simulations), while the red curves are from simulations with 1600 grids in the same domain range. The blue curves in the middle panel are the analytical solutions of this problem solved with the two-stream approximation:

$$\begin{eqnarray*} T(\tau)^4=\frac{3}{4}T_{\rm eff}^4\left(\tau \left(1-\frac{\tau }{\tau _{\rm tot}}\right)+\sqrt{\frac{1}{3}}\right), \end{eqnarray*}$$

(35)

where |$\sigma T_{\rm eff}^4$| is the flux emerging from one side of the disc and τ_tot is the total optical depth from both sides of the disc. Clearly, when the opacity is low (e.g. the upper panels), the simulations with different resolutions agree with the analytical solution very well, even if the opacity has sharp changes among grids. On the other hand, when the opacity is high (e.g. the bottom panels), the optical depth can jump more than one order of magnitude from one grid to another grid. As expected, this jump leads to large errors in the calculations. Unfortunately, even with the resolution that is 10 times higher (red curves in the lower panels), we still cannot recover the analytical solution accurately. One way to overcome this problem in future is using adaptive mesh-refinement for those grid cells having high optical depths. Overall, this test shows that, with our current set-up, we may underestimate the temperature of some extremely optically thick grid cells by a factor of 2.

Since FU Ori’s disc temperature can change dramatically before and during the outburst, we also need to test if the code can capture the time evolution of the disc’s temperature accurately. Especially, our adoption of the reduced speed of light approach may delay the escape of the radiation energy. This is a particular concern when the disc is very optically thick (Skinner & Ostriker 2013), since the diffusion time-scale Lτ/c can now be longer than the dynamical time-scale. For a typical size scale of 0.1 au and an optical depth of 1000, the radiation diffusion time-scale is ∼1 d. Naively, we would think that decreasing the speed of light by 1000 will increase the diffusion time-scale to 1000 d, which is even longer than the total simulation time-scale. On the other hand, it can be shown that the formulation in Zhang et al. (2018) guarantees that the radiative diffusion flux is the correct flux when the thermal energy of the gas dominates over the radiation energy. Thus, we should expect a correct diffusion time-scale for our set-up where the thermal energy of the gas always dominates. However, one could also argue that the optically thick region is joined by the optically thin region, and the escape of the total energy will be controlled by the optically thin region so that the disc will still cool/heat slower with the reduced speed of light approach.

To resolve these concerns, we carry out a test with a suddenly increased heating rate. We fix the absorption opacity to be 0.1 cm² g⁻¹ in this test. Initially, the disc is heated at the heating rate of C = 0.000 2316 for a period of 2 T₀ so that the disc reaches to a steady state. Then, we suddenly increase the heating rate by a factor of 100 and watch the subsequent disc evolution. As shown in Fig. 3, the reduced speed of light approach indeed slows down the heating of the disc. On the other hand, the temperature structure at 0.5 T₀ after the heating event for the disc using the reduced speed of light approach (the black dashed curve) overlaps with the temperature structure at 0.1 T₀ after the heating event for the disc using the normal speed of light (the red dotted curve). Thus, the reduced speed of light approach increases the diffusion time-scale by a factor of ∼ 5. This is larger than 1, but it is also much smaller than 1000 so that the diffusion time-scale is still much shorter than the simulation time-scale. Nevertheless, since the reduced speed of light approach increases the diffusion time-scale to ∼T₀, we cannot trust short time-scale variations of the radiation field in the simulations, and we can only study the state when the disc is relatively steady for the orbital time-scale. Thus, in this paper, we only focus on the disc at the steady state with a constant accretion rate instead of discussing the outburst stage when the disc suddenly brightens by orders of magnitude within a short period of time.

Figure 3.

Plane-parallel atmosphere tests similar to Fig. 2 but with a sudden increase of the heating rate. With a normal heating rate, the disc reaches to a steady state after 2 T₀ (the solid black and red curves). Then, the heating rate suddenly jumps to a value that is 100 times higher. After another 0.1 T₀, the disc thermal structures are shown as the dotted curves. Then, after another 0.4 T₀, the disc thermal structures are shown as dashed curves. The adopted absorption opacity is 0.1 cm² g⁻¹. Clearly, using the reduced speed of light approach increases the time-scale of radiation escaping the disc.

4 RESULTS

The temperature and density structures of our fiducial model (V1000) at 50T₀ are shown in Fig. 4. We can see that the disc atmosphere at z ∼ R still has a significant density, which is similar to the disc structure in Zhu & Stone (2018). With the radiative transfer included in our simulations, we can now study the disc’s temperature structure. The disc’s temperature is quite high (≳5000 K) close to the central star (≲0.15 au). There is a sharp temperature jump around 0.15 au, indicating that the inner disc is at the upper branch of the equilibrium ‘S’ curve which is dominated by the bound-free and free–free opacity while the outer disc is at the lower branch of the equilibrium ‘S’ curve (≲2500 K) that is dominated by the molecular opacity. We also use ρκ_R × 0.1 au ∼10 to illustrate the disc’s photosphere. Clearly, the photosphere is hotter at the inner disc than at the outer disc, and the photosphere is not smooth having noticeable structures. Fig. 5 shows the density structure at the disc surface and the mid-plane. At the mid-plane, we clearly see spiral arms similar to those found in Mishra et al. (2020). On the other hand, the disc surface has filamentary structure due to surface accretion, as found in Zhu & Stone (2018) and Suriano et al. (2018). Due to these large-scale structures at the photosphere, we expect that FU Ori has short time-scale variations that have been implied by observations (Kenyon et al. 2000; Herbig, Petrov & Duemmler 2003; Powell et al. 2012; Siwak et al. 2013).

Figure 4.

The poloidal slice of the temperature (the upper half) and density (the lower half) from the V1000 case at 50 T₀. This illustrated region represents FU Ori disc within 0.5 au from the central star. For the upper half of the image, the disc’s photosphere is illustrated with the iso-surface having ρκ_R × 0.1 au = 10.

Figure 5.

The contours of log₁₀ρ at the θ = 1 plane (the left-hand panel) and the mid-plane (the right-hand panel) at 50 T₀. The colour range represents three orders of magnitude change of density in both panels.

After running for 50 T₀, our fiducial model has reached to a steady state within R ∼0.5 au, i.e. the inner factor of ∼ 20 in radius, as evident in Fig. 6. From the mass accretion rate panel (the upper right-hand panel), we can see that the region that is accreting inwards expands with time since the outer disc region takes more time for MRI to grow. At 50 T₀, the region within 0.5 au, i.e. the inner factor of ∼20 in radius, accretes inwards at a steady rate. Such constant accretion rates are also consistent with the stress profiles shown in the middle left panel. The vertically integrated Rϕ stress follows R^−1.5 and this leads to a constant accretion rate based on equation (10). Such accretion and stress structures are very similar to the global MHD simulations with the locally isothermal equation of state (compared with fig. 3 in Zhu & Stone 2018).

Figure 6.

The disc mid-plane density, surface density, mass accretion rate (upper panels), stresses (the solid curves are rϕ stresses at the mid-plane while the dashed curves are the vertically integrated Rϕ stresses), mid-plane α, vertically integrated α (middle panels), temperature, mid-plane Rosseland mean opacity, and 〈B²〉/2P_mid,0 (lower panels) at different times. α_total and α_int are calculated with the rϕ and Rϕ stresses, respectively. In the temperature and 〈B²〉/2P_mid,0 (where P_mid,0 is the mid-plane pressure from the initial condition) panels, the solid curves are the mid-plane quantities and the dashed curves are the quantities along r at θ = 0.78 (where the photosphere roughly sits). The black dotted line in the temperature panel is from equation (37) with an accretion rate of 4|$\times 10^{-4}\rm M_{\odot }\, yr^{-1}$| around a 0.3 M_⊙ star.

However, other quantities shown in Fig. 6 are drastically different from those in fig. 3 of Zhu & Stone (2018). For example, the surface density in Fig. 6 is almost flat, which is different from R^−0.6 in Zhu & Stone (2018). The mid-plane α is also flat compared with R^0.5 in Zhu & Stone (2018). Such differences are likely due to the temperature structure at the mid-plane. In the viscous heating dominated disc presented here, the mid-plane temperature follows ∼R^−3/4 (the lower left-hand panel), while, in the locally isothermal simulations, the mid-plane temperature follows R^−1/2. Another evidence that the mid-plane temperature affects the α profile is that, at R ∼ 0.15 au where the mid-plane temperature jumps down, the α_total,mid there jumps up so that the total stress T_total is still smooth. It is quite surprising that the accretion and stress profiles are smooth despite the jump of disc temperature. Considering that most stress is from the magnetic stress, this implies that the disc’s accretion structure is mainly controlled by the global geometry of magnetic fields and is insensitive to the disc local temperature. The magnetic fields at the mid-plane and θ = 0.78 are shown in the lower right-hand panel, and we can see that the field strength changes smoothly in the disc despite the temperature jump at R ∼0.15 au.

4.1 Accretion structure

The flow structure in MHD discs is tightly coupled with the magnetic field geometry. Magnetic fields determine the accretion structure, while the accretion process drags and alters the magnetic fields. We plot the azimuthally averaged temperature, density, and magnetic field structures for our fiducial run in Fig. 7.

Figure 7.

The azimuthally averaged temperature (the left-hand panel), density (the middle panel), and B_ϕ (the right-hand panel) for the V1000 case at 50 T₀. The green lines in the middle panel are the streamlines for the poloidal velocity fields, while the green lines in the right-hand panel are the streamlines for the poloidal magnetic fields (the direction of the magnetic fields at the upper boundary is pointing upwards). The white contours in all these panels are the β = 1 surfaces. The purple curves in the left-hand panel are the contours where T = 4000, 7000, and 10 000 K. The blue curves in the three panels are the τ_R = 1 surfaces. The dashed curves in the middle and right-hand panels are the Alfven surfaces.

The velocity and magnetic field structures are remarkably similar to the ‘surface accretion’ picture in locally isothermal discs with net vertical fields (Zhu & Stone 2018). Although we called such surface accretion as ‘coronal accretion’ in Zhu & Stone (2018) following Beckwith, Hawley & Krolik (2009), the accreting surface may not be as hot as Sun’s ‘corona’ that exceeds 10⁶ K (as shown in this work and Jiang et al. 2019b). On the other hand, the accreting surface is more associated with the strong magnetic fields (β ≲1, or called magnetically elevated in Mishra et al. 2020). Thus, in this work, we call this structure as ‘surface accretion’ instead. The flow structure can be separated into three regions from the mid-plane upwards: the disc region which is dominated by MRI turbulence, the surface accreting region which is above the β = 1 surface and extends all the way to z ∼ R, and the disc wind region (with v_r > 0) at z ≳ R. The accretion flow mainly occurs at the surface, as shown in the middle panel of Fig. 7 where the velocity streamlines are towards the star in the surface accreting region. Such surface inflow drags magnetic fields inwards so that the fields are pinched at the disc surface (the right-hand panel of Fig. 7). Due to the increase of the Keplerian rotation speed towards the inner disc, these dragged-in magnetic fields are sheared azimuthally, leading to fields with opposite B_ϕ between the lower and higher surface regions. Such surface accretion has been seen as early as Stone & Norman (1994) and recently in several simulations (Beckwith et al. 2009; Suzuki & Inutsuka 2009; Suriano et al. 2018; Takasao et al. 2018; Zhu & Stone 2018; Mishra et al. 2020; Jiang et al. 2019b). Analytical works by Guilet & Ogilvie (2012), Guilet & Ogilvie (2013) have also seen such surface accretion when the turbulent viscosity and diffusivity are considered in their analytical works.

On the other hand, our radiation MHD simulations reveal new information on the disc thermal structure, especially the position of the disc photosphere. The left-hand panel of Fig. 7 shows that the thermal radiation field is very smooth except at the sharp jump ∼0.15 au separating the two states that reside at the upper and lower branches of the ‘S’ curve. If we integrate the Rosseland mean opacity along the z-direction (starting from 20^o off the axis to avoid the coarse grids at the pole), the derived τ_R = 1 surface is plotted as the blue curves in all three panels. We can see that the τ_R = 1 surface is at the wind base or upper surface accreting region at the inner disc (≲0.07 au) and within the lower surface accreting region at the outer disc (≳0.07 au). Thus, B_ϕ derived from the atomic lines at the photosphere could have opposite directions depending on where these lines are produced. This has important implications for the B field measurements of FU Ori, which will be discussed in greater detail in Section 5.1. This transition radius ∼ 0.07 au, which roughly corresponds to the filamentary structure shown in Fig. 5, may also be related to the periodic variability at 10–15 d found in Herbig et al. (2003), Powell et al. (2012), Siwak et al. (2013), and Siwak et al. (2018).

To understand the disc’s accretion structure quantitatively, we plot the vertical profiles of various quantities at 0.1 au in Fig. 8. The yellow shaded region is the surface accreting region. We see that the density flattens out in the surface accreting region, and the radial accretion velocity can reach 20 km s⁻¹ there (the v_R panel). Considering that the Keplerian velocity is 50 km s⁻¹ at 0.1 au, the surface inflow velocity is ∼40 per cent of the Keplerian velocity. Due to the high speed, most disc mass is accreted through this surface accreting region despite its low density (the ρv_r panel). The azimuthal velocity also deviates from the Keplerian velocity. In the surface accreting region, the lowest azimuthal velocity can reach to |$60{{\ \rm per\ cent}}$| of the Keplerian velocity (the v_ϕ panel). Such low azimuthal velocity and high radial velocity can be understood as magnetic breaking by the mid-plane so that the surface loses angular momentum and falls inwards. The mid-plane is very hot with a high opacity. Here, at R = 0.1 au, the disc’s photosphere (τ_R = 1) is within the surface accreting region (the τ panel).

Figure 8.

Density, velocities, temperature, mass flux, opacity, and optical depth along the z direction at R = 0.1 au at 50 T₀. The quantities have been averaged azimuthally. The dashed curves in the velocity panels show the velocity components in the spherical-polar coordinates (V_r and −V_θ). The yellow curves are from the initial condition. The yellow shaded region labels the surface accreting region. Note the fast inward flow at the disc surface.

The magnetic field structure at R = 0.1 au is shown in Fig. 9. The surface inflow drags the initial vertical magnetic fields inwards, pinching the magnetic fields at the disc surface. The radial component of the magnetic fields in the surface accreting region has been sheared by the Keplerian rotation to produce a strong azimuthal component. The azimuthal B component can reach to 100 G, which is ∼5 times the radial B component. The combination of B_z and B_ϕ produces positive ∂T_ϕz/∂z at the base of the surface accreting region. Using equation (10), we can see that this T_ϕz leads to the inward accretion of the surface. In other words, the mid-plane is magnetically breaking the surface region. On the other hand, the internal T_ϕz stress will only transfer angular momentum from the surface to the disc mid-plane, and thus it won’t lead to the overall disc accretion. The overall disc accretion is led by the T_Rϕ stress within the disc and the T_ϕz stress at the disc atmosphere (e.g. the magnetocentrifugal wind). The detailed analysis on the surface accretion can be found in Zhu & Stone (2018). The accretion mechanisms are very similar. The only difference we notice by comparing Fig. 9 in this work with Fig. 7 in Zhu & Stone (2018) is that T_ϕz plays a more important role in FU Ori discs which are thicker than discs in Zhu & Stone (2018). We have verified that the radiation viscosity is not important here. It is at least 5 orders of magnitude lower than the magnetic stress, which is different from the sub-eddington accretion discs around supermassive black holes (Jiang et al. 2019b).

Figure 9.

Similar to Fig. 8 but for quantities that are related to magnetic fields. The dashed curves are B components in the spherical-polar coordinates (B_r and −B_θ). The blue curves in the T_Rϕ and T_ϕz panels are magnetic stresses that are calculated using the mean fields, −|$\overline{B_{R}}\times \overline{B_{\phi }}$| and −|$\overline{B_{z}}\times \overline{B_{\phi }}$|⁠. The mean fields are azimuthally averaged before being used to calculate the stress.

Although it is mainly the magnetic field that determines the accretion process, the radiation pressure in FU Ori plays a role in supporting the disc. The lower panel of Fig. 10 shows the force balance with various terms in the vertical momentum equation (equation 1). In a steady state, the stress tensor divergence and the vertical gradient of the total pressure are balanced by the vertical component of the gravitational force and the radiation pressure force. For a slowly moving fluid, the radiation pressure force is −S_r(P) = σ_tF_{r, 0}/c. Close to the disc midplane (the white region around z = 0), it is mainly the gradient of the gas pressure (the red curve) that balances the vertical forces (the black curves). The magnetic pressure gradient (the blue curve) has the same strength as the radiation pressure (the black dotted curve, |$\sim 30{{\ \rm per\ cent}}$| of the gas pressure), and thus they balance each other. The stress tensor also contributes to compressing the disc. In the surface accreting region, It is mainly the gradient of the magnetic pressure that balances the gravity. Both the radiation pressure and the gradient of the gas pressure are negligible at the surface in comparison. This again suggests that the surface accretion occurs in the magnetically dominated region.

Figure 10.

The upper panel: The vertical energy flux at R = 0.1 au due to radiation (F_{r, z}) and convection (<E_gv_z >). The bottom panel: The force balance in the vertical direction, including the gravitational force (F_gra), the radiation force (σ_tF_{r0, z}/c), the gas pressure gradient (dP/dz), and the magnetic pressure gradient (dP_mag/dz). All quantities are averaged over both the azimuthal direction and time (45 to 50 T₀ with a 0.1 T₀ interval).

4.2 Energy budget

Angular momentum transport and energy transport are the two most important aspects of accretion discs. In Zhu & Stone (2018), we have done analyses on the angular momentum budget of accretion discs threaded by net vertical magnetic fields. With the radiative transfer included in this work, we will do similar analyses for the disc’s energy budget. The formulas are laid out in Section 2. Since the energy budget is related to the angular momentum budget, we will first repeat the angular momentum analysis as we did in Zhu & Stone (2018).

The angular momentum budget is shown in the upper panel of Fig. 11. Four different terms in the angular momentum equation (equation 19) are plotted. The m_rϕ term is the radial gradient of the r–ϕ stress (the first term on the right-hand side of equation 19). After the integration over a volume in the disc, this term represents the transport due to the internal stress exerted at the face that is perpendicular to the disc mid-plane, either from the turbulent stress or the stress due to the large scale organized magnetic fields. The m_θϕ term is the θ gradient of the θ–ϕ stress (the third term on the right-hand side of equation 19). After the integration over a volume, it is the stress that is exerted at the disc surface. That is normally due to the magnetocentrifugal disc wind. The other two terms (the |$\dot{m}_r$| term, which is the second term on the right-hand side of equation 19, and the |$\dot{m}_\theta$| term, which is the fourth term on the right-hand side of equation 19) are the momentum transport due to the radial and poloidal mass flux. In the thin disc theory, the poloidal mass flux term is small enough to be ignored so that the radial mass flux is balanced by the m_rϕ and m_θϕ terms during the steady state.

Figure 11.

The angular momentum (the upper panel) and energy (the lower panel) budget for our fiducial run (V1000). Various components of the budget have been averaged over time (from t = 42T₀ to 52T₀) and integrated over space (θ from 0.59 to 2.55 to include both the accreting surface and the mid-plane region). The averaged quantities have also been multiplied by r^3.5 so that these quantities are almost flat in radii. The green dashed curve in the lower panel is −E_pot/2 for comparison. The black curve in each panel is the addition of all the four components.

In Fig. 11, these terms are integrated over θ from θ = 0.59 to 2.55 covering both the surface accreting region and the mid-plane region. Similar to the results in Zhu & Stone (2018), the wind stress (m_θϕ) plays a less important role in accretion than the r–ϕ stress. The m_θϕ term is ∼1/4 of the m_rϕ term around R ∼1. Thus, only 20 per cent of accretion is due to the θ–ϕ stress. On the other hand, this value is larger than 5 per cent in the simulation of Zhu & Stone (2018). Considering that this disc is thicker than the disc in Zhu & Stone (2018), it implies that wind plays a more important role for accretion in thicker discs. Nevertheless, most accretion is still due to the internal r–ϕ stress within the disc, as in Zhu & Stone (2018).

On the other hand, the disc wind seems to play a much more important role in the energy transport. Flux_r, Flux_θ, E_cool, and E_pot in the lower panel of Fig. 11 are the four terms on the right-hand side of equation (24). The traditional thin disc theory (equation 29) suggests that, far away from the inner boundary, the energy transport in the radial direction (the first two terms on the right-hand side of equation 29) actually adds the disc energy by an amount that is equal to half the released gravitational energy. The energy gain/loss in the poloidal direction is normally ignored. Thus, the total cooling rate is 1.5 times the released gravitational potential energy. However, our particular simulation suggests that energy transport in the radial direction (the red curve) is small compared with the energy loss in the poloidal direction by the wind (the blue curve). The wind carries half of the gravitational potential energy (the green curve) so that only the rest half gravitational potential energy needs to be radiated away (the cyan curve). Thus, the cooling rate is

$$\begin{eqnarray*} \langle Q_{\rm cool}\rangle =\frac{\dot{\widetilde{M}}v_{K}^2}{2r^2}, \end{eqnarray*}$$

(36)

which is roughly 1/3 of the value in the thin disc theory. This cooling rate is plotted as the green dashed curve in the lower panel of equation (24), and it agrees with simulations very well (even at the inner disc close to the inner boundary). Thus, the disc’s effective temperature in the simulation can be approximated by

$$\begin{eqnarray*} \sigma T_{eff}^4=\frac{GM\dot{M}}{8\pi R^3}. \end{eqnarray*}$$

(37)

Based on our simulations, such temperature estimate indeed agrees with the measured temperature at the τ_R ∼ 1 surface. The disc vertical structure at R = 0.1 au is shown in Fig. 12. At τ_R = 1 (the dotted line in the right panels), the value calculated using equation (37) (the blue curve in the temperature panel) agrees with the measured temperature very well.

Figure 12.

The disc vertical structure along R = 0.1 au with respect to τ starting from the disc surface (left-hand panels) or z starting from the mid-plane (right-hand panels). In the bottom panels, the crosses with the solid black curves are F_{r, θ}, while the crosses with the dashed black curves are F_r,z. All quantities are averaged over both the azimuthal direction and time (45–50 T₀ with a 0.1 T₀ interval). The red and blue curves in the temperature and F panels are the analytical solutions using equation (35) with two different fluxes. The red one uses the flux that is calculated with equation (18) and the measured |$\dot{M}$|⁠; the blue dashed curve uses the flux that is calculated with equation (37) and the measured |$\dot{M}$|⁠. The black dotted line labels where τ_R = 1 in simulations.

However, except for the similar T_eff, the temperature structure along z in simulations is very different from the temperature structure based on the analytical theory. First, the radiation flux in the θ-direction deviates significantly from the flux in the z direction when τ_R ≲1 (the bottom panels in Fig. 12). This is because the radiation from the inner disc (R < R₀) is so strong that the flux measured in the optically thin region at R₀ consists of a significant contribution from the disc inside R₀. Thus, we use the measured flux at τ_R ∼1 to represent the flux emitted by the local annulus at R. Secondly, the measured flux in either the z direction or the θ-direction rises much slower from the mid-plane to the τ_R = 1 surface than the models (red and blue solid curves) where the heating rate is proportional to the disc local density (equation 35). The measured radiative flux only rises quickly beyond one disc scale height. This difference is due to (1) energy transport by turbulence is as important as the radiative energy transport within the disc so that less temperature gradient is needed to radiate the thermal energy, as shown in the upper panel of Fig. 10; (2) both heating and accretion processes becomes more efficient at high above the disc mid-plane in our MHD simulations. Even with the similar emergent flux, the mid-plane temperature of the analytical α disc model is hotter than the measured mid-plane temperature in simulations by a factor of ≳3. This result is consistent with previous local radiation MHD simulations (Turner 2004; Hirose, Krolik & Stone 2006; Jiang, Stone & Davis 2014b), suggesting that, towards the disc surface, MHD heating becomes more efficient compared with heating in viscous models. Thirdly, the emergent flux at τ_R = 1 is significantly lower than the flux (red curves) estimated based on the traditional accretion disc theory (equation 17) using the measured disc accretion rate of |$4\times 10^{-4}\rm M_{\odot }\, yr^{-1}$|⁠. This is mostly due to the energy lost in the poloidal direction as discussed in Fig. 11. Equation (37) which has accounted for the energy loss in the poloidal direction agrees with the measured F_z at τ = 1 much better. We note that equation (37) only stands at the inner disc. As shown in the temperature panel of Fig. 6, the measured disc temperature is higher than the dotted line beyond R ∼0.2 au. This is probably due to the fact that the outer disc is irradiated by the inner disc so that it gets heated up.

4.3 Different field strengths and geometries

Since the disc temperature structure is self-consistently determined by the radiative transfer process in these simulations, the only major disc parameters that we can vary are the initial field geometry and strength. Thus, we carry out two additional simulations (V1e4 and T100) to explore how a weaker field or a toroidal field can affect the disc accretion.

The disc temperature, density, velocity, and magnetic field structures are shown in Fig. 13. Although these two simulations have similar temperature structures, one major difference which is quite noticeable in the middle panels is that disc wind fails to be launched in the net toroidal field simulations. In T100, disc material high above the atmosphere falls to the disc (green curves) instead of leaving the disc. Furthermore, the surface accreting region in T100 is much thinner if it exists at all. In the right-hand panels, V1e4 shows an extended surface accreting region with high B_ϕ and B_r values due to the surface accretion mechanism, while T100 only shows a thin region at the disc surface with noticeable B_ϕ and very weak fields above that. There is no large-scale organized fields in T100 either. The disc is dominated by turbulent fields in T100.

Figure 13.

Similar to Fig. 7 but for the V1e4 case at t = 55T₀ (upper panels) and the T100 case at t = 60T₀ (lower panels).

This lack of surface accretion in net toroidal field simulations is also evident in Fig. 14 where the radial profiles of various quantities are shown. In the T_total and α panels, the two simulations have similar values at the disc mid-plane for both T_Rϕ and α, while the vertically integrated T_Rϕ and α are significantly higher for V1e4. This indicates that V1e4 has a higher stress level at the disc atmosphere than that in T100. The magnetic field panel also shows that, while B² at the mid-plane is similar between two simulations, V1e4 has much stronger fields at the disc atmosphere. This leads to a higher accretion rate for V1e4 even though these two simulations have very similar turbulent levels at the disc mid-plane.

Figure 14.

Similar to Fig. 6 but for the V1e4 case at t = 55T₀ (black curves) and the T100 case at t = 60T₀ (red curves). In the temperature and 〈B²〉/2P_{mid, 0} panels, the solid curves are the mid-plane quantities and the dashed curves are the quantities along r at θ = 1.1 where the photosphere is. The black dotted line in the temperature panel is from equation (37) with an accretion rate of 10|$^{-4}\rm\, M_{\odot }\, yr^{-1}$|⁠.

The difference in disc wind is clearly shown in the vertical profiles of various quantities (Fig. 15). At the wind region above Z ∼ R, V1e4 has a much higher density than T100. The outflow nature of this region in V1e4 is clearly shown in the velocity panels, while this region in T100 is falling back to the disc. The magnetic fields and stresses are also very weak in the wind region of T100. Although there are some hints of surface accretion for T100 at z/0.1 au∼1 shown in the v_R panel, the density there is more than 5 orders of magnitude lower than the disc mid-plane (the ρ panel) so that the radial accretion of this surface is negligible in T100.

Figure 15.

Similar to Figs 8 and 9 but for the V1e4 case (black curves) and the T100 case (red curves).

Since net poloidal magnetic fields are essential for wind launching, it is important to understand how FU Ori’s inner disc acquires such strong poloidal fields (tens to hundreds of Gauss). Current disc theory suggests that net poloidal magnetic fields can be either from the central star’s magnetosphere, or inherited from the natal molecular cloud core. Königl, Romanova & Lovelace (2011) have carried out MHD simulations to study how the kG magnetosphere of FU Ori’s central star can interact with the fast accreting inner disc. They found that the magnetosphere truncation radius is pushed close to the central star, but the wind that is launched at the truncation radius is still largely consistent with the observed outflow properties (e.g. mass loss rate and speed). On the other hand, the detailed modelling for wind lines (Calvet et al. 1993; Milliner et al. 2019) suggests that the wind is launched from a much larger scale (disc wind). Thus, detailed synthetic observations for the simulations of Königl et al. (2011) are needed to test if these simulations are consistent with the observed line profiles. For the second scenario, inheriting magnetic fields from molecular cloud cores has been studied extensively for discs controlled by both ideal MHD and non-ideal MHD processes (Rothstein & Lovelace 2008; Guilet & Ogilvie 2012, 2013; Okuzumi, Takeuchi & Muto 2014; Bai & Stone 2017). Based on the simple field diffusion equation, the thin disc can lose the magnetic fields outwards quickly (Lubow, Papaloizou & Pringle 1994). But recent MHD simulations by Zhu & Stone (2018) found that the disc is quite thick for the perspective of the magnetic field structure, and the disc can actually transport field inwards slowly with time. Thus, FU Ori may gain poloidal magnetic fields from the outer disc during the low accretion state while material is piling up at the inner disc. When MRI is triggered at the disc mid-plane (Armitage et al. 2001; Zhu et al. 2009b), such strong fields lead to strong accretion.

5 DISCUSSION

After studying the disc’s physical structure, we will compare the simulations with existing observations regarding the disc temperature, magnetic fields, and disc wind.

5.1 Photosphere properties

Previous FU Ori SED modelling from Zhu et al. (2007) suggests that the disc’s effective temperature follows the standard viscous disc model and the disc’s maximum effective temperature is ∼6420 K. This temperature profile is plotted against the photosphere temperature (at τ_R = 1) in our simulations, shown in Fig. 16. Our fiducial model (V1000) has a similar maximum disc temperature as the observations, although its accretion rate (⁠|$\sim 5\times 10^{-4}\rm M_{\odot }\, yr^{-1}$|⁠) is twice the accretion rate used in Zhu et al. (2007) (2.4|$\times 10^{-4}\rm M_{\odot }\, yr^{-1}$|⁠). Considering that most disc luminosity comes from the hottest region, our fiducial model has a similar luminosity as the observation. All our simulations have flatter profiles compared with observations, which is due to the irradiation from the inner disc to the outer disc as discussed above. Thus, our simulations may need to be combined with a slightly different extinction curve from Zhu et al. (2007) to explain all the observations at different wavelengths. The photospheres in our simulated discs have densities of 10⁻¹⁰–10⁻⁹ g cm⁻³, and almost rotate at the local Keplerian speed.

Figure 16.

The temperature, density, and azimuthal velocity at the disc photosphere (τ_R = 1) along R for three simulations. The dashed curve in the left-hand panel is the effective temperature derived from FU Ori’s SED modelling (Zhu et al. 2007).

5.2 Comparison with magnetic Field Zeeman observations

Donati et al. (2005) use a high-resolution spectropolarimeter to measure circularly polarized light (Stokes V) from thousands of spectral lines for FU Ori. The circular polarized light is produced by Zeeman splitting that depends on both the field geometry and strength. The measured polarization signal corresponds to the line-of-sight magnetic field of ∼32 G. Together with some additional constraints on the disc parameters (e.g. 60° inclination) and theoretical disc wind models (Ferreira 1997), the detailed decomposition of the Stokes V into antisymmetric and symmetric components has put a much more stringent constraint on the magnetic fields of FU Ori. To summarize the findings (1) comparing the polarized light with the unpolarized light reveals that strong magnetic fields occupy |$\sim 20{{\ \rm per\ cent}}$| of the disc surface, and the magnetic plasma rotates ∼2–3 times slower than the local Keplerian velocity; (2) the vertical component of the magnetic fields (leaving the disc surface) is pointing towards us with a strength of ∼1 kG at 0.05 au; (3) the toroidal fields in the disc point to a direction which is opposite to the disc’s orbital rotation with a strength of ∼500 G at 0.05 au.

Although these measurements are consistent with previous resistive MHD simulations (Ferreira 1997) where the MRI turbulence is simplified by the resistivity parameters, we can now compare these observations directly with our first-principle radiation MHD simulations. We thus measure the magnetic field direction and strength at the τ_R = 1 surface in our simulations. The magnetic fields at R = 0.05 and 0.1 au are shown in Fig. 17. Please note the direction of the magnetic field in this figure. a_z is a parameter that equals 1 if B_z at the τ_R = 1 surface is pointing in a direction that is leaving the disc midplane and it is −1 if B_z is pointing towards the disc mid-plane. |$\hat{V}_{\phi }$| is the unit vector in the disc’s rotational direction. The reason that we express B_ϕ in this |$\vec{B}_{\phi }\cdot \hat{V}_{\phi }a_z$| form is due to the facts that we can view the disc from either the top or bottom side of the disc in Fig. 7 and the disc’s B_z can also be either aligned or anti-aligned with the angular momentum vector of the disc’s rotation. Let’s take the V1000 case as an example. As shown in the upper left-hand panel of Fig. 17, |$\vec{B}_{\phi }\cdot \hat{V}_{\phi }a_z$| (the solid black curve) is negative. If we observe the disc downwards from the upper side of the disc in Fig. 7, B_z is pointing to us so that a_z = 1. In this case, |$\vec{B}_{\phi }\cdot \hat{V}_{\phi }$| is negative implying that B_ϕ is in the opposite direction from the disc rotation. This can be seen in Fig. 7 where B_ϕ has negative values in the wind region. If we view the disc from the bottom and |$\vec{B}_{z}$| is pointing towards the disc mid-plane, a_z = −1 so that B_ϕ at the τ_R = 1 surface on this side of the disc is in the same direction as the disc rotation (as shown with the positive B_ϕ values at the bottom side of the wind region in Fig. 7). On the other hand, since we do not know if the rotational axis of the disc is aligned or anti-aligned with the magnetic fields (e.g. both Sun and Earth have magnetic reversals), we can reverse the field direction in simulations and the disc velocity structure will be unchanged. In that case, if we look at the disc downwards from the upper side of Fig. 7, a_z = −1 and B_ϕ at the wind region will be positive (in the same direction as the disc rotation) so that |$\vec{B}_{\phi }\cdot \hat{V}_{\phi }a_{z}$| is still negative.

Figure 17.

Upper panels: the vertical (blue curves) and azimuthal (black curves) components of magnetic fields measured at the τ_R = 1 surface at R = 0.05 au (solid curves) and 0.1 au (dashed curves) for three simulations (from left-hand to right-hand panels). a_z equals 1 if the B_z field at the τ_R = 1 surface is pointing in a direction that is leaving the disc mid-plane and equals −1 if the B_z field is pointing towards the mid-plane. |$\hat{V}_{\phi }$| is the unit vector in the disc’s rotational direction, and |$\vec{B}_{\phi }$| is the projection of the magnetic field vector to the disc’s rotational direction. Lower panels: the radial (blue curves) and azimuthal (black curves) velocity at the τ_R = 1 surface at R = 0.05 au (solid curves) and 0.1 au (dashed curves) for three simulations.

Our fiducial case (V1000) roughly reproduces the velocity and field geometries inferred from Donati et al. (2005). At R = 0.05 au, the τ_R = 1 surface is at z ∼ R which is the top of the surface accreting region or the bottom of the wind region (Fig. 7). At z ∼ R, the disc rotates with ∼60 per cent of the mid–plane Keplerian velocity (the lower left-hand panel of Fig. 17), while the disc becomes Keplerian slightly deeper in the disc (the V_ϕ panel in Fig. 8). Considering that the photospheres in other two cases are slightly deeper and they are Keplerian rotating, this ∼60 per cent of Keplerian rotation speed sensitively depends on the photosphere position and can be quite uncertain. At the τ_R = 1 surface of R = 0.05 au, the field strength is quite strong with B_z ∼ 150 G. If B_z is pointing to us, B_ϕ will be in a direction that is opposite to the disc rotation, which is consistent with observations. B_ϕ is half of B_z, which is also consistent with observations. At deeper regions in the disc, both B_ϕ and B_z decreases significantly. In the surface accreting region and down towards the disc mid-plane, B_ϕ changes from negative to zero and to positive. Thus, the 20 per cent covering factor from observations could be that 20 per cent light comes from the strong B and sub-Keplerian region, while the rest 80 per cent comes from the deeper Keplerian and weaker B region. The only difference between our simulations and the observations is that the field strength measured in simulations is weaker than the observed inferred kG strength by a factor of ∼5. On the other hand, we note that the first-order moment of the observed Zeeman signature is only ∼32 G. The kG strength is inferred from matching models considering the 60° inclination and the assumed filed geometry and filling factor. As will be shown in Section 5.3, the assumed inclination is too high compared with recent ALMA observations. Overall, the relatively good agreement regarding the field and velocity structure is very encouraging.

Our model also predicts that new observations by SpIROU at near-IR may reveal a different field structure than earlier results using optical lines from Donati et al. (2005) since near-IR lines come from further out in the disc (e.g. 0.1 au). The simulation indicates that the τ_R = 1 surface has very different field geometries and strengths at R = 0.1 au (the dashed curves in Fig. 17) compared with those at R = 0.05 au. From Fig. 7, we can see that, further away from the central star, the τ_R = 1 surface is closer to the disc mid-plane due to the lower disc surface density there. The upper left-hand panel in Fig. 8 shows that both B_z and B_ϕ at the τ_R = 1 surface change their signs moving from 0.05 to 0.1 au and the field strength gets a lot weaker. Furthermore, unlike at 0.05 au, B_ϕ is stronger than B_z at the photosphere of 0.1 au since the photosphere is at the bottom of the surface accreting region and closer to the disc mid-plane.

The surface accreting regions in our other two simulations, V1e4 and T100, have much lower density so that the τ_R = 1 surface is close to the disc mid-plane even at R = 0.05 au (Fig. 13). Thus, B_ϕ is always stronger than B_z at the photosphere as shown in the right two panels of Fig. 17. If B_z is pointing towards us, B_ϕ will be in the same direction as the disc rotation in these cases.

Various possible scenarios for B_z and B_ϕ measurements are summarized in Fig. 18. Under the surface accretion picture, B_z becomes quite strong at the upper surface/the base of the wind region at R ∼ z, and B_ϕ changes sign there. Thus, if the disc has a very high density so that the photosphere is only in the wind region or at the wind-base region (the thin dashed curve is the photosphere under this scenario), we are expecting to measure strong B_z and B_ϕ at all disc radii. On the other hand, the disc normally has a lower density at the outer cooler region and the opacity there is lower, it is more likely that the photosphere changes from the wind-base region at the inner disc to the lower surface/disc region at the outer disc (e.g. V1000 case). In this case, the B_z at the photosphere decreases dramatically at the outer disc and B_ϕ changes sign from the inner photosphere to the outer disc photosphere, indicating observations at different wavelengths may reveal different field and velocity geometries. For the third scenario that the photosphere is always closer to the disc (e.g. V1e4 and T100 cases), B_z will be significantly smaller than B_ϕ at all radii and observations at different wavelengths may reveal similar field and velocity geometries. We note that the signs of various B components can change depending on our viewing angle and the orientation between the fields and the rotational axis (as described in Fig. 18).

Figure 18.

The schematic plot showing B_z and B_ϕ measured at different radii under different scenarios (the thin red curve: the photosphere is always at the wind region; the middle thick red curve: the photosphere of the inner disc is at the wind region while the photosphere at the outer disc is within the disc; the lower thick red curve: the photosphere is always within the disc region). The signs of B follow the right-hand rule with respect to the angular momentum vector of the disc.

We want to caution that we use the τ_R = 1 surface to represent both the photosphere and where the magnetic fields are measured. In reality, the magnetic fields are measured by Donati et al. (2005) using a subset of G0 line list. These lines are likely to trace disc region that is above the photosphere. Detailed radiative transfer modeling with lines is needed to compare our simulations with observations.

5.3 Comparison with disc wind observations

FU Ori shows evidence of strong winds in P Cygni profiles, especially in the Na i resonance lines (Bastian & Mundt 1985; Croswell, Hartmann & Avrett 1987). The blue-shifted line absorption implies a disc outflow with a typical velocity of 100–300 km s⁻¹ and a mass-loss rate of |$\sim 10^{-5}\rm M_{\odot }\, yr^{-1}$| (Calvet et al. 1993). Recent work by Milliner et al. (2019) suggests that the wind may be turbulent.

We have plotted the gas radial velocity and mass-loss rate at different poloidal directions in Fig. 19. As long as the disc is threaded by net vertical fields, the magnetic fields accelerate the gas flow along the radial direction, reaching ∼400 km s⁻¹ terminal velocity. The integrated outflow rate at a distance r from the central star is

$$\begin{eqnarray*} \dot{M}_{\rm wind}(r)&=&\int _{0}^{2\pi }{\rm d}\phi \int _{0}^{\pi }{\rm d}\theta r^2\sin (\theta)\rho v_{r}\\ &=&\int _{0}^{\pi }2\pi r^2\sin (\theta)\langle \rho v_{r}\rangle {\rm d}\theta , \end{eqnarray*}$$

(38),(39)

where 〈〉 means that the quantities have been averaged over the azimuthal direction. The lower left-hand panel of Fig. 19 shows that 2|$\pi$|r²sin (θ)〈ρv_r〉 is around |$10^{-5}\rm M_{\odot }\, yr^{-1}$|⁠. Thus, the integrated wind loss rate from the pole to 30° (0.52 in Radian) away from the pole is |$\sim 10^{-5}\rm M_{\odot }\, yr^{-1}*0.52*2\sim 10^{-5}\rm M_{\odot }\, yr^{-1}$| where 2 comes from both sides of the disc. Thus, our fiducial simulation can reproduce both the observed outflow velocity and outflow rate.

Figure 19.

The radial velocity (upper panels) and mass-loss rate (lower panels) at 0.2 and 1 au along the θ direction in our three simulations (from left to right). The quantities have been averaged over both time (the last 2T₀ of each simulation) and azimuthal direction.

If the disc is threaded by net toroidal fields, wind cannot be launched, as shown in the right-hand panel of Fig. 19. Thus, the existence of disc wind in FU Ori implies that the disc is threaded by net vertical magnetic fields.

5.4 New FU Ori parameters

While we are preparing this manuscript, the distance to FU Ori is more precisely constrained by Gaia. The new distance is 416 ± 9 pc (Gaia Collaboration et al. 2018) instead of 500 pc assumed in Zhu et al. (2007). The disc inclination is also better constrained to be 35° by ALMA (Pérez et al. 2020) instead of 55° assumed in Zhu et al. (2007). With these updated parameters, Pérez et al. (2020) derive that the central star mass is updated to be 0.6 M_⊙ instead of 0.3 M_⊙, and the disc accretion rate is 3.8|$\times 10^{-5}\rm M_{\odot }\, yr^{-1}$| instead of 2.4|$\times 10^{-4}\rm\, M_{\odot }\, yr^{-1}$|⁠. The disc accretion rate now is only 1/6 of the earlier estimate due to the fact that both the closer distance and more face-on configuration reduce the disc accretion rate estimate. In the Appendix (Fig. A1), we have shown the SED fitting using the new parameters.

To be consistent with these new parameters, we have carried out a simulation which is similar to the V1e4 case but with M_* = 0.6M_⊙. The results are shown in Fig. 20. The overall ‘surface accretion’ picture still stands. But due to the short duration of this simulation (only to 31.5 T₀), the field structure at the surface accreting region is not fully established. The high disc accretion rate and the high central star mass release a significantly amount of gravitational energy so that the disc is significantly hotter than the V1e4 case with M_* = 0.3M_⊙. The real FU Ori system may have weaker net vertical fields or a lower surface density than those we assumed in Fig. 20.

Figure 20.

Similar to Figs 7 and 6 but for the new FU Ori parameters at t = 31.5T₀. In the lower left-hand panel, the dashed curve is the temperature at θ = 0.9 where the photosphere is.

6 CONCLUSIONS

We have carried out 3D global ideal MHD simulations to study the inner outbursting disc of FU Ori. Since the accretion disc outshines the central star, the radiation field of the disc plays an important role in the disc accretion dynamics. The radiative transfer is also crucial for connecting with observations. Thus, we self-consistently solve the radiative transfer equations along with the fluid MHD equations. We have carried out simulations where the disc is threaded by either net vertical or net toroidal magnetic fields.

We find that, when the disc is threaded by net vertical fields, most accretion occurs in the magnetically dominated atmosphere at z ∼ R, very similar to the ‘surface accretion’ mechanism in previous simulations with the simple locally isothermal equation of state. This implies that the ‘surface accretion’ is a general feature of accretion discs threaded by net vertical fields. The disc mid-plane shows spiral arms while the disc surface has filamentary structures. With radiative transfer included, we can study the accretion disc’s temperature structure. The radiation pressure is |$\sim 30{{\ \rm per\ cent}}$| of the gas pressure at the inner disc (e.g. 0.1 au). The disc mid-plane has a sharp temperature transition at ∼0.15 au separating the inner and outer discs that are at the higher and lower branches of the equilibrium ‘S’ curve. But the accretion and stress profiles are smooth despite the jump of disc temperature. This implies that the global accretion structure is mainly controlled by the global geometry of magnetic fields and is insensitive to the disc local temperature.

Compared with the simulations for thinner discs in Zhu & Stone (2018), the simulations here have stronger disc wind. 20 per cent of disc accretion is due to the wind θ–ϕ stress, which is higher than 5 per cent in Zhu & Stone (2018). The wind mass loss rate from the disc surface spanning one order of magnitude in radii is 1–10 per cent of the disc accretion rate, which is also higher than 0.4 per cent in Zhu & Stone (2018). Thus, the disc wind seems to be stronger in thicker discs. The mass loss rate of |$\sim 10^{-5}\rm M_{\odot }\, yr^{-1}$| in our FU Ori simulations is consistent with observations. The wind’s terminal speed is ∼300–500 km s⁻¹. This speed is also consistent with the observed wind speed and is several times the Keplerian speed at the launching point (V_K at the inner boundary is 100 km s⁻¹). On the other hand, no disc wind is launched when the disc is threaded by net toroidal fields, implying that net vertical fields are crucial for launching the disc wind. The net toroidal field simulation also shows weaker accretion and smaller vertically integrated stresses due to the lack of the surface accretion at the disc surface.

The moderate disc wind also carries half of the accretion gravitational potential energy so that only the rest half of gravitational potential energy needs to be radiated away. The emergent flux is only ∼1/3 of the traditional value with the same disc accretion rate (comparing equation 37 with equation 18). Thus, the disc photosphere temperature is lower than that predicted by the thin α-disc model having the same accretion rate. Thus, using the observed flux, the previously inferred disc accretion rate may be lower than the real disc accretion rate by a factor of ∼2–3. The disc mid-plane is also much cooler than that predicted by viscous models due to the energy transport by turbulence at the mid-plane and the efficient heating at the disc surface. With the surface accretion, the disc is heated up at the surface and the energy there can be more easily radiated away.

We have compared the magnetic fields at the photosphere in our simulations with Zeeman observations from Donati et al. (2005). The disc’s τ_R = 1 photosphere can be either in the wind launching region or the accreting surface region, depending on the accretion rates and the disc radii. Magnetic fields have drastically different directions and magnitudes between these two regions. It is very encouraging that the photosphere in our fiducial model, which is at the base of the wind launching region, agrees with previous Zeeman observations regarding both the magnetic field direction and magnitude. On the other hand, we suggest that the magnetic fields probed by future Zeeman observations at different wavelengths (e.g. near-IR) or for different systems (e.g. with lower accretion rates) can be quite different from the existing measurements in Donati et al. (2005) since the photosphere can be deep into the surface accreting region.

Overall, we find excellent agreements between the first-principle MHD simulations having net vertical fields and existing observations regarding both the wind and magnetic field properties. This strongly supports that accretion discs in FU Orionis systems are threaded by net vertical magnetic fields and MHD processes are important for the accretion process. More comparisons between simulations and future observations will allow us to probe the 3D structures of magnetic fields and gas flow in accretion systems.

ACKNOWLEDGEMENTS

All simulations are carried out using computer supported by the Texas Advanced Computing Center (TACC) at the University of Texas at Austin through XSEDE grant no. TG-AST130002 and from the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. ZZ acknowledges support from the National Science Foundation under CAREER grant no. AST-1753168. The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation.

Footnotes

1

In poorly ionized discs where the non-ideal MHD effects become important, hydrodynamical processes may also play an important role in disc accretion (Turner et al. 2014).

2

Current theory includes fragmented clumps (Vorobyov & Basu 2006), spiral arms from gravitational instability (Armitage, Livio & Pringle 2001; Zhu, Hartmann & Gammie 2009a; Martin et al. 2012; Bae et al. 2014; Kadam et al. 2019), or binary interaction (Bonnell & Bastien 1992).

APPENDIX: SED FITTING FOR FU ORI

With the updated FU Ori inclination, Pérez et al. (2020) use the disc atmospheric radiative transfer model (Zhu et al. 2007) to update FU Ori’s parameters. The best-fitting SED is shown in Fig. A1.

Figure A1.

With the new FU Ori distance from Gaia and disc inclination from ALMA, FU Ori’s disc parameters have changed moderately (Pérez et al. 2020). This shows the new SED fit using the updated FU Ori parameters (Pérez et al. 2020).