Time-dependent global simulations of a thin accretion disc: the effects of magnetically driven winds on thermal instability

ABSTRACT

1 INTRODUCTION

Observations on X-ray binaries have suggested that except for GRS 1915+105 and IGR J17091−3624, the high/soft state of X-ray binaries is quite stable in the luminosity (L) range of 0.01–0.5 Eddington luminosity (L_Edd; Belloni et al. 2000; Done, Wardziński & Gierliński 2004; Altamirano et al. 2011). This is inconsistent with the standard thin disc theory, which predicts that the radiation-pressure-dominated inner region of a thin disc is both thermally and viscously unstable when the accretion rate reaches 0.06L_Edd (Shakura & Sunyaev 1973; Piran 1978). For GRS 1915+105 and IGR J17091−3624, they are known to exhibit a limit-cycle light curve (Belloni et al. 2000; Janiuk, Czerny & Siemiginowska 2002; Altamirano et al. 2011), which is believed to arise from the thermal instability present in the thin disc (Belloni et al. 1997; Szuszkiewicz & Miller 2001; Li, Xue & Lu 2007). Although the time-scale of the limit-cycle light curve in active galactic nuclei (AGNs) is too long to be directly observed, the outbursts in some young radio galaxies appear to be triggered by the thermal instability (Czerny et al. 2009; Wu 2009). Additionally, the periodic outbursts observed in the repeating changing-look (CL) AGNs are also thought to be a consequence of the thermal instability in the intermediate zone between an inner advection-dominated accretion flow and outer thin disc (Sniegowska et al. 2020; Pan, Li & Cao 2021). These studies at least indicate that it is common for the thin disc of L ≳ 0.01L_Edd to be in a thermally stable state (Gierliński & Done 2004). Regarding the limit-cycle behaviour, the intermittent outbursts in some young radio galaxies, and the periodic outbursts in the repeating CL AGNs, two possibilities can be considered: (1) there might be other physical processes for the outbursts; (2) the occurrence of thermal instability in the thin disc is contingent upon specific conditions.

The thermal instability of a thin disc is still an open issue. In theory, eliminating the thermal instability of thin discs is very necessary. The thermal instability arises from an imbalance between heating and cooling in the radiation-pressure-dominated inner region (Pringle 1976; Shakura & Sunyaev 1976). Therefore, an in-depth investigation of the mechanisms responsible for heating and cooling is vital for understanding and resolving this issue. The heating processes in thin discs are believed as releasing the gravitational energy through the turbulence driven by the magnetorotational instability (MRI; Balbus & Hawley 1991). In the standard thin disc theory, the α-viscosity is employed to mimic the turbulence (Shakura & Sunyaev 1973). Some works have studied that when the viscous stress in the disc is solely proportional to the gas pressure, the disc is stable (Lightman & Eardley 1974); when the viscous stress is proportional to the total pressure of gas and radiation, the disc is unstable (Shakura & Sunyaev 1976). This point has been confirmed by Ohsuga (2006) through radiation hydrodynamical (RHD) simulations. Furthermore, radiation magnetohydrodynamical (RMHD) simulations conducted by Hirose, Blaes & Krolik (2009) supported that the viscous stress is proportional to the total pressure. Jiang, Stone & Davis (2013) have also found through the RMHD simulations that the disc is thermally unstable. In theory, for the disc to be stable, it is necessary to consider processes that cool the disc or decrease the relative importance of radiation pressure compared to the total pressure. For instance, Zheng et al. (2011) suggest that the magnetic pressure contributes to the vertical support force, reducing the relative importance of radiation pressure and thus stabilizing the disc. This finding has been confirmed by numerical simulations implemented by Sądowski (2016). Li & Begelman (2014) propose that the magnetically driven winds cool the disc, promoting stability. Jiang, Davis & Stone (2016) find that the presence of an iron opacity ‘bump’ can absorb heat in thin discs of AGNs, thereby maintaining disc stability.

Previous studies have suggested that magnetic pressure and magnetically driven flows play a important role in eliminating thermal instability (Merloni 2003; Zheng et al. 2011; Li & Begelman 2014; Habibi & Abbassi 2019). However, these studies primarily focused on analysing the instability and did not investigate the time-dependent evolution of the thin disc with poloidal and toroidal magnetic fields. The toroidal magnetic fields can support the thin disc by replacing radiation pressure (Zheng et al. 2011), while the poloidal magnetic fields can produce winds/outflows that cool the thin disc and transfer angular momentum (Li & Begelman 2014). In this work, we aim to simulate the time-dependent evolution of the thin disc with magnetic fields to study the effects of magnetically driven winds on thermal instability.

The paper is organized as follows. In Section 2, we describe basic equations. In Section 3, we present our results and related discussions. In Section 4, we give conclusions and discussion.

2 BASIC EQUATIONS

In the standard thin disc around a black hole (BH), thermal instability may occur. Compared to the previous works on the time-dependent global simulations of thermal instability, our global simulations consider the effect of magnetic fields. The magnetic fields include the toroidal (B_ϕ) and poloidal (B_p) components. At the disc surface, the poloidal component is B_p ≡ B_re_r + B_ze_z and |${B}_{\rm p}\equiv \sqrt{B_{\mathrm{ r}}^2+B_{z}^2}$|⁠, where B_r is the radial component at the surface and B_z is the vertical component, respectively. In order to study the effect of the magnetocentrifugal-force-driven winds, the inclination angle of the field lines at the disc surface is requested to be less than 60° with respect to the cold disc surface (Blandford & Payne 1982). Therefore, |$B_{\rm z}=\sqrt{3}B_{\rm r}$| at the disc surface is employed. In the accretion disc, the total pressure (p) at the mid-plane includes the gas pressure (p_g), the radiation pressure (p_r), and the magnetic pressure (⁠|$p_{\rm m}=\frac{B^2}{8 \pi }$|⁠), and reads

Here, |$p_{\rm g}=k\rho T/(\overline{\mu } m_{\rm H})$|⁠, where k, m_H, |$\overline{\mu }$|⁠, ρ, and T are the Boltzmann constant, the hydrogen atom mass, the mean molecular weight, the density at mid-plane, and the temperature at mid-plane, respectively. We set |$\overline{\mu }=0.62$|⁠. The radiation pressure is determined by

where c is the light speed, τ_R the Rosseland mean optical depths from the mid-plane to the disc surface, and F⁻ is the radiative flux away from the disc surface, respectively. The radiative flux reads

where σ is the Stefan–Boltzmann constant, and τ_p the Plank mean optical depths from the mid-plane to the disc surface, respectively. Equations (2) and (3) are valid for both optically thick and thin regimes (Szuszkiewicz & Miller 1998). When the accretion disc is extremely optically thick, i.e. τ_R and τ_p are much larger than unity, equations (2) and (3) can reduce to the standard expressions in the standard thin disc. For τ_R and τ_p, we follow Szuszkiewicz & Miller (1998) and have

and

Here, Σ is the surface density of the accretion disc and reads Σ = 2ρH, where H is the half-thickness of the accretion disc.

The global evolution of the thin disc can be described using a set of hydrodynamical equations in cylindrical polar coordinates (r, ϕ, z), where r is the disc radius, z is the symmetric axis of the disc, and ϕ is the azimuth, respectively. The mid-plane of the disc is the plane at z = 0. For simplicity, we assume the disc to be axisymmetic (i.e. |$\partial /\partial \phi \equiv 0$|⁠) and vertically average the hydrodynamical equations. Therefore, the basic quantities of the disc properties, such as surface density (Σ), radial velocity (v_r), specific angular momentum (l), and temperature (T), are functions of time and radius. According to conservation of mass, the evolution of the surface density reads

where |$2\dot{m}_{\rm w}$| represents the mass-loss rate on unit area of disc surface due to the magnetocentrifugal-force-driven winds (Pan et al. 2021) and is given by

where μ is the mass loading parameter and Ω is the angular velocity, respectively. The evolution of radial velocity reads

where l_k is the Kepler angular momentum. In this work, we employ pseudo-Newton potential to mimic the general relativistic effects (Paczyńsky & Wiita 1980) and then |$l_\mathrm{ k}={\sqrt{GM_{\mathrm{ BH}}r^3}}/{(r-r_{\rm s})}$|⁠, where G is the gravitational constant, M_BH is the BH mass, and |$r_{\rm s}\equiv \frac{2GM_{\mathrm{ BH}}}{c^2}$| is the Schwarzschild radius, respectively. The evolution of specific angular momentum reads

where c_s (⁠|$\equiv \sqrt{p/\rho }$|⁠), α, and T_m are the local sound speed, the α parameter in the viscosity stress tensor, and the magnetic torque carried by the magnetic winds, respectively. T_m is given by

The evolution of temperature reads

where we define β₁ = (p_g + p_r)/p_m and β₂ = p_g/(p_g + p_r) to measure the strength of magnetic fields and the ratio of gas pressure to the sum of gas and radiation pressure, respectively. Since the strength of magnetic fields cannot be determined inherently here, β₁ is a free parameter at t = 0. β₂ depends on density and temperature. Q⁺ is the viscous heating rate and is given as

Equations (6), (8), (9), and (11) are consistent with conservation of mass, radial momentum, angular momentum, and energy, respectively. To describe the evolution of H, we follow Li et al. (2007) and add two additional equations, as follows:

where v_z is the vertical velocity of the disc surface. When v_z > 0, the disc expands vertically. When v_z < 0, the disc shrinks vertically. The assumption of the vertical hydrostatic equilibrium is abandoned in this work. Therefore, equations (6), (8), (9), (11), (13), and (14) make up a set of equations to be solved. When magnetic fields are neglected, the equations are consistent with the equations in Li et al. (2007).

The evolution of magnetic fields in the disc is essentially three-dimensional. The turbulence driven by MRI plays an important role not only in angular momentum transfer but also in the diffusion of magnetic fields (Balbus & Hawley 1991; Fromang & Stone 2009; Guan & Gammie 2009; Cao & Spruit 2013). In our one-dimensional models, we cannot accurately describe MRI, while we can phenomenologically describe the advection and diffusion of magnetic fields. Phenomenologically, the change in the magnetic field strength is caused by the change in the surface density and the half-thickness. According to Li & Begelman (2014)’s assumption, we have

where the subscript t and t+△t indicate that the values are physical quantities at time t and t+△t, respectively, and ϵ is a parameter depending on the diffusion and advection time-scales of magnetic fields. When ϵ = 0, corresponding to the case where the diffusion time-scale is far smaller than the advection time-scale, the strength of magnetic fields is a constant with the change of the surface density; when ϵ = 1, corresponding to the case where the advection time-scale is far smaller than the diffusion time-scale, all the magnetic flux (Φ = B/Σ) is effectively advected inward, and the flux is a constant (Li & Begelman 2014). Therefore, 0 ≤ ϵ ≤ 1 is adopted. On the other hand, when the half-thickness (H) becomes larger (i.e. the disc expands vertically), the magnetic fields in the disc become weaker. This is supported by the MHD numerical simulations of a hot accretion flow (Machida, Inutsuka & Matsumoto 2006). As in Zheng et al. (2011), the toroidal component of magnetic fields is taken as a function of H, and given as

where γ is a free parameter. Equation (16) can keep B_ϕH^γ constant with time. It is assumed that γ ≥ 1. This implies that as the disc expands vertically, the toroidal magnetic field becomes weaker. When γ = 1, the toroidal magnetic flux per unit radius is conserved.

3 NUMERICAL METHODS AND INITIAL CONDITION

The dynamical quantities of the accretion disc, Σ, v_r, l, T, H, and v_{|$\mathrm{ z}$|}, are obtained by solving a set of dynamical equations, which consist of the equations (6), (8), (9), (11), (13), and (14). The evolution of magnetic fields can be obtained from equations (15) and (16). For the time integration, a third-order total variation diminishing Runge–Kutta scheme is adopted to integrate the first two time steps only in the initial time (t = 0) and then the third-order backward differentiation explicit scheme is used for the remaining computations to speed up the calculation. The standard Chebyshev pseudo-spectral method is implemented in discretizing the spatial difference. Computational domain is set to be between 2r_s and 10⁴r_s. We have modified the code in Li et al. (2007) to include the effect of magnetic fields and implement the present simulations. More details on numerical methods are referred to in Li et al. (2007).

Our model parameters include M_BH, the accretion rate at the outer boundary (⁠|$\dot{M}$|⁠), α, μ in equation (7), ϵ in equation (15), γ in equation (16), and β_1,0, where β_1,0 is the initial value of β₁. In the present simulations, we set M_BH = 10 M_⊙, |$\dot{M}=0.1\dot{M}_{\rm Edd}$| (⁠|$\dot{M}_{\rm Edd}\equiv 10 L_{\rm Edd}/c^2$| is the Eddington accretion rate), α = 0.1, ϵ = 0.4, and γ = 1, respectively. Then, μ and β_1,0 become two main free parameters in our models. A transonic steady thin disc around a BH with |$0.1\dot{M}_{\rm Edd}$| is taken as our initial condition. β_1,0 determines the strength of initial magnetic field. Moreover, we calculate two cases of B_ϕ,0 = 0.1B_p,0 and B_ϕ,0 = 10B_p,0, which represents the two cases where the polar and toroidal component dominate, respectively. For the initial poloidal component, we assume |$B_{\mathrm{ z},0}=\sqrt{3}B_{\mathrm{ r},0}$|⁠. In this way, we can determine the initial structure of magnetic fields. Table 1 gives the basic parameters of our models.

4 RESULTS

As the initial model for numerical simulations, we employ a transonic thin disc without magnetic fields (Matsumoto et al. 1984; Abramowicz et al. 1988). We firstly calculate the evolution of the thin disc without magnetic fields and then calculate the evolution of the thin disc with weak magnetic fields.

4.1 Evolution of the thin disc without magnetic fields

In model M1, magnetic fields are ignored. Fig. 1 shows the evolutions of the luminosity (L), the half-height (H), the surface density (Σ), and the temperature (T), respectively. This model exhibits the limit-cycle behaviour, which is consistent with Li et al. (2007). All physical quantities evolve periodically, with one cycle lasting for 680 s. When the disc evolves to an approximate thin disc in a low-luminosity state, we set t = 0. The material is slowly accreted by a BH from the outer zone to the inner zone with time, making the surface density and temperature in the inner zone increase. This increase in radiation pressure (∝ T⁴) triggers thermal instability in the inner zone. At around 452 s, the temperature rises rapidly, the disc height significantly expands, and the density peak drifts outward, carrying a large amount of material moving outward. The disc remains in a high-luminosity state for about 40 s.

Figure 1.

Evolutions of the disc luminosity [panel (A)], half-height [panel (B)], surface density [panel (C)], and temperature [panel (D)] in model M1 without magnetic fields.

Open in new tabDownload slide

4.2 Effects of magnetically driven winds

In models M2–M33, the magnetic fields are dominated by a poloidal magnetic field and we set B_ϕ,0 = 0.1B_p,0. In models M34–M45, the magnetic fields are dominated by a toroidal magnetic field and we set B_ϕ,0 = 10B_p,0. The large-scale poloidal magnetic field can produce winds, and the mass outflow rate is described by equation (7), which indicates that a stronger poloidal field can produce more powerful magnetic winds, and larger μ enhances the mass outflow rate of winds. To vary the initial strength of the magnetic fields in our models, we adjust the value of β_1,0.

In our simulations, we initially superimpose a weak magnetic field on the initial model. Superimposing the weak magnetic field is to avoid disrupting the structure of the initial model caused by the superimposed magnetic field. Moreover, when we superimpose a stronger magnetic field, our calculation becomes unstable. Therefore, β_1,0 is taken from 100 to 1000.

Fig. 2 shows the evolution of luminosity in some typical models. Panel (A) displays three models that exhibit limit-cycle behaviour, while panel (C) shows four models that do not display any outburst. In panel (B), two models experience an outburst and then keep calming down. This indicates that these models are in a transitional state from being completely unstable to completely stable. We think that these models are ultimately stable. Generally, when considering the same β_1,0, a larger value of μ is needed to make the disc stable. For example, the thin disc remains stable and does not transition to a high luminous state when μ exceeds 1.5 for β_1,0 = 100, and μ exceeds 10 for β_1,0 = 1000. If μ is small enough to neglect the effect of magnetic winds, the disc is still thermally unstable and exhibits the limit-cycle behaviour. In panel (A), the limit-cycle behaviour is displayed in models M1, M2, and M30. For the model M2 with β_1,0 = 100 and μ = 0.1, the model M30 with β_1,0 = 1000 and μ = 1, the limit-cycle behaviour takes place. Compared to model M1 without magnetic fields, the limit-cycle periods in these two models are shorter. Due to the weaker magnetic fields in model M30 compared to model M2, the limit-cycle period in model M30 is slightly longer than that in M2.

Figure 2.

Evolutions of luminosity in the models M1–M5 and M30–M33 with the magnetic field structure of B_ϕ,0 = 0.1B_p,0. Panel (A) shows that the three models exhibit the limit-cycle behaviour. Panel (B) shows that the two models experience an outburst and then calm down the next outburst. Panel (C) shows that the four models can calm down the outburst triggered by the thermal instability.

Open in new tabDownload slide

Fig. 3 further shows the dependence of the limit-cycle period on the parameter β_1,0: the models with the stronger magnetic fields have slightly shorter periods. In particular, the limit-cycle period in model M2 has decreased by about 20 per cent compared to model M1, as shown in Table 1. The limit-cycle period in the B_ϕ = 10B_p case is shorter than that in the B_ϕ = 0.1B_p case due to the stronger magnetic torque produced by the magnetic field structure in B_ϕ = 10B_p case (equation 10). The variation of period agrees with that in Pan et al. (2021) qualitatively, who numerically studied a local radiation-pressure-dominated thin disc with a poloidal magnetic field for understanding repeating CL AGNs. The limit-cycle period is reduced by the poloidal magnetic field because the magnetically driven winds take away the angular momentum of the disc and then reduce the time interval for the material to fall back to the outburst region.

Figure 3.

The dependence of the limit-cycle period on the parameter β_1,0. The red line corresponds to the models of B_ϕ = 0.1B_p, while the blue line corresponds to the models of B_ϕ = 10B_p.

Open in new tabDownload slide

Fig. 4 shows the dependence of the luminosity variation amplitude on the parameter β_1,0. In our models, the highest luminosity is almost the same for different μ and β_1,0, and the luminosity variation also is not significant. In the B_ϕ = 10B_p case, the luminosity variation amplitude is slightly smaller than that in B_ϕ = 0.1B_p case. However, in Pan et al. (2021), the highest luminosity and variation amplitude are different for different μ and β_1,0. This is because our calculations are based on global simulations instead of the local area in thin disc, called the transition zone in Pan et al. (2021). In our simulations, when an outburst takes place in the disc, the height and temperature of the inner region increase significantly, while the surface density decreases by outflow. This leads to an extremely weak magnetic field in the inner region (see equations 15 and 16). Consequently, the temperature and area of the hottest inner region are not significantly affected by the magnetic field during the outburst. Therefore, the maximum and minimum luminosity is almost the same in our models with different parameters.

Figure 4.

The dependence of the luminosity variation amplitude on the parameter β_1,0. The red line corresponds to the models of B_ϕ = 0.1B_p, while the blue line corresponds to the models of B_ϕ = 10B_p.

Open in new tabDownload slide

For the two models (M3 and M31) shown in panel (B) of Fig. 2, the initial introduction of magnetic fields does not immediately have a significant effect on cooling. Therefore, the disc firstly experiences an outburst, and at the end of the outburst, the disc’s height rapidly shrinks and strengthens the magnetic fields to produce strong winds, so that the disc is effectively cooled by the winds over time and finally becomes thermally stable like the models [shown in panel (C) of Fig. 2] with larger μ.

For the models in a thermally stable state [panel (C) of Fig. 2], the disc luminosity gradually increases, reaches a maximum value, and then gradually decreases. As an example, we show the evolution of the disc structures of model M4 in Fig. 5. For the other thermally stable models, their disc structures are similar in evolution. As shown in Fig. 5, during the first 1000 s of disc accretion, a large amount of material is carried to move into inside ∼200r_g [as shown in panel (A)], the disc height and temperature both increase, shown in panels (B) and (C).

Figure 5.

Evolutions of the disc structures in model M4 with β_1,0 = 100 and μ = 1.5. Panels (A)–(D) show the ratio of the surface density, half-height, temperature, and magnetic pressure to their initial values at different time, respectively. Panels (E) and (F) show the ratio of the poloidal and toroidal fields to their initial values at different time, respectively.

Open in new tabDownload slide

However, in model M4, the increase of the disc temperature is significantly suppressed compared to model M1. With the accumulation of material inside ∼200r_g, the magnetic fields in the inner region become stronger [as shown in panel (D)]. Panel (E) indicates that the poloidal field becomes stronger than the initial fields. Panel (F) indicates that the toroidal field within ∼15r_s becomes stronger, while the toroidal field beyond ∼15r_s becomes weaker. According to equations (15) and (16), the change of the toroidal field is determined by the surface density and disc height. Within ∼15r_s, the toroidal field enhancement is dominated by the surface density increase. Beyond ∼15r_s, the toroidal field weakening is dominated by the disc height increase. For the poloidal field, its strength change is only determined by the surface density. Therefore, the magnetic field enhancement is caused by the poloidal field enhancement. This results in the magnetically driven winds becoming so powerful that a significant portion of thermal energy is dissipated. Consequently, the disc temperature is constrained from further increasing, so the thermal instability is suppressed in model M4. After about 1000 s, more and more material is lost through the magnetically driven winds, resulting a reduction in the surface density and subsequent weakening of the magnetic fields. Simultaneously, the disc temperature and luminosity decrease gradually. This is attributed to the diminishing viscous heating rate, which is proportional to the surface density (see equation 12).

Fig. 6 shows the boundary between stable and unstable areas in the μ–β_1,0 parameter space for the B_{ϕ, 0} = 0.1B_p,0 case. The area above the line is stable for models, while the areas below the line are unstable. For the same strength of magnetic fields, increasing the mass loading parameter (μ) is helpful to suppress thermal instability. If the mass loading parameter is small enough, the disc becomes thermally unstable and exhibits the limit-cycle behaviour. In this case, the role of a poloidal magnetic field is to transfer the angular momentum and then change the period of the limit-cycle behaviour. In the B_ϕ,0 = 10B_p,0 case, the μ value on the boundary is 100 times higher than that in the B_ϕ,0 = 0.1B_p,0 case. As an example, model M36 becomes thermally stable when the μ value is 100 times higher than that in model M4.

Figure 6.

Boundary between stable and unstable areas in the μ–β_1,0 parameter space for the B_ϕ,0 = 0.1B_p,0 case. In the B_ϕ,0 = 10B_p,0 case, the μ value on the boundary is 100 times higher than that in the B_ϕ,0 = 0.1B_p,0 case.

Open in new tabDownload slide

5 CONCLUSION AND DISCUSSION

In this work, we calculate the global time-dependent evolution of a thin accretion disc with magnetic fields to investigate the effects of magnetically driven winds. The magnetically driven winds have two important effects on the thin disc. First, they facilitate the angular momentum transfer, which is related to the strength of magnetic fields. Secondly, these winds contribute to the cooling of the disc, which is related to the mass outflow rate of winds. In our models, the mass outflow rate is proportional to a mass loading parameter (μ). When μ is small, the mass outflow rate of winds is low, and although the winds play an important role in transferring angular momentum, they cannot be helpful to eliminate the thermal instability. In this case, due to the transfer of angular momentum, the limit-cycle period is shortened. When the initial magnetic field strength is one per cent of the sum of gas and radiation pressure, the limit-cycle period is shortened by about 20 per cent compared to the case without magnetic fields. When μ is large, the mass outflow rate of winds is high and then the winds carry away enough energy from the disc to suppress the outburst driven by the thermal instability. Therefore, the thin disc becomes thermally stable.

The mass loading parameter μ is an important parameter to determine whether the thin disc is thermally stable or not, while it is taken as a free parameter in our models. μ is proportional to ρv_p/B_p along magnetic field lines from the definition in Spruit (1996), and therefore it depends on local gas and field properties. In theory, the value of μ can be obtained based on the cold approximation of the Weber–Davis model (Weber & Davis 1967; Spruit 1996), where the thermal pressure is neglected. When models deviate from the cold approximation, what is the exact value of μ is an open issue. Unfortunately, numerical simulations have not given the range of μ. Under the cold approximation of the Weber–Davis model, the magnetic torque exerted by magnetically driven winds on the thin disc can be written as |$T_{\rm m}=3r B_{\rm p}^2\mu (1+\mu ^{-2/3})/4\pi$| (e.g. Cao & Spruit 2013). Combining with equation (10), we have |${B_{\phi }}/{B_{\mathrm{ p}}}=\frac{3}{2}\mu (1+\mu ^{-2/3})$|⁠. This is used to determine μ based on B_ϕ and B_p in Pan et al. (2021), and then μ = 2.9 × 10⁻⁴ for B_ϕ = 0.1B_p, while μ = 5 for B_ϕ = 10B_p. However, these values of μ are much smaller than the critical value required for the thermal stability of the disc, as determined in our calculations. There are probably factors to affect the loading of material on the magnetic field lines. Above the disc surface, it is often suggested that there is a corona. The corona is helpful to load material along the magnetic lines (Cao & Spruit 1994). Moreover, when the disc expands, the vertical velocity at the disc surface is helpful to push the material along the magnetic lines. Furthermore, as shown in the panel (A) and panel (B) of Fig. 5, the disc surface density increases faster compared to the disc height, suggesting an increase in density within the disc, so more material can be pushed along the magnetic lines. At the state of disc expanding, the μ value is expected to be larger than that given in the cold Weber–Davis model. However, the exact value of μ is unknown now, which is really worth using MHD simulations to further investigate this issue.

Fig. 5 [panel (B)] shows that a thermally stable disc experiences the processes of expanding and shrinking. In our simulations, the disc does not restore to its initial state after the disc shrinks, although it is expected that the disc should restore to its initial state. This is because strong winds always exist, whether the disc expands or shrinks. In fact, the mass loading parameter μ should be a function of radius and time. It increases when the disc expands, while it decreases when the disc shrinks. In this work, we set it constant with radius and time. This means that when the disc shrinks, the strong winds still exist, and then the disc in our simulations cannot be restored to its initial state.

ACKNOWLEDGEMENTS

XH was supported in part by the Natural Science Foundation of China (grant 11973018) and Natural Science Foundation of Chongqing (grant CSTB2023NSCQ-MSX0093). XL was supported by the Natural Science Foundation of Fujian Province (grant 2023J01008). SL was supported in part by the Natural Science Foundation of China (grant 12273089).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

References