Dynamics of magnetic flux tubes in accretion discs of T Tauri stars

TL;DR

This study models the behavior of magnetic flux tubes in T Tauri star accretion discs, revealing oscillation regimes, rise velocities, and potential impacts on disc mass loss and observed variability.

Contribution

It introduces detailed simulations of magnetic flux tube dynamics including buoyancy, drag, and heat exchange, highlighting new oscillation behaviors and their observational implications.

Findings

MFT oscillate and rise with velocities up to 15 km/s.

Approximately 20% of disc mass and flux can escape via buoyancy.

MFT oscillations may explain T Tauri stars' IR variability.

Abstract

Dynamics of slender magnetic flux tubes (MFT) in the accretion discs of T Tauri stars is investigated. We perform simulations taking into account buoyant, aerodynamic and turbulent drag forces, radiative heat exchange between MFT and ambient gas, magnetic field of the disc. The equations of MFT dynamics are solved using Runge-Kutta method of the fourth order. The simulations show that there are two regimes of MFT motion in absence of external magnetic field. In the region $r<0.2$ au, the MFT of radii $0.05 \leq a_0 \leq 0.16\,H$ ($H$ is the scale height of the disc) with initial plasma beta of 1 experience thermal oscillations above the disc. The oscillations decay over some time, and MFT continue upward motion afterwards. Thinner or thicker MFT do not oscillate. MFT velocity increases with initial radius and magnetic field strength. MFT rise periodically with velocities up to 5-15 km s$^{-1}$ and periods of $0.5-10$ yr determined by the toroidal magnetic field generation time. Approximately 20 % of disc mass and magnetic flux can escape to disc atmosphere via the magnetic buoyancy over characteristic time of disc evolution. MFT dispersal forms expanding magnetized corona of the disc. External magnetic field causes MFT oscillations near the disc surface. These magnetic oscillations have periods from several days to 1-3 months at $r < 0.6$ au. The magnetic oscillations decay over few periods. We simulate MFT dynamics in accretion discs in the Chameleon I cluster. The simulations demonstrate that MFT oscillations can produce observed IR-variability of T Tauri stars.

Tables

Table 1. Table 1: The characteristics of the accretion disc

$r$ [au]	${\rho }_m$ [g cm-3]	$T_m$ [K]	$T_a$ [K]	$H$ [au]	$B_e$ [G]
(1)	(2)	(3)	(4)	(5)	(6)
0.027	$2.0\times {10}^{-6}$	$4830$	$1700$	$6.2\times {10}^{-4}$	$170$
0.15	$4.1\times {10}^{-8}$	$2015$	$715$	$5.3\times {10}^{-3}$	$29.5$
0.2	$2.1\times {10}^{-8}$	$1840$	$625$	$7.7\times {10}^{-3}$	$22.4$
0.4	$3.7\times {10}^{-9}$	$1430$	$445$	$2.0\times {10}^{-2}$	$9.9$
0.6	$1.4\times {10}^{-9}$	$1240$	$360$	$3.3\times {10}^{-2}$	$6.3$

Table 2. Table 2: Parameters of runs

run	$r$ [au]	${\beta }_0$	$a_0$ [$H$]	$M$ [$M_{☉}$]	$\Phi$ [Mx]	$v_{\mathrm{surf}}$ [$\mathrm{km}{s}^{-1}$]	$a_{\mathrm{surf}}$ [$H$]	$E_m$, [erg]	$t_{\mathrm{surf}}$ [yr]	$t_{\mathrm{gen}}$ [yr]	${\dot{M}}_b,\left[M_{☉}{\mathrm{yr}}^{-1}\right]$	$T_{\min }/T_e$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
1	0.6	0.01	0.01	$2.8\times {10}^{-11}$	$2.9\times {10}^{21}$	0.6	0.034	$2.2\times {10}^{34}$	2.0	10.6	$2.6\times {10}^{-12}$	1.0
2	0.6	0.1	0.01	$2.5\times {10}^{-10}$	$2.7\times {10}^{21}$	0.5	0.035	$1.7\times {10}^{34}$	2.2	10.2	$2.5\times {10}^{-11}$	1.0
3	0.6	1	0.01	$1.4\times {10}^{-9}$	$2.0\times {10}^{21}$	0.2	0.046	$5.5\times {10}^{33}$	3.5	7.6	$1.8\times {10}^{-10}$	1.0
4	0.6	10	0.01	$2.5\times {10}^{-9}$	$8.6\times {10}^{20}$	0.05	0.059	$4.8\times {10}^{34}$	12.0	3.2	$1.8\times {10}^{-10}$	1.0
5	0.6	0.01	0.1	$2.8\times {10}^{-9}$	$2.9\times {10}^{23}$	5.2	0.336	$2.1\times {10}^{36}$	0.2	10.6	$2.6\times {10}^{-10}$	1.0
6	0.6	0.1	0.1	$2.5\times {10}^{-8}$	$2.7\times {10}^{23}$	3.6	0.355	$1.7\times {10}^{36}$	0.2	10.2	$2.5\times {10}^{-9}$	0.95
7	0.6	1	0.1	$1.4\times {10}^{-7}$	$2.0\times {10}^{23}$	1.5	0.463	$5.6\times {10}^{35}$	0.4	7.6	$1.8\times {10}^{-8}$	0.78
8	0.6	10	0.1	$2.4\times {10}^{-7}$	$8.6\times {10}^{22}$	0.3	0.578	$6.6\times {10}^{34}$	1.4	3.2	$7.9\times {10}^{-8}$	0.91
9	0.6	0.01	0.2	$1.1\times {10}^{-8}$	$1.1\times {10}^{24}$	9.0	0.673	$8.4\times {10}^{36}$	0.1	10.6	$1.0\times {10}^{-9}$	1.0
10	0.6	0.1	0.2	$1.0\times {10}^{-7}$	$1.1\times {10}^{24}$	5.5	0.709	$7.0\times {10}^{36}$	0.1	10.2	$9.8\times {10}^{-9}$	0.83
11	0.6	1	0.2	$5.6\times {10}^{-7}$	$8.1\times {10}^{23}$	2.2	0.925	$2.2\times {10}^{36}$	0.3	7.6	$7.4\times {10}^{-8}$	0.58
12	0.027	0.01	0.01	$6.3\times {10}^{-13}$	$7.0\times {10}^{19}$	0.9	0.028	$2.8\times {10}^{33}$	0.02	0.69	$9.1\times {10}^{-13}$	1.0
13	0.027	0.1	0.01	$5.8\times {10}^{-12}$	$6.8\times {10}^{19}$	0.8	0.029	$2.4\times {10}^{34}$	0.02	0.66	$8.8\times {10}^{-12}$	1.0
14	0.027	1	0.01	$3.2\times {10}^{-11}$	$5.0\times {10}^{19}$	0.4	0.036	$8.7\times {10}^{32}$	0.03	0.49	$6.5\times {10}^{-11}$	1.0
15	0.027	10	0.01	$3.2\times {10}^{-11}$	$2.3\times {10}^{19}$	0.07	0.047	$1.0\times {10}^{32}$	0.10	0.20	$2.9\times {10}^{-10}$	0.98
16	0.027	0.01	0.1	$6.3\times {10}^{-11}$	$7.0\times {10}^{21}$	8.6	0.276	$2.8\times {10}^{35}$	0.002	0.69	$9.1\times {10}^{-11}$	0.98
17	0.027	0.1	0.1	$5.8\times {10}^{-10}$	$6.8\times {10}^{21}$	6.5	0.290	$2.3\times {10}^{35}$	0.002	0.66	$8.8\times {10}^{-10}$	0.34
18	0.027	1	0.1	$3.2\times {10}^{-9}$	$5.0\times {10}^{21}$	2.8	0.356	$8.4\times {10}^{34}$	0.004	0.49	$6.5\times {10}^{-9}$	0.48
19	0.027	10	0.1	$5.8\times {10}^{-9}$	$2.3\times {10}^{21}$	0.5	0.431	$1.0\times {10}^{34}$	0.13	0.2	$2.9\times {10}^{-8}$	0.89
20	0.027	0.01	0.2	$2.5\times {10}^{-10}$	$2.8\times {10}^{22}$	15.6	0.554	$1.1\times {10}^{36}$	0.001	0.69	$3.6\times {10}^{-10}$	0.58
21	0.027	0.1	0.2	$2.3\times {10}^{-9}$	$2.7\times {10}^{22}$	10.3	0.576	$9.4\times {10}^{35}$	0.0011	0.66	$3.5\times {10}^{-9}$	0.35
22	0.027	1	0.2	$1.3\times {10}^{-8}$	$2.0\times {10}^{22}$	4.3	0.708	$3.4\times {10}^{35}$	0.0024	0.49	$2.7\times {10}^{-8}$	0.39

Table 3. Table 3: Parameters of YSO in Chameleon.

No	$M_{\star }$ [$M_{☉}$]	$\dot{M}$ [$M_{☉}{\mathrm{yr}}^{-1}$]	$L_{\star }$ [$L_{☉}$]	$T_{\mathrm{eff}}$ [K]	$R_{\star }$ [$R_{☉}$]	$P$ [d]	$r$ [au]	${\rho }_m$ [g cm-3]	$T_m$ [K]	$T_{\mathrm{irr}}$ [K]	$H$ [au]	$B_e$ [G]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
439	0.6	$4.8\times {10}^{-8}$	0.8	3669	2.2	32	0.2 (i)	$1.3\times {10}^{-8}$	1590	590	0.007	18.0
							0.42 (o)	$2.0\times {10}^{-9}$	1200	410	0.019	5.0
530	0.63	$2\times {10}^{-9}$	0.66	3955	1.7	35	0.07 (i)	$1.6\times {10}^{-8}$	1460	950	0.0013	20.4
							0.1 (o)	$5.6\times {10}^{-9}$	1220	800	0.0024	4.0

Equations

$$P+\frac{B^2}{8\pi }=P_{\mathrm{e}}.$$

$$P=\frac{R_{\mathrm{g}}}{\mu }\rho T$$

$$\Delta \rho ={\rho }_{\mathrm{e}}-\rho =\frac{B^2}{8\pi v_{\mathrm{s}}^2},$$

$$v_{\mathrm{s}}=\sqrt{\frac{R_{\mathrm{g}}{T}_{\mathrm{e}}}{\mu }}$$

$$f_{\mathrm{b}}=-\Delta \rho g_z,$$

$$\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}$$

$$M_{\mathrm{l}}$$

$$\Phi$$

$$\mathrm{d}Q$$

$$P+\frac{B^2}{8\pi }$$

$$\frac{\mathrm{d}{P}_{\mathrm{e}}}{\mathrm{d}z}$$

$$U$$

$$\frac{\mathrm{d}v}{\mathrm{d}t}$$

$$\frac{\mathrm{d}z}{\mathrm{d}t}$$

$$g_z=-z\frac{G{M}_{\star }}{r^3}{\left(1+\frac{z^2}{r^2}\right)}^{-3/2}.$$

$$f_{\mathrm{d}}=-\frac{{\rho }_{\mathrm{e}}{v}^2}{2}\frac{C_{\mathrm{d}}}{\rho \pi a},$$

$$f_{\mathrm{d}}=-\frac{\pi {\rho }_{\mathrm{e}}{\left({\nu }_{\mathrm{t}}a{v}^3\right)}^{1/2}}{\rho \pi a^2},$$

$${\nu }_{\mathrm{t}}=\alpha v_{\mathrm{s}}H,$$

$$a$$

$$B$$

$$\frac{\mathrm{d}U}{\mathrm{d}t}-\frac{P_{\mathrm{e}}}{{\rho }^2}\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(P+\frac{B^2}{8\pi }\right)$$

$$h_{\mathrm{c}}=\frac{\mathrm{d}Q}{\mathrm{d}t}$$

$$U_T\frac{\mathrm{d}T}{\mathrm{d}t}+\left(U_{\rho }-\frac{P_{\mathrm{e}}}{{\rho }^2}\right)\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\left(P_{\rho }+C_{\mathrm{m}}\rho \right)\frac{\mathrm{d}\rho }{\mathrm{d}t}+P_T\frac{\mathrm{d}T}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}T}{\mathrm{d}t}$$

$${\rho }_{\mathrm{e}}(z)$$

$$T_{\mathrm{e}}(z)$$

$${\Omega }_{\mathrm{k}}=\sqrt{\frac{G{M}_{\star }}{r^3}}$$

$$T_{\mathrm{a}}=T_{\mathrm{eff}}=280{\left(\frac{L}{L_{\text{sun}}}\right)}^{1/4}{\left(\frac{r}{1\mathrm{au}}\right)}^{-1/2}\mathrm{K},$$

$$z_{\mathrm{s}}=H\sqrt{\frac{2k}{k-1}\left(1-\frac{T_{\mathrm{a}}}{T_{\mathrm{m}}}\right)}.$$

$$h_{\mathrm{c}}=\frac{1}{\rho }\mathrm{div}\text{ q},$$

$$\text{ q}=-\kappa \nabla T.$$

$$\kappa =\frac{4{\sigma }_{\mathrm{R}}{T}^3}{3{\kappa }_{\mathrm{R}}\rho },$$

$$\mathrm{div}\text{q}=\frac{1}{r^{\prime }}\frac{\partial }{\partial r^{\prime }}\left(r^{\prime }q\right)\approx \frac{q_{\mathrm{ext}}-q_{\mathrm{in}}}{a},$$

Full text

††pagerange: Dynamics of magnetic flux tubes in accretion discs of T Tauri stars–References††pubyear: 2018

1 Introduction

A number of observations show that young stellar objects (YSO) have large-scale magnetic field. Investigations of Zeeman splitting and broadening of spectral lines of classical T Tauri stars have shown that the stars have magnetic field with strength of $1-3$ kG at their surface (Guenther et al., 1999; Johns-Krull, 2007). Donati et al. (2005) reported the registration of a magnetic field with strength of $\sim 10^{3}$ G near the inner edge of the accretion disc in FU Orionis using Zeeman splitting of spectral lines.

It is now possible to make polarization maps of the accretion discs at millimeter and (sub)millimeter wavelengths with ALMA. Interpretation of the polarization maps in frame of Davis-Greenstein mechanism allows to determine geometry of the magnetic field in the disc (Stephens et al., 2014; Li et al., 2016, 2018). Although there are other possible interpretations of the polarization maps (see Lazarian & Hoang, 2007; Tazaki et al., 2017; Kataoka et al., 2017; Stephens et al., 2017).

According to modern theory of star formation, young stars with the accretion discs form as a result of gravitational collapse of rotating molecular cloud cores with magnetic field (see reviews by Inutsuka, 2012; Li et al., 2014). Numerical simulations show that the magnetic flux of the molecular cloud cores is partially conserved during the collapse, so it is natural to assume that the magnetic field of the accretion discs of young stars is the fossil one (see reviews of the theory of the fossil magnetic field by Dudorov, 1995; Dudorov & Khaibrakhmanov, 2015). The magnetic field of accretion discs can also be a result of dynamo (see, for example, Brandenburg et al., 1995; Gressel & Pessah, 2015; Moss et al., 2016, and references therein).

Evolution of magnetic field in the accretion discs was usually investigated in kinematic approximation for prescribed uniform diffusivity (Bisnovatyi-Kogan & Ruzmaikin, 1976; Lubow et al., 1994; Agapitou & Papaloizou, 1996; Guilet & Ogilvie, 2014; Okuzumi et al., 2014). Dudorov & Khaibrakhmanov (2014) and Khaibrakhmanov et al. (2017) developed a magnetohydrodynamic (MHD) model of the accretion discs taking into account Ohmic, ambipolar diffusion and the Hall effect. They have shown that the magnetic field geometry varies through the disc. Inside the region of low ionization fraction (‘dead’ zone (Gammie, 1996)), Ohmic diffusion hinders amplification of the magnetic field, so that the magnetic field has poloidal geometry. Ambipolar diffusion operates in outer regions of the accretion discs, where the magnetic field asquires quasi-radial or quasi-azimuthal geometry depending on the intensity of ionization and grain parameters. The Hall effect operates near the borders of the ‘dead’ zones and leads to redistribution of the poloidal and toroidal components of the magnetic field. The magnetic field is frozen in gas near the inner edge of the accretion discs, where thermal ionization operates. Effective generation of the toroidal magnetic field is possible in this region. Dudorov & Khaibrakhmanov (2014) have suggested that magnetic buoyancy instability can solve the problem of runaway growth of toroidal magnetic field.

The magnetic buoyancy instability leads to formation of separate magnetic flux tubes (MFT) from the regular magnetic field (see Parker, 1979). The MFT can float from the region of their formation to the surface under the action of buoyant force. The instability has been found in numerical simulations of unstable gas layers with planar magnetic field (Cattaneo & Hughes, 1988; Matthews et al., 1995; Wissink et al., 2000; Fan, 2001). Vasil & Brummell (2008) have shown that the instability also arises in the gas layer with the magnetic field that is generated out of perpendicular magnetic field by a shear flow in the plane of the layer. Takasao et al. (2018) detected formation of MFT in 3D MHD simulations of the inner regions of the magnetized accretion discs of young stars.

Magnetic buoyancy instability develops and MFT form if the magnetic pressure is of the order of the gas pressure. Nonlinear evolution of the instability can lead to formation of the MFT with strong magnetic field (plasma $\beta<1$). For example, Machida et al. (2000) performed MHD simulations of magnetic buoyancy instability in differentially rotating magnetized disks. They reported about the formation of a filamentary-shaped intermittent magnetic structures ($\beta<1$) inside the disc. Inside the disc, magnetic buoyancy instability will lead to formation of magnetic rings of the toroidal magnetic field. Major radius of the rings will be equal to the distance to the star, while minor radius will be limited by the pressure scale height of the disc.

The MFT rise upwards in the disc under the action of buoyant force, because their density is less than the density of ambient gas. The MFT dynamics in the discs has been investigated numerically in frame of slender magnetic flux tube approximation (Sakimoto & Coroniti, 1989; Torkelsson, 1993; Chakrabarti & D’Silva, 1994; Schramkowski, 1996; Achterberg, 1996). Sakimoto & Coroniti (1989) investigated the dynamics of MFT inside the radiation pressure dominated regions of accretion discs of quasi-stellar objects taking into account aerodynamic drag, MFT shear, and heat exchange between the MFT and external gas in radiative diffusion approximation. Torkelsson (1993) considered similar problem for the case of Stokes’ drag law. Schramkowski (1996) studied MFT dynamics in optically thick radiation pressure dominated discs paying special attention to the role of shear and magnetic tension. He have found that strong shear can lead to formation of coronal loops from the initially horizontal toroidal MFT. Chakrabarti & D’Silva (1994) investigated dynamics of toroidal MFT inside the geometrically thick radiation pressure dominated discs around black holes. They have shown that magnetic tension leads to MFT collapse. Deb et al. (2017) extended model of Chakrabarti & D’Silva (1994) to take into account time-dependent evolution of the flow around the MFT.

Dynamics of slender MFT in accretion discs of young stars, and effect of external magnetic field on the dynamics of the MFT have not been investigated yet. Some of the results of simulations in slender flux tube approximation were confirmed in 2D and 3D MHD simulations of the magnetic flux escape from the accretion discs (Matsumoto et al., 1988; Shibata et al., 1990; Ziegler, 2001) and stratified gas layer (Martínez-Sykora et al., 2015).

After rising from the disc, the MFT can lead to various effects, such as outflows, variability and bursts. Dudorov (1991) proposed that MFT can be a part of the molecular outflows in the star formation regions. Chakrabarti & D’Silva (1994) and Deb et al. (2017) argued that the MFT can play a role in acceleration and collimation of jets from accretion discs of black holes. Formation of looplike magnetic structures above the discs due to magnetic buoyancy instability and magnetic reconnection can lead to burst acitivity and heating of the region above the disc (Galeev et al., 1979; Stella & Rosner, 1984; Schramkowski, 1996; Miller & Stone, 2000; Hirose & Turner, 2011; Uzdensky, 2013). Mazets & Bykov (1993) proposed that the MFT can carry away angular momentum from the accretion discs.

Present paper concerns mass and magnetic flux escape from the accretion discs of young stars due to magnetic buoyancy, formation of expanding magnetized ‘corona’ above the disc, and connection between rising MFT and IR-variability of the accretion discs. We investigate MFT dynamics in the accretion discs of young stars taking into account the radiative heat exchange between MFT and ambient gas, effects of the aerodynamic and turbulent drag, and external magnetic field. Some particular aspects of MFT dynamics in the accretion discs of young stars were considered by Dudorov & Khaibrakhmanov (2016) and Khaibrakhmanov et al. (2018). Structure of the accretion disc is calculated with the help our MHD model of the accretion discs (Dudorov & Khaibrakhmanov, 2014; Khaibrakhmanov et al., 2017).

The paper is organized as follows. In section 2.1, we discuss approximations of the model. In section 2.2, the governing equations are derived. The governing equations are written in terms of non-dimensional variables in section 2.3. Section 2.4 is devoted to solution methods of the model equations. The model of the disc is described in section 2.5. Fiducial results are presented and discussed in section 3.1. In section 3.2, we investigate influence of the model parameters on the MFT dynamics. We estimate mass loss rates due to rising MFT in section 3.3. The MFT dynamics taking into account magnetic pressure of the disc is investigated in section 3.4. We apply our model for interpretation of the observational data on IR-variability in section 3.5. Section 4 summarize results and conclusions.

2 Model

2.1 Problem statement

We consider geometrically thin, optically thick accretion disc of a young star (see Figure 1). The disc has pressure $P_{{\rm e}}$, density $\rho_{{\rm e}}$, temperature $T_{{\rm e}}$, and magnetic field with strength $B_{{\rm e}}$. Disc mass is small compared to the mass of the star $M_{\star}$, and therefore self-gravity of the disc can be neglected. We use the cylindrical system of coordinates $(r,\,\varphi,\,z)$. The vertical axis $z$ is directed along the angular velocity vector of the disc ${\bf\Omega}=(0,\,0,\,\Omega)$. The accretion disc is considered to be in hydrostatic equilibrium in the $z$-direction. Vertical coordinate of the surface of the disc is $z_{{\rm s}}$.

We assume that magnetic buoyancy instability leads to formation of a MFT in the form of a torus out of toroidal magnetic field in the disc. In the case of axial symmetry, we investigate the dynamics of the unit length cylindrical element of this magnetic torus. The MFT is located at the distance $r$ from the rotation axis at a coordinate $z_{0}$ in the beginning. It has velocity ${\bf v}$, cross-section radius $a$, pressure $P$, density $\rho$, temperature $T$, and magnetic field strength $B$. Schematic problem statement is shown in Figure 1.

The MFT is in pressure equilibrium with the external gas,

[TABLE]

Magnetic pressure outside the MFT is not considered in this equation. The effect of external magnetic field is considered in Section 3.4. We assume that temperatures inside and outside the MFT are equal to each other initially. For the ideal gas with equation of state

[TABLE]

the density difference equals

[TABLE]

where

[TABLE]

is the sound speed, $R_{{\rm g}}$ is the universal gas constant, $\mu=2.3$ is the mean molecular weight of the gas. The gas inside the MFT has smaller density comparing to the surrounding gas, as $P<P_{{\rm e}}$. Due to positive density difference $\Delta\rho>0$, the buoyant force

[TABLE]

causes rise of the MFT in the $z$-direction, where $g_{z}<0$ is the vertical component of stellar gravity. Drag force counteracts the motion of the MFT. We consider drag forces of two types, turbulent and aerodynamic drag.

2.2 Basic equations

Our model of MFT dynamics is based on the slender flux tube approximation. Similar models have been used by Sakimoto & Coroniti (1989), Torkelsson (1993), Chakrabarti & D’Silva (1994), Schramkowski (1996). We write the equations of the MFT dynamics following Dudorov & Kirillov (1986),

[TABLE]

In equation of motion (5) $\textbf{{ f}}_{{\rm d}}$ is the drag force per unit mass of the flux tube, (6, 7, 8) are the equations defining the velocity ** v**, mass $M_{{\rm l}}$ per unit length and magnetic flux $\Phi$ of the MFT, (9) is the first law of thermodynamics ($Q$ is the quantity of heat per unit mass, $U$ is the energy of MFT per unit mass, $V=1/\rho$ is the specific volume), (10) is the pressure balance equation, (11) is the equation of the hydrostatic equilibrium of disc in the $z$-direction, (1) is the equation of state, $\gamma$ is the adiabatic index.

We consider MFT motion in the $z$-direction, $\textbf{{ v}}=\left(0,\,0,\,v\right)$, $\textbf{{ f}}_{{\rm d}}=\left(0,\,0,\,f_{{\rm d}}\right)$. In this case, Equations (5, 6) are reduced to

[TABLE]

The vertical component of stellar gravity acceleration

[TABLE]

In the case of aerodynamic drag (see Parker (1979))

[TABLE]

where $C_{{\rm d}}$ is the drag coefficient $\sim 1$. The turbulent drag force can be calculated as (Pneuman & Raadu, 1972)

[TABLE]

Turbulent viscosity $\nu_{{\rm t}}$ can be estimated as (Shakura, 1972; Shakura & Sunyaev, 1973)

[TABLE]

where $\alpha$ is non-dimensional parameter characterizing the turbulence efficiency, $H$ is the scale height of the disc. Drag force is evaluated using formula (17) inside the disc, $z\leq z_{{\rm s}}$, and using formula (16) above the disc, $z>z_{{\rm s}}$.

Equations (7) and (8) give

[TABLE]

where $a_{0}$, $\rho_{0}$ and $B_{0}$ are the initial radius, density and magnetic field strength of the MFT, respectively.

Taking time derivative of Equations (9) and (10) we obtain

[TABLE]

where

[TABLE]

is the rate of heat exchange. Energy and pressure are the thermodynamic functions of density and temperature, $U=U\left(\rho,\,T\right)$ and $P=P(\rho,\,T)$. Then equations (21) and (22) can be written as

[TABLE]

where the superscript $\left(...\right)_{T}$ means derivative with respect to $T$ (with constant $\rho$) and the superscript $\left(...\right)_{\rho}$ means derivative with respect to $\rho$ (with constant $T$). Term $C_{{\rm m}}\rho$ is the derivative of magnetic pressure $B^{2}/8\pi$ with respect to $\rho$, and $C_{{\rm m}}=\dfrac{B_{0}^{2}}{4\pi\rho_{0}^{2}}$.

Solving equations (24, 25) for time derivatives of $\rho$ and $T$ and using (11), we derive equations describing the evolution of density and temperature of the MFT

[TABLE]

Now consider equations of the accretion disc structure. We solve equation of the hydrostatic equilibrium (11) using polytropic dependence of pressure on density together with equation of state (1), and get the vertical profiles of density and temperature inside the disc

[TABLE]

where $\rho_{{\rm m}}=\rho_{{\rm e}}(z=0)$, $T_{{\rm m}}=T_{{\rm e}}(z=0)$ are the density and temperature in the midplane of the disc, $k=1+1/n$, $n$ is the polytropic index, scale height $H=v_{{\rm s}}/\Omega_{{\rm k}}$,

[TABLE]

is the Keplerian angular velocity.

We set temperature $T_{{\rm a}}$ of the gas above the disc, $z>z_{{\rm s}}$, equal to the effective temperature of the disc

[TABLE]

where $L$ is the luminosity of the star. Formula (31) is derived under the assumption that the gas in the photosphere of the disc is heated by stellar radiation (Hayashi, 1981). We determine the coordinate of the disc surface from the equality of (29) and effective temperature of the disc (31)

[TABLE]

We assume that temperature is constant $T_{{\rm a}}$ above the disc surface, $z>z_{{\rm s}}$, and density falls down with $z$ according to Equation (11) to the point where $\rho_{{\rm e}}$ becomes equal to the density of molecular cloud core $\rho_{{\rm ism}}=3.8\times 10^{-20}\,{\rm g}\,{\rm cm}^{-3}$.

Heating rate of the MFT can be evaluated from heat transfer equation (see Zel’dovich & Raizer, 1967)

[TABLE]

where ** q** is the vector of the heat flux density. We consider the heat flux driven by the radiative heat conductivity (Mihalas, 1978)

[TABLE]

[TABLE]

where $\kappa_{{\rm R}}$ is the Rosseland mean opacity, $\sigma_{{\rm R}}$ is the Stefan-Boltzmann constant. The heat exchange occurs through the surface of the MFT. Let us introduce cylindrical coordinates $(r^{\prime},\,\varphi^{\prime},\,z^{\prime})$, where $r^{\prime}$ is the distance from the axis of the MFT, $\varphi^{\prime}$ is the azimuthal angle, $z^{\prime}$ is the coordinate along the axis of the MFT. Then

[TABLE]

where $q_{{\rm ext}}$ and $q_{{\rm in}}$ are external and internal heat flux densities, correspondingly. Then

[TABLE]

We determine $\kappa_{{\rm R}}$ as the power-law function of gas density and temperature following Dudorov & Khaibrakhmanov (2014).

2.3 Non-dimensional variables

Let us introduce following non-dimensional variables:

[TABLE]

where $v_{{\rm a}}$ is the Alfvén speed, $t_{{\rm A}}=H/v_{{\rm a}}$ is the Alfvén crossing time, $h_{{\rm m}}={\varepsilon}_{{\rm m}}/t_{{\rm A}}$, ${\varepsilon}_{{\rm m}}$ is the energy density of magnetic field, $B_{{\rm e}}$ is the magnetic field strength, $f_{{\rm a}}=v_{{\rm a}}/t_{{\rm A}}$. All scales are defined at the midplane of the disc.

Equations (13, 14, 26, 27, 19, 20) in the non-dimensional variables (tilde signs are omitted):

[TABLE]

where $\beta$ is the plasma beta calculated for the parameters of the disc in the midplane,

[TABLE]

[TABLE]

2.4 Method of solution and initial conditions

The system of dynamic equations (39-42) is solved with the help of the explicit Runge-Kutta method of the fourth order of accuracy. Automatic step selection is used, and adopted relative accuracy equals $10^{-4}$. At each time step, the radius and magnetic field strength of the MFT are calculated from Equations (43-44).

At the initial moment of time, the MFT with coordinates $r$ and $z_{0}$ has zero velocity $u_{0}=0$. The MFT is in thermal equilibrium with surrounding gas, $T(z_{0})=T_{{\rm e}}$. We specify the initial magnetic field strength of the MFT $B_{0}$ with the help of plasma beta

[TABLE]

where $\beta_{0}$ is the initial plasma beta inside the MFT. The initial density is determined from the condition of the pressure equilibrium (10) using plasma beta,

[TABLE]

2.5 Model of the disc

We use our MHD model of the accretion discs to calculate the structure and magnetic field of the accretion disc. Let us describe briefly the features of the model (see Dudorov & Khaibrakhmanov (2014) and Khaibrakhmanov et al. (2017) for details).

The model is MHD-generalization of Shakura & Sunyaev (1973) model. We solve MHD equations in the approximation of a geometrically thin and optically thick stationary disc. It is assumed that the turbulence is the main mechanism of the angular momentum transport in the disc. Turbulent viscosity is estimated according to expression (18). The temperature of the disc is calculated from the balance between turbulent ‘viscous’ heating and radiative cooling following Shakura & Sunyaev (1973). We use low-temperature opacities from Semenov et al. (2003). The heating of the outer parts of the disc by stellar radiation and cosmic rays is also taken into account following D’Alessio et al. (1998). These two mechanisms determine the temperature of the disc in the regions, where the turbulence can be weak. The model has two main parameters: $\alpha$ and accretion rate $\dot{M}$.

In addition to equations of Shakura & Sunyaev (1973), we solve the induction equation taking into account Ohmic diffusion, magnetic ambipolar diffusion, magnetic buoyancy and the Hall effect. Ionization fraction is determined from the equation of collisional ionization (see Spitzer, 1978) considering the ionization by cosmic rays, X-rays and radioactive decay, radiative recombinations and the recombinations on dust grains. The evaporation of dust grains and thermal ionization of hydrogen and metals are included in the model following Dudorov & Sazonov (1987).

Inner radius of the disc is assumed to be equal to the radius of stellar magnetosphere. Outer radius of the disc, $r_{\rm{out}}$, is determined as the contact boundary, where the disc pressure equals the pressure of the external medium.

3 Results

We consider accretion disc of classical T Tauri star with mass $M_{\star}=1\,M_{\sun}$, radius $R_{\star}=2\,R_{\odot}$, surface magnetic field strength $B_{\star}=2$ kG, liminocity $L_{\star}=1\,L_{\odot}$. The disc is characterized by turbulence parameter $\alpha=0.01$ and mass accretion rate $\dot{M}=10^{-7}\,M_{\sun}\,{\rm yr}^{-1}$. In this case, the inner radius of the disc lies at $r_{{\rm in}}=0.027$ au from the star, and the outer radius $r_{{\rm out}}=320$ au. Ionization fraction $x$ in the disc is calculated for dust grain radius $0.1\,\mu$m, cosmic rays ionization rate $\xi=10^{-17}\,{\rm s}^{-1}$ and attenuation length $\Sigma_{{\rm CR}}=100\,{\rm g}\,{\rm cm}^{-2}$.

In Figure 2a, we plot radial profiles of gas surface density and ionization fraction of the disc calculated using our model for adopted parameters. Surface density decreases with $r$ from $\approx 2.5\times 10^{4}\,\mathrm{g}\,\mathrm{cm}^{-2}$ at the inner boundary of the disc to $\approx 5\,\mathrm{g}\,\mathrm{cm}^{-2}$ near its outer boundary. Typical slope of $\Sigma(r)$ dependence is $-0.7$ in the region $r=[1,\,100]$ au. The radial profile of ionization fraction is non-monotonic. In the innermost part of the disk, $r<0.6$ au, the ionization fraction is high, $x>10^{-10}$, due to thermal ionization. A plateau in $x(r)$ profile at $r=[0.03,\,0.2]$ au reflects total ionization of Potassium. The ionization fraction is minimum at $r\approx 0.8$ au. At further distances, ionization fraction increases with $r$ due the decrease of surface density and corresponding more efficient ionization by cosmic rays.

In Figure 2b, we plot radial profiles of $B_{z}$ and corresponding plasma beta calculated as $\beta_{z}=8\pi\rho_{{\rm m}}v_{{\rm s}}^{2}/B_{z}^{2}$. Magnetic field strength decreases with distance $r$. Near the inner edge of the disk, $B_{z}\approx 170$ G, which is nearly equal to the stellar magnetic field at this distance. In the region of thermal ionization, $r=[0.027,\,0.6]$ au, the magnetic field is frozen into gas, and the strength of $B_{z}$ is proportional to gas surface density. Plasma beta decreases from 300 to 30 in this region. Magnetic field strength abruptly decreases by two orders of magnitude at $r\simeq 0.6$ au. This is the transition to the ‘dead’ zone, region of low ionization fraction where Ohmic and ambipolar diffusion prevent magnetic field amplification. The ‘dead’ zone occupies region from 0.6 to 30 au for considered parameters. Plasma beta is $10^{3}-10^{5}$ inside the ‘dead’ zone. Further $B_{z}$ and $\beta_{z}$ decrease with $r$ to the values $5\times 10^{-3}$ G and $6$, respectively, near the outer edge of the disc.

In the following we consider dynamics of the MFT that can form in the inner region of the disc, $r<0.6$ au. We carry out simulations of the MFT dynamics for the various initial radii $a_{0}=\left[0.01,\,0.1,\,0.2,\,0.4\right]\,H$, plasma betas $\beta_{0}=[0.01$, $0.1$, $1$, $10]$, coordinates $r$ in range $0.027\div 0.6$ au and coordinates $z_{0}=[0.5,\,1]\,H$. Adiabatic index of the gas $\gamma=7/5$.

Density, temperature and magnetic field strength of the accretion disc are listed in Table 1 (column 1: $r$-distance, column 2: midplane density, column 3: midplane temperature, column 4: effective temperature of the disc, column 5: scale height, column 6: magnetic field strength $B_{{\rm e}}=B_{z}$ in the disc, column 7: $z$-coordinate of the surface of the disc). Table 1 shows that density, temperature and magnetic field strength decrease with $r$-distance in the disc. Scale height of the disc increases with $r$. We choose polytropic index $n=3$. In this case, the coordinate of the disc surface varies from $z_{{\rm s}}=2.28\,H\approx 0.0014$ au at $r=0.027$ au to $z_{{\rm s}}=2.38\,H\approx 0.08$ au at $r=0.6$ au.

3.1 Fiducial run

In this section, we present and discuss the simulations of MFT dynamics for the following representative parameters: $r=0.15$ au, $\beta_{0}=1$, $a_{0}=0.1\,H$, and $z_{0}=0.5\,H$. We perform two sets of simulations to study the role of radiative heat exchange. First, we simulate the dynamics of MFT evolving in thermal equilibrium with ambient gas (Section 3.1.1). In this case, Equations (41) and (42) are excluded, equality $T=T_{{\rm e}}$ is adopted, and density is determined from pressure balance Equation (10). Second, we carry out the simulations taking into account the radiative heat exchange according to Equations (41) and (42) (Section 3.1.2).

3.1.1 Thermal equilibrium

In this section, we discuss general features of the MFT dynamics in the case of thermal equilibrium. In Figures 3(a-c), we plot the dependences of velocity, temperature and radius of the MFT on its $z$-coordinate. Corresponding dependences of $z$-coordinate, drag and buoyant forces, internal end external densities on time are depicted in Figures 3(d-f). Absolute value of the drag force is plotted for convenience.

Figure 3(a) shows the MFT begins to rise with high acceleration. The velocity of the MFT grows very fast from 0 to $0.5\,\mathrm{km}\,\mathrm{s}^{-1}$, because the buoyant force $f_{{\rm b}}$ is much stronger than the drag force $f_{{\rm d}}$ in the beginning of motion (see Figure 3(e)). After that, the buoyancy and drag forces become nearly equal to each other, $|f_{{\rm d}}|\lesssim f_{{\rm b}}$, and the velocity monotonically increases up to $\approx 1.8\,\mathrm{km}\,\mathrm{s}^{-1}$ at the surface of the disc $z_{{\rm s}}\approx 2.27\,H$. The MFT rises to the surface during time $\approx 0.8\,t_{{\rm A}}\approx 18$ d, as Figure 3(d) shows. Absolute values of the buoyant and drag forces become small, acceleration of the MFT vanishes and it acquires nearly steady velocity $v\approx 1.8\,\mathrm{km}\,\mathrm{s}^{-1}$ further.

The MFT is in thermal equilibrium with ambient gas, so that $T=T_{{\rm e}}$ during its motion, as Figure 3(b) demonstrates. Figure 3(c) shows that the MFT expands in the course of the rise, i.e. its radius increases with $z$. Correspondingly, density of the MFT decreases, as Figure 3(f) shows. Internal density always stays less than the external one in thermal equilibrium. Ultimately radius of the MFT exceeds the thickness of the disc, $a>z_{{\rm s}}$, at $z\approx 2.8\,H$. We do not simulate further motion of the MFT, as the slender tube approximation violates under such circumstances. The process of substantial expansion of MFT above the disc can be interpreted as a transformation of rising slender flux tubes into non-uniform expanding magnetized corona of the disc. Further we will call this process as a dispersal of the MFT. Effect of magnetic flux escape from the disc with subsequent formation of magnetized corona has been found by Miller & Stone (2000); Machida et al. (2000); Johansen & Levin (2008); Turner et al. (2010); Romanova et al. (2011); Takasao et al. (2018) in the MHD simulations of the accretion discs. Takasao et al. (2018) found formation of the flux tubes from the regular magnetic field of the disc, that confirms our assumptions.

3.1.2 Role of heat exchange

Consider dynamics of the MFT in the case, when radiative heat exchange is taken into account self-consistently (Figure 4). Dynamics of the MFT inside the disc, $z<z_{{\rm s}}$ ($t<1\,t_{{\rm A}}$), is similar to the one discussed in Section 3.1.1 (Figure 3). The MFT rises with increasing speed, expands and its density decreases. It accelerates to $v\approx 1.8\,{\rm km}\,{\rm s}^{-1}$ near the surface of the disc. Temperature $T$ is nearly equal to $T_{{\rm e}}$, and its vertical profile matches the polytropic profile (29). Approximate equality of the temperatures inside and outside the MFT reflects fast heat exchange at $z<z_{{\rm s}}$.

Figure 4(a) shows that the MFT starts to decelerate after rising from the disc. The velocity goes to zero, when the MFT reaches point $z\approx 2.7\,H$ at $t\approx 1\,t_{{\rm A}}$. After that, the velocity becomes negative, and the MFT starts to move downwards, i.e. it ‘sinks’. When the MFT reaches point $z\approx 2.2\,H$, its velocity changes sign again, and the MFT starts to move upwards in the second time. Such an up and down motion is observed within the time period from $t\approx 1\,t_{{\rm A}}$ to $t\approx 2\,t_{{\rm A}}$. Velocity periodically changes its sign, while its absolute value decreases at this part of the trajectory. Further, at $t>2\,t_{{\rm A}}$, the MFT rises from the disc and continues monotonic upward motion with small nearly steady velocity of about $0.04\,{\rm km}\,{\rm s}^{-1}$.

The oscillatory motion of the MFT near the surface of the disc in the time interval from $t\approx 1\,t_{{\rm A}}$ to $t\approx 2\,t_{{\rm A}}$ is explained by the following. The internal temperature of the MFT decreases with respect to constant temperature outside, when the MFT rises from the disc (see Figure 4(b)). The temperature takes minimum value $T_{{\rm min}}\approx 380$ K, when the MFT rise to the point $z\approx 2.7\,H$. The cooling of the MFT is caused by its practically adiabatic expansion. The radiative heat exchange appears to be inefficient to compensate adiabatic cooling at this part of the trajectory. Decrease of internal temperature with respect to the external one causes the MFT to expand more slowly than in the case of thermal equilibrium. The density of the MFT decreases with $z$ more slowly than the external density, and $\rho$ exceeds $\rho_{{\rm e}}$, i.e. the MFT loses buoyancy at the point $z\approx 2.4\,H$. The MFT moves by inertia upwards for some time, until it stops at $z\approx 2.7\,H$. Then the MFT starts to ‘sink’ due to negative buoyancy. The MFT contracts a little, and its temperature and density grow in a process of downward motion (see Figures 4(b, c, f)). The buoyancy restores ($\rho<\rho_{{\rm e}}$), when the internal and external temperatures become nearly equal to each other. This leads to deceleration of the MFT and to change of the velocity sign near the point $z=2.2\,H$. The radiative heat exchange leads to equalization of internal and external temperatures in a process of further periods of up and down motion with respect to the point of zero buoyancy ($\rho=\rho_{{\rm e}}$). Ultimately, temperatures $T$ and $T_{{\rm e}}$ become nearly equal to each other, and the MFT oscillations change onto monotonic upward motion with steady velocity above the disc surface at $t>2\,t_{{\rm A}}$. The oscillations decay due to equalization of the temperatures, and corresponding decrease of density difference and buoyant force. The MFT is optically thick during its motion. We call the oscillations discussed in this section as the thermal ones.

3.2 Dependence on parameters

We investigate the dynamics of MFT for various initial radii $a_{0}$, distance $r$ and plasma beta $\beta_{0}$ in this section. As it was mentioned in the introduction, MFT are likely form with $\beta_{0}\sim 1$. In this section we consider MFT dynamics for plasma beta in range $[0.01,\,10]$ to study the dependence of MFT characteristics on the initial magnetic field strength. Khaibrakhmanov et al. (2018) have shown that the characteristics of MFT near disc surface practically do not depend on the initial position $z_{0}$. In this section, we present calculations for $z_{0}=0.5H$. Radiative heat exchange is taken into account in all considered runs. Internal and external temperatures are equal to each other initially.

In Figure 5, we plot the dependences of MFT $z$-coordinate and temperature on time, and velocity versus $z$-coordinate , for $r=0.15$ au, $\beta_{0}=1$, $z_{0}=0.5\,H$ and various initial radii $a_{0}$. Black lines correspond to the fiducial case.

Figure 5(a) shows that rise time increases with decreasing radius of MFT. For example, the MFT with $a_{0}=0.01\,H$ (blue line) rises to the surface of the disc over time of $7\,t_{{\rm A}}$, while the MFT with $a_{0}=0.4\,H$ (magenta line) reaches disc surface within time of $0.4\,t_{{\rm A}}$ ($t_{{\rm A}}$ is the Alfvén crossing time). MFT with larger radius moves faster, because the buoyant force is proportional to MFT volume ($\propto a^{3}$), while the drag force is proportional to MFT surface area ($\propto a^{2}$) (see Figure 5(c)).

MFT with initial radius $a_{0}=0.1\,H$ in Figure 5 exhibits thermal oscillations. Our simulations show that the MFT with initial radius $a_{0}>0.17\,H$ (see green and magenta lines in Figure 5) disperse fast above the disc before reaching the point of zero buoyancy. The MFT with initial radius $a_{0}<0.06\,H$ (see orange and blue lines in Figure 5) move slowly, so that the radiative heat exchange effectively equalizes internal and external temperatures (see Figure 5(b)), preventing the loss of buoyancy.

In Figure 6, we plot dependence of MFT $z$-coordinate on time for $a_{0}=0.1\,H$, $\beta_{0}=1$, $z_{0}=0.5\,H$ and various $r$-coordinates. Black lines correspond to the fiducial run discussed in Section 3.1.2.

Figure 6 shows that time of MFT rise to the surface of the disc increases with $r$-distance. The MFT float up to the surface over time of $2$ d at $r=0.027$ au, and $140$ d at $r=0.6$ au, and velocity of the MFT decreases from $2.8\,{\rm km}\,{\rm s}^{-1}$ to $1.5\,{\rm km}\,{\rm s}^{-1}$.

MFT experience the thermal oscillations only in the region $r\leq 0.2$ au (blue and black lines in Figure 6). Oscillation period increases with $r$-distance. It is approximately equal to $1$ d at the distance $r=0.027$ au (blue line), and $10$ d at the distance $r=0.15$ au (black line). The simulations of MFT dynamics for various initial plasma beta show that only the MFT with $\beta_{0}=1$ experience thermal oscillations in the region $r\leq 0.2$ au.

In the region $r>0.2$ au, the MFT rise so slowly that the radiative heat exchange is able to equalize internal and external temperatures, preventing the loss of buoyancy.

In Table 2, we present initial parameters and some results of the simulations of the MFT dynamics at the distances $r=0.027$ and $r=0.6$ au. Initial coordinate $z_{0}$ is $0.5\,H$ in all presented runs. Column 1 shows number of run. Column 2 contains the values of $r$-distance. Initial plasma beta $\beta_{0}$ and radius of the MFT $a_{0}$ are listed in columns 3 and 4, respectively. The MFT mass, $M$, and magnetic flux, $\Phi$, are given in columns 5 and 6, respectively. We list velocity, $v_{{\rm surf}}$, radius, $a_{{\rm surf}}$, and magnetic energy, $E_{{\rm m}}$, of the MFT in the moment, when it crosses the disc surface, in columns 7, 8 and 9. Time of rise to the surface of the disc, $t_{{\rm surf}}$, is given in column 10. Column 11 gives time of the toroidal magnetic field generation (see discussion in Sections 3.3 below). Column 12 contains the rate of mass loss due to buoyancy (see discussion in Section 3.3 below). In column 13, we list ratios of the minimal MFT temperature during its rise with respect to corresponding external temperature. Typically, the minimal temperature is achieved at $z\approx 2.3-3\,H$.

Table 2 shows that masses of MFT range from $6.3\times 10^{-13}\,M_{\sun}$ (run 12) to $1.0\times 10^{-7}\,M_{\sun}$ (run 10). The MFT mass increases with its radius. The density and, therefore, MFT mass increase with $\beta_{0}$ (it follows from Equation (2)). The MFT transfer magnetic fluxes in range from $2.3\times 10^{19}$ Mx (run 15) to $1.1\times 10^{24}$ Mx (run 9). Velocity of the MFT with $\beta_{0}=1$ has the values in interval $0.2-4.3\,{\rm km}\,{\rm s}^{-1}$ that is comparable with local sound speed of $1-2.5\,{\rm km}\,{\rm s}^{-1}$. The MFT with smaller $\beta_{0}$ accelerate to higher velocities and can have supersonic speed. Maximum velocity $\sim 15.6\,{\rm km}\,{\rm s}^{-1}$ is achieved in run 20. In this case, bow shocks will probably form during motion of the MFT, which can change drag law and lead to additional heating. In the following sections we will present and discuss the simulations, in which the velocity of the MFT does not exceed significantly the sound speed.

When the MFT crosses the disc surface, its radius is several times smaller than $z_{{\rm s}}$. For example, $a(z=z_{{\rm s}})=0.035\,H$ for run 1. The MFT continue to expand, when they move above the disc (see Figure 3(c)). The MFT radius becomes larger than the disc height at $z\sim 2.7-3.5\,H$ in performed runs. For instance, the MFT radius exceeds the disc height at $z=2.7\,H$ in run 8, and at $z=3.5\,H$ in run 9. As it was stated in Section 3.1, we interpret the dispersal of the rising MFT as a formation of expanding magnetized corona above the surface of the disc.

Thinnest MFT with $a_{0}\leq 0.01\,H$ stay in thermal equilibrium with external gas, $T\approx T_{{\rm e}}$, as column 13 of Table 2 shows. Thicker MFT cool down in comparison to the external gas during their motion. For example, the MFT with $\beta_{0}=1$, $a_{0}=0.2$ and $r=0.6$ au cools down to minimal temperature $\approx 0.58\,T_{{\rm e}}$ (run 11). Among the runs presented in Table 2, prominent thermal oscillations are found only in runs 18 and 19, in agreement with the discussion of the Figure 6. Oscillation behaviour is also observed in run 8 with $\beta_{0}=10$, $a_{0}=0.2$, $r=0.6$, but the oscillations rapidly decay within three periods in this case. In the other runs, the MFT disperse fast above the disc.

3.3 Mass and magnetic flux loss due to buoyancy of flux tubes

As it has been shown in sections 3.1 and 3.2, the MFT rise from the disc to its atmosphere carrying away mass and magnetic flux. The rising MFT can be the seed for the formation of jets and outflows from accretion discs. Similar idea was considered by Chakrabarti & D’Silva (1994) and Deb et al. (2017) in application to the accretion discs around black holes. The rate of vertical mass transport via buoyancy can be estimated as a mass of the flux tubes rising from the disc per unit of time. We consider times of the MFT formation and rise as characteristic time scales. The characteristic time of MFT formation is comparable to the time of the toroidal magnetic field generation. To estimate efficiency of the mass transport form disk interior to its atmosphere via buoyancy, we compare the characteristic times in this section.

The time of the azimuthal magnetic field generation can be estimated from the $\varphi$-component of induction equation in the approximations of the accretion disc model,

[TABLE]

This equation shows that $B_{\varphi}$ is generated from $B_{z}$ by the differential rotation of the disc. The velocity $v_{\varphi}$ is determined from the balance between the gravity and centrifugal force in the $r$-direction. In our case,

[TABLE]

so Equation (49) in the case $z^{2}/r^{2}\ll 1$ turns to

[TABLE]

Therefore, the time of the azimuthal magnetic field generation up to a given value $B_{\varphi}$

[TABLE]

where $P_{{\rm k}}$ is the Keplerian period. Typical time scale of $B_{\varphi}$ amplification up to the value $B_{\varphi}=B_{z}$ is nearly two Keplerian periods for $z=0.05\,r$. Keplerian period increases with distance as $P_{{\rm k}}\propto r^{3/2}$, therefore fastest generation of $B_{\varphi}$ takes place in the innermost region of the disc.

MFT rise and generation times for considered parameters are listed in columns 10 and 11 of Table 2, respectively. We plot dependence of $P_{{\rm k}}$, $t_{{\rm gen}}$ and $t_{{\rm surf}}$ on the $r$-coordinate in Figure 7(a). Table 2 and Figure 7(a) show that characteristic times increase with distance. The rise time $t_{{\rm surf}}$ increases from 0.5 d at $0.027$ au to $200$ d at $r=0.6$ au. It is nearly equal to Keplerian period and order of magnitude less than generation time $t_{{\rm gen}}$ (52). This means that the dynamics of the MFT will occur in two stages. At the first stage having duration $t_{{\rm gen}}$, toroidal magnetic field is generated and MFT form due to magnetic buoyancy instability. At the second stage, the MFT rise from the disc over the time $t_{{\rm surf}}$ and carry away some mass and magnetic flux to disc atmosphere. After that, this two-stage process repeats. Since $t_{{\rm surf}}\ll t_{{\rm gen}}$, the process of mass and magnetic flux transport from disc to its atmosphere is periodic with typical period $t_{{\rm gen}}$. Table 2 and Figure 7(a) show that $t_{{\rm gen}}\approx(0.5-0.7)$ yr at $r=0.027$ au and $t_{{\rm gen}}\approx(8-10)$ yr at $r=0.6$ au.

The MFT can form inside the region $r=0.027\div 0.6$ au. Results presented in Figure 7(a) indicate that MFT will form inside the disc and rise from it first at small $r$ and then at farther distances from the star.

The rate of vertical mass transport via the buoyancy, $\dot{M}_{{\rm b}}$, can be estimated by division of the MFT mass by the characteristic time of the magnetic field amplification $t_{{\rm gen}}$. Values of $\dot{M}_{{\rm b}}$ lie in range $\sim 10^{-12}-10^{-7}\,M_{\sun}\,{\rm yr}^{-1}$ (column 12 of Table 2). This value increases with plasma beta and radius of MFT. For example, maximum rate $\dot{M}_{{\rm b}}=7.9\times 10^{-8}\,M_{\sun}\,{\rm yr}^{-1}$ is found in run 8 with $a_{0}=0.1\,H$ and $\beta_{0}=10$. Rate of magnetic flux transport via buoyancy can be estimated in similar way, $\dot{\Phi}_{{\rm b}}=\Phi/t_{{\rm surf}}$.

In Figure 7(b), we plot the dependences of $\dot{M}_{{\rm b}}$ and $\dot{\Phi}_{{\rm b}}$ on the $r$-coordinate for the MFT with initial parameters $a_{0}=0.1\,H$, $\beta_{0}=1$, $z_{0}=0.5\,H$. Figure 7(b) shows that $\dot{M}_{{\rm b}}$ increases with distance, from $\approx 6\times 10^{-9}\,M_{\odot}\,{\rm yr}^{-1}$ at $r=0.027$ au to $\approx 2\times 10^{-8}\,M_{\odot}\,{\rm yr}^{-1}$ at $r=0.6$ au. Therefore, the periodic process of mass and magnetic flux transport from the disc to its atmosphere caused by magnetic buoyancy is the most efficient near the outer edge of thermal ionization zone, where the MFT with large radius form. The rate of vertical mass transport via the buoyancy at $r=0.6$ au is five times smaller than the mass accretion rate in the disc $\dot{M}$, i.e. $20\%$ of the accreted mass can be transported from the disc to its atmosphere via buoyancy.

It should be noted that ‘dead’ zone is situated in the region between $r=0.6$ au and $r=33$ au for the considered parameters of the disc. The rate $\dot{M}_{{\rm b}}$ will decrease rapidly with $r$ beyond $r=0.6$ au, since amplification of the toroidal magnetic field is hindered by Ohmic diffusion and MFT cannot form inside the ‘dead’ zone. MFT can form only in the surface layer of the disc above the ‘dead’ zone. Typical surface density of this layer is $5-10\,{\rm g\,cm}^{-2}$ (e.g., Dudorov & Khaibrakhmanov (2014)), and the MFT forming inside this layer have small radius of $\sim 0.01\,H$ and carry away tiny mass.

Figure 7 shows that rate of magnetic flux transport via buoyancy also increases with $r$, from $\dot{\Phi}_{{\rm b}}\approx 1\times 10^{22}\,{\rm Mx\,yr}^{-1}$ at $r=0.027$ au to $\dot{\Phi}_{{\rm b}}\approx 2.7\times 10^{22}\,{\rm Mx\,yr}^{-1}$ at $r=0.6$ au. Minimum and maximum magnetic fluxes are $2.3\times 10^{19}$ and $1.1\times 10^{24}$ Mx (runs 15 and 9), and magnetic energies are $3.0\times 10^{33}$ and $3\times 10^{37}\,{\rm erg}$ (see columns 6 and 9 of Table 2). Total magnetic flux of the disc equals $5\times 10^{29}$ Mx for considered parameters. Therefore, nearly $20\%$ of the disc magnetic flux can be lost via the magnetic buoyancy over $\sim 1$ Myr.

3.4 Effect of external magnetic field

In sections 3.1-3.3, we investigated the dynamics of the MFT in the disc without external magnetic field. The magnetic field of the disc can influence the MFT dynamics mainly through the magnetic pressure. In this section, we discuss the MFT dynamics in the case, when the magnetic pressure outside the MFT is taken into account in pressure equilibrium Equation (10), i.e.

[TABLE]

First, we discuss the dynamics of the MFT for the fiducial parameters $r=0.15$ au, $\beta_{0}=1$, $a_{0}=0.1H$, $z_{0}=0.5\,H$ (see Section 3.1.1). In Figure 8, we plot the vertical profiles of velocity, temperature and radius of the MFT (panels (a), (b) and (c), respectively), as well as dependence of its $z$-coordinate, drag and buoyant forces, internal and external densities on time (panels (d), (e) and (f), respectively). Insets in panels (a-c) show dependences of $v$, $T$ and $a$ on time. Simulation is performed for the case of thermal equilibrium to eliminate effects of thermal oscillations.

Figure 8(a) shows that the MFT floats with increasing velocity at the initial part of rise, $z<2\,H$ ($t<0.75\,t_{{\rm A}}$). Velocity of the MFT reaches the value of $1.35\,{\rm km}\,{\rm s}^{-1}$ at $z=2\,H$. After that the MFT decelerates, and its velocity goes to zero at $z=2.35\,H$. Then the velocity becomes negative, and the MFT starts to move downwards. Its velocity changes sign again at $z=2\,H$, and the MFT starts to move upwards. In the following, the MFT experiences oscillatory motion up and down with respect to the point $z\approx 2.15\,H$. Maximum value of the velocity decreases during these oscillations, i.e. the oscillations decay. Figures 8(b, c) show that the MFT pulsates in the process of the oscillations, i.e. its temperature and radius periodically change with respect to the values of $860$ K and $0.25\,H$, respectively. Period of the oscillations, $P_{{\rm osc}}$, equals $0.4\,t_{{\rm A}}$.

The oscillations of the MFT are explained as follows. Figures 8(e, f) demonstrate that at the initial part of the trajectory, $t<0.75\,t_{{\rm A}}$, the MFT density is less than the external density, and the buoyant force is positive. Figures 8(c, d) show that the MFT expands and its magnetic field weakens. Strength of the external magnetic field, $B_{{\rm e}}$, is assumed to be constant in our simulations. Internal magnetic field strength $B$ becomes less than $B_{{\rm e}}$ at the point $z=2.18\,H$ ($t=0.8\,t_{{\rm A}}$). Increasing magnetic pressure outside the MFT causes it to expand slower than in the case of rise without external magnetic field (see Figure 3(f)). Density of the MFT is greater than the external density, and the buoyant force is negative, since $B$ is less than $B_{{\rm e}}$ in this point. The MFT rise to the point $z\approx 2.35\,H$ by inertia, and then start to move downwards. The MFT contracts during the downward motion, i.e. its radius decreases, while the temperature, density and magnetic field strength increase. The buoyancy become positive again, and the MFT starts to move upwards after $B$ exceeds $B_{{\rm e}}$. In the following, the MFT oscillates near the point of zero buoyancy, that is determined by equality $B=B_{{\rm e}}$ in this case. We call these oscillations as magnetic ones.

Drag force reduces the kinetic energy of the MFT, and leads to the decay of oscillations. Magnitude of the velocity oscillations reduces by factor of ten over five periods of oscillations. The MFT radius approaches the value of $0.25\,H$ in this case.

The magnetic oscillation has one significant difference form the thermal ones discussed in Section 3.1.2. The thermal oscillations are caused by departure from thermal equilibrium. Therefore, they cease, and the MFT start to move upwards monotonically after the temperatures inside and outside the MFT become equal to each other. The magnetic oscillations are caused by the action of the external magnetic pressure. In the considered case, constant $B_{{\rm e}}$ prevents further rise of the MFT. Dynamics of the MFT is characterized by decaying oscillations near the point of zero buoyancy. In the case, when external magnetic field $B_{{\rm e}}$ decreases with height, the MFT would float farther.

We performed the simulation for the same parameters as discussed above, but taking into account the radiative heat exchange. The simulations show that the picture of the MFT dynamics is practically the same, as in Figure 8, i.e. departure from thermal equilibrium does not influence the magnetic oscillations.

Figure 9 demonstrates dynamics of the MFT taking into account external magnetic field for various initial radii and plasma betas at $r=0.2$ au.

Figure 9 shows that the magnetic oscillations are observed for all considered parameters. The MFT with stronger magnetic field (small $\beta_{0}$) experience oscillations at higher altitude. For example, the MFT with $\beta_{0}=0.1$ and $a_{0}=0.2\,H$ (orange lines) oscillate above the surface of the disc, near the point $z=2.5\,H$. Oscillation period in this case is less than in the case $\beta_{0}=1$ (black lines), while magnitude of velocity oscillations is higher, and magnitude of radius variations is smaller.

The MFT with a smaller initial radius move slower, but oscillate at about the same height as the tubes of a larger radius at the same initial field strength (compare black and green lines in Figure 9).

We carried out the simulations of MFT dynamics taking into account external magnetic pressure for a full set of parameters $a_{0}$ and $\beta_{0}$ considered in Table 2 and for various $r$-coordinates in range $[0.027,\,0.6]$ au. The simulations show that the altitude, at which MFT oscillate, increases with initial magnetic field strength of the MFT. The oscillations of the MFT with $\beta_{0}=1$ take place under the surface of the disc, while the MFT with $\beta_{0}<1$ oscillate above the surface. Typical radii of the MFT are of order of $0.5\,H$, so that the oscillations of the MFT will lead to periodical changes of the disc structure near its surface, in the region from $z\approx 1.5\,H$ to $z\approx 3\,H$.

In Figure 10, we plot the periods of MFT oscillations $P_{{\rm o}}$ for various $r$-coordinates. The results are obtained for the MFT with initial radius $a_{0}=0.1\,H$, coordinate $z_{0}=0.5\,H$, plasma beta $\beta_{0}=1$ (black line) and $\beta_{0}=0.1$ (grey line). Figure 10 shows that oscillation period increases with the radial distance from $1$ d $\approx 0.6\,P_{{\rm k}}$ at $r=0.027$ au to $100$ d $\approx 0.6\,P_{{\rm k}}$ at $r=0.6$ au in the case $\beta_{0}=1$. The oscillation period of the MFT with $\beta_{0}=0.1$ are smaller and vary form $0.5$ d at $r=0.027$ au to $45$ d at $r=0.6$ au.

Comparison of toroidal magnetic field generation and MFT rise times depicted in Figure 7(a) and oscillation periods in Figure 10 shows that $P_{{\rm osc}}\ll t_{{\rm gen}}$. We conclude that the dynamics of the toroidal magnetic field in the considered region proceeds in two stages: slow generation of the toroidal magnetic field ($t_{{\rm gen}}=0.5-10$ yr) with subsequent fast rise and oscillations of the MFT ($P_{{\rm osc}}=1-100$ d).

3.5 Comparison with observations

As it has been shown in Sections 3.1.2 and 3.4, thin MFT oscillate near the surface of the disc after rising from interior. Density, radius and temperature periodically change during the oscillations. The oscillations can lead to variability of the accretion disc radiation. MFT contain both gas and refractory dust particles in the considered region. Therefore, the oscillations can cause the IR-variability of the disc radiation. Based on the resistive MHD simulations of the inner region dynamics of the minimum mass solar nebula, Turner et al. (2010) proposed that magnetic activity can lift dust grains into the disc atmosphere and cause the IR-variability of YSO (see also discussion in Flaherty et al. (2011)). In this paper, we study similar effect in the frame of slender magnetic flux tube approximation.

IR-variability has been observed in many YSO. Flaherty et al. (2016) have found that the stars in the Chameleon I cluster exhibit variability on time-scales of months (20-200 days). Magnitude of the fluctuations ranges from 0.05 to 0.5 mag. In order to test the hypothesis that rising MFT cause the IR-variability of YSO, we perform the simulations of MFT dynamics in two discs of classical T Tauri stars from the sample presented in Flaherty et al. (2016). We have chosen stars J11092266-7634320 and J11100369-7633291 (indexed as 439 and 530, respectively, in Table 1 from Flaherty et al. (2016)). Masses, accretion rates, luminosities and effective temperatures of these stars are given in columns 2-5 of Table 3. Stellar radii are estimated from the relation $L_{\star}=\sigma_{{\rm R}}T_{{\rm eff}}^{4}4\pi R_{\star}^{2}$ (column 6), where $T_{{\rm eff}}$ is the stellar effective temperature, $R_{\star}$ is the stellar radius. The periods of IR-variations, $\Delta t$, measured by Flaherty et al. (2016) are listed in column 7. Column 8 gives values of considered $r$-distances from the star. For each star, we simulate the MFT dynamics at two distances: $r_{{\rm i}}$ (values with symbol ‘i’ in brackets in column 8), where temperature is nearly equal to temperature of silicate dust grains evaporation $\sim 1500$ K (see Pollack et al., 1994), and $r_{{\rm o}}$ (values with symbol ‘o’ in brackets in column 8), determining the outer boundary of the region of thermal ionization. Radii $r_{{\rm i}}$ and $r_{{\rm o}}$ bound the region of the efficient toroidal magnetic field amplification, where formation of the MFT with both gas and dust is possible. We take magnetic field strength at the stellar surface to be 2 kG for both stars (see Yang & Johns-Krull, 2011). Using these stellar parameters, we calculate the structure of the accretion discs using accretion disc model of Dudorov & Khaibrakhmanov (2014). Midplane density, temperature, scale height and magnetic field strength of the discs of stars 439 and 530 are presented in columns 9-12 of Table 3.

The simulations are carried out for $\beta_{0}=[0.1,1]$ and $a_{0}=0.01\,H\div 0.4\,H$. The simulations show that oscillation periods $P$ weakly depend on the MFT radius. As an example, in Figure 11 we plot dependences of $z$-coordinate , radius $a$ and temperature $T$ on time for the MFT in the disc of star 439, at $r_{{\rm i}}$ and $r_{{\rm o}}$. Initial MFT parameters: $z_{0}=0.5\,H$, $a_{0}=0.1\,H$, $\beta_{0}=1$. Figure 11 shows that dynamics of the MFT is similar to that considered in Section 3.4. The MFT rise and then oscillate near the surface of the disc. Oscillation period at $r_{{\rm i}}$ (20 d) is less than the oscillation period at $r_{{\rm o}}$ (70 d). The MFT experience temperature pulsations around value $T\sim 700$ K at $r_{{\rm i}}$ and around $T\sim 500$ K at $r_{{\rm o}}$. With these temperatures, maximum of emission peaks at wavelengths $\lambda\approx(4-6)\,\mu$m, according to Wien’s displacement law. Maximal temperature variations are $100$ K at $r_{{\rm i}}$ and $200$ K at $r_{{\rm o}}$.

The dependences of oscillation periods $P_{{\rm osc}}$ on the $r$-distance for each star are shown in Figure 12. The results are obtained for $a_{0}=0.1\,H$, $z_{0}=0.5\,H$, plasma beta $\beta_{0}=1$ (star 439) and $\beta_{0}=5$ (star 530). Star 530 has small accretion rate, and the region of thermal ionization is situated close to the star, $r<0.1$ au. The disc of star 439 has more extended region of thermal ionization, $r<0.42$ au. In both cases, the oscillation periods increase with distance from the star.

For the adopted MFT parameters, the range of oscillation periods is in agreement with the observational values of variability times. Observed period for star 439, $P_{439}=32$ d corresponds to oscillation of the MFT at $r\approx 0.25$ au for $\beta_{0}=1$. Period of star 530, $P_{530}=35$ d, corresponds to the MFT oscillating at $r=0.1$ au for $\beta_{0}=5$. Thus, our simulations confirm the hypothesis that rising magnetic fields can be a source of the IR-variability of YSO.

To determine whether the MFT can contribute to IR radiation of considered T Tauri stars or not, we estimate optical depth of the MFT with respect to IR radiation as $\tau=2a\rho\kappa$, where opacity $\kappa$ is determined as a function of density and temperature according to Dudorov & Khaibrakhmanov (2014). We get that $\tau$ has the values between 20 and 70 during oscillations at $r=r_{{\rm o}}$ for star 439. Optical thickness for star 530 is of the same order. We conclude that the MFT are very optically thick.

There may be two possible processes causing the variability. First, MFT temperature fluctuations (like in Figure 11) can cause the radiation variability. Second, the oscillating MFT can intercept stellar radiation and periodically cast a shadow on the outer disc regions. Stellar radiation reprocessed by the dust in outer disc regions will have an oscillating IR-component.

4 Conclusion and discussion

We investigated numerically the dynamics of slender MFT in the accretion discs of T Tauri stars. We considered the MFT forming in the region of intense generation of the toroidal magnetic field due to magnetic buoyancy instability. We formulate the equations of slender MFT dynamics taking into account the aerodynamic and turbulent drags. Khaibrakhmanov et al. (2018) have investigated the MFT dynamics taking into account the heat exchange with the external medium of constant temperature. The radiative flux has been estimated in the diffusion approximation. In present work, we investigated additionally influence of the magnetic field of the disc on the MFT dynamics. The vertical structure of the disc was calculated using the polytropic equation of state. In particular, the thermal and magnetic oscillations of the MFT was investigated in this paper.

Simulations were performed for initial radii of the MFT $a_{0}=0.01\,H\div 0.4\,H$ and radial distances $r=0.027\div 0.6$ au in the accretion disc of solar mass T Tauri star. Slender flux tube approximation allows us to investigate the dynamics of ‘soft’ MFT with weak magnetic field ($\beta\ga 1$) and ‘stiff’ MFT with strong magnetic field ($\beta<1$). Latter case cannot be considered using multidimensional numerical simulations in frame of classical MHD because of strict limitations on the time step. For comparison, we investigated dynamics of the MFT with plasma betas $\beta_{0}=0.01\div 10$. The structure of the accretion disc was calculated with the help of the model of Dudorov & Khaibrakhmanov (2014).

First of all, we considered the dynamics of the MFT in absence of external magnetic field. In this case, we found two regimes of the MFT dynamics. Thin MFT with initial radius $a_{0}\leq 0.05\,H$ and thick MFT with $a_{0}\geq 0.16\,H$ rapidly accelerate, then rise with slowly increasing velocity, decelerate a little and acquire nearly steady velocity above the surface of the disc. The MFT of intermediate radii $a_{0}\sim 0.1\,H$ experience thermal oscillations during some time after rising from the disc. The oscillations are due to adiabaticity and slow heat exchange with the external gas. These oscillations are found only for the MFT with $\beta_{0}=1$ formed at $r$-distances less than $0.2$ au. After radiative heat exchange equalizes internal and external temperatures, the oscillations decay and the MFT continue to move upwards above the disc.

We studied dependence of the MFT velocity, mass and magnetic flux on its initial parameters. Velocity of the MFT increases with initial radius and magnetic field strength. Typical rise velocities are of several km s*-1*. The MFT with weak magnetic field ($\beta_{0}=10$) have velocities of $0.05-0.5\,{\rm km}\,{\rm s}^{-1}$. The MFT with plasma $\beta=1$ accelerate to velocity $\approx 0.2-4$ km s*-1* comparable with sound speed, in agreement with analytical estimates of Khaibrakhmanov et al. (2018). The MFT with strong initial magnetic field (plasma beta $<1$) can reach supersonic velocity up to $\sim 10-15$ km s*-1*. In this case, accurate investigation of MFT supersonic motion with bow shocks can be done by including dependence of the aerodynamic drag on the Mach number. We plan to do this in future works.

The dynamics of the MFT consists in two stages. At the first stage, toroidal magnetic field is generated and MFT form due to magnetic buoyancy instability over time $t_{{\rm gen}}$. At the second stage, the MFT rise from the disc to its suface over the time $t_{{\rm surf}}$. In a process of further upward motion the MFT carry away some mass and magnetic flux to disc atmosphere. Our calculations show that $t_{{\rm surf}}\ll t_{{\rm gen}}$. Therefore, the process of mass and magnetic flux transport from the disk to its atmosphere due to buoyancy is periodic with typical period $t_{{\rm gen}}\sim 0.5-10$ yr. The vertical mass transport rate due to buoyancy $\dot{M}_{{\rm b}}\sim 10^{-12}-10^{-7}\,M_{\sun}\,{\rm yr}^{-1}$. Our calculations show that approximately $20$ % of disc mass flux can come out from the disc via buoyancy. The rising MFT can be the seed for the formation of jets and outflows from accretion discs. The MFT carry magnetic fluxes $\Phi_{{\rm b}}\sim 10^{19}-10^{24}$ Mx, so that magnetic flux of $\Phi_{{\rm b}}\approx 10^{30}$ Mx, that comprises $20$ % of total disc magnetic flux, can be carried away from the disc by rising MFT during period of $1$ Myr. We assume that magnetic buoyancy is the mechanism responsible for the efficient magnetic flux escape from the accretion discs of young stars. Khaibrakhmanov & Dudorov (2017) have shown that formation and rise of the MFT with initial radius $a_{0}=0.1\,H$ stabilize the strength of the toroidal magnetic field at the level of the poloidal magnetic field strength.

The buoyancy decreases and MFT acquires steady speed above the disc. This is explained by the fact that the external density $\rho_{{\rm e}}$ and the density of MFT decrease with $z$-coordinate. Correspondingly, the buoyant and drag forces also reduce, acceleration of MFT approaches to zero, MFT moves by inertia and significantly expands above the disc. We interpret this process as a formation of an expanding magnetized ‘corona’ above the disc. Expansion of the MFT can be reduced by the internal magnetic tension. We will consider this possibility in future works.

Second important case investigated by us concerns the effect of the magnetic field of the disc on the dynamics of MFT. In this case, MFT rise from the disc and start to oscillate at height $z_{{\rm osc}}\approx 2-2.5H\,$ above the midplane of the disc. The magnetic oscillations take place near the point where strengths of internal, $B$, and external, $B_{{\rm e}}$, magnetic fields are nearly equal. Above this point, the buoyant force is negative since $B<B_{{\rm e}}$ and the MFT has greater density than the external gas. The oscillation periods $P_{{\rm osc}}$ increase with the distance from the star, $P_{{\rm osc}}=(1-10)$ days at $r=0.027$ au and $P_{{\rm osc}}=(1-3)$ months at $r=0.6$ au. The MFT with stronger initial magnetic field (smaller plasma beta) rise to higher altitudes and experience oscillations with smaller periods. The magnetic oscillations decay with time, and occur in the disc over time of the toroidal magnetic field generation. The oscillations of the MFT found in our simulations indicate that the mass carried by the MFT rsing from the disc can come back to the disc. Accumulation of the magnetic flux near the surface of the discs can further lead to the burst release of the magnetic energy as it was suggested by Shibata et al. (1990). Magnetic energy of the MFT of $\sim 10^{33}-10^{37}$ erg can be released due to such bursts. Magnetic field $B_{{\rm e}}$ can decrease with height, and the MFT probably can rise to higher altitudes before reaching the point of zero buoyancy in this case.

Periodic formation, rise and thermal or magnetic oscillations of the MFT can be interpreted as periodical changes in the accretion disc structure. This process is similar to periodic ejection of buoyant magnetic fields from the disc that has been found in frame of numerical MHD simulations by (e.g. Turner et al., 2010; Takasao et al., 2018). At the distance $r=0.6$ au, temperature is $\sim 1000$ K, and the refractory dust grains are present in this region. Therefore, buoyant MFT contain dust, and the oscillations can produce IR-variability of the disc.

IR-variability is the common feature of YSO (see Flaherty et al. (2013, 2016) and references therein). Flaherty et al. (2016) have found that the low-mass YSO in the Chameleon I star forming region exhibit variability over time scales of months (20-200 days). Classical T Tauri stars also exhibit variability in optical wavelengths. For example, Rigon et al. (2017) reported the detection of the optical variability with periods of $20-60$ d of classical T Tauri stars in Taurus-Auriga region.

We simulated MFT dynamics in the discs of two classical T Tauri stars J11092266-7634320 and J11100369-7633291 (439 and 530) from the sample presented in Flaherty et al. (2016). The magnetic oscillation periods found in our simulations are in good agreement with the observed periods for considered stars. Our simulations have shown that variability period of star 439 (32 days) is consistent with the oscillations of the MFT with initial $\beta_{0}=1$ and $a_{0}=0.1\,H$. Observed period for star 530 (35 days) can be explained by the oscillations of the MFT with $a_{0}=0.1\,H$ and $\beta_{0}=5$.

Our estimations have shown that the MFT are optically thick with respect to IR-radiation. The variability can be caused by MFT temperature fluctuations, as well as by periodical screening the outer disc regions from stellar radiation by the MFT. It should be noted that Takasao et al. (2018) found formation of the MFT in 3D MHD simulations of the inner regions of the accretion discs and also argued that periodically rising MFT may contribute to the variability of discs radiation.

The MFT can form in the region $r_{{\rm i}}<r<r_{{\rm o}}$ au ($r_{{\rm i}}=r(T=1500\,K)$, $r_{{\rm o}}$ is the outer boundary of the zone of the thermal ionization), where magnetic field is efficiently amplified and temperature is lower than 1500 K, so that the MFT contain the refractory dust particles. For stars 439 and 530 radii $r_{{\rm i}}=0.2$ au, $0.07$ au and $r_{{\rm o}}=0.42$ au, $0.1$ au, respectively. Therefore, the oscillations occur in a range of radii, which is consistent with conclusions of Flaherty et al. (2016). Magnetic oscillations decay with time. For example, magnitude of temperature oscillations reduces by $100$ K over $1.3$ months at $r_{{\rm i}}=0.2$ and over $4.7$ months at $r_{{\rm o}}=0.42$ in the disc of star J11092266-7634320. Thus, significant temperature variations of oscillating MFT, $\Delta T\sim 100$ K, can be detected within time of $200$ d corresponding to observation time in Flaherty et al. (2016).

Tambovtseva & Grinin (2008) suggested that inhomogeneities in the disc centrifugal winds containing both gas and dust can cause variations in circumstellar extinction observed in T Tauri stars. Our simulations show that magnetic buoyancy can transport mass from the accretion disc to its atmosphere. We propose that such this process can be responsible for the variations in circumstellar extinction observed in T Tauri stars. Similar hypothesis was proposed by Miyake et al. (2016), who simulated the time evolution of the dust grain distribution in the vertical direction inside the minimum mass solar nebula.

As Schramkowski (1996) have shown, longitudinal perturbations of the magnetic flux tubes can lead to the formation of the arcs above the disc. We assume that such orbiting arcs can cast shadows on the outer regions of the accretion discs. This problem needs further investigation.

Toroidal magnetic field in the disc is generated from the poloidal one over the time scale of rotation period. Tension of the poloidal magnetic field lines can slow down the disc rotation, and, as a consequence, hinder the generation of the toroidal magnetic field in the considered region. This effect takes place if the energy of the poloidal magnetic field is comparable with the rotational energy. This condition for the keplerian disc can be written as

[TABLE]

i.e. the magnetic field with $\beta\sim 10^{-4}$ will influence significantly the disc rotation. Kinematic approximation $\beta>1$ is used in our model of the accretion disc, and the poloidal magnetic field cannot be strong enough to stop generation of the toroidal magnetic field and formation of the MFT.

We did not consider the effect of the azimuthal velocity shear on the MFT dynamics. Shibata et al. (1990) have found that rise velocity of the MFT in the presence of velocity shear is smaller than in the case of no shear. MHD simulations of MFT rising from the upper convection zone of the Sun to the solar atmosphere have shown that the MFT dynamics is sensitive to the twist of the MFT (see Fan et al., 1998; Magara, 2001; Martínez-Sykora et al., 2015). Strongly twisted MFT retain their coherent structures during the rise, while the MFT without twist splits into a vortex pair and lose significant amount of their magnetic flux. Up to date, there are no detailed simulations of the MFT formation in the accretion discs. It is hard to make conclusions about degree of their twist. In our simulations we implicitly assumed that the MFT retain their coherence and do not lose the magnetic flux. Probably, loss of the magnetic flux would lead to decrease of the MFT velocity.

The magnetic pressure acts only in the direction perpendicular to the magnetic lines. Uniform magnetic field of the disc $B_{{\rm e}}$ considered in this work will lead to flattening of the MFT in the direction of motion. Generally speaking, the toroidal magnetic field is generated not only in the region of MFT formation, but in the entire volume of the disc. Therefore, toroidal magnetic field of the disc will influence MFT dynamics together with the poloidal magnetic field. Assuming that intensities of the toroidal and poloidal magnetic fields are comparable, we adopted in this work that the magnetic pressure is isotropic, and the cross-section of the MFT remains nearly round. Effects of non-uniform MFT expansion can be investigated in two-dimensional model of the MFT dynamics, that we aim to elaborate in future. We plan to develop more detailed model of the vertical structure of the disc in the future. Interesting task is investigation of the dynamics of the magnetic rings in the accretion discs of young stars, i.e. study the evolution of major radius of the toroidal MFT.

Acknowledgements

The work of AED is supported by the Russian Foundation for Basic Research (project 18-02-01067). The work of SAKh and AMS is supported by the Ministry of Science and High Education (the basic part of the State assignment, RK No. AAAA-A17-117030310283-7) and by the Act 211 Government of the Russian Federation, contract No. 02.A03.21.0006. The authors thank prof. Dmitry Bisikalo, Dr. Sergey Parfenov, Dr. Anna Evgrafova and Dr. Vitaly Akimkin for useful comments. We are also grateful to Lyudmila Lapina for checking the English language in the paper. We thank anonymous referee for a detailed review and useful comments.