Influence of Ohmic and ambipolar heating on thermal structure of accretion discs

TL;DR

This paper models the impact of Ohmic and ambipolar heating on the thermal structure of accretion discs around young stars, revealing significant temperature increases near the 'dead' zone boundaries due to dissipative MHD effects.

Contribution

It introduces a comprehensive MHD model including Ohmic and ambipolar effects, analyzing their influence on accretion disc temperature structure around T Tauri stars.

Findings

Temperature increases by ~1000 K near the inner 'dead' zone boundary.

Temperature increases by ~100 K near the outer 'dead' zone boundary.

Dissipative effects are significant in low ionization regions.

Abstract

We investigate dynamics of accretion discs of young stars with fossil large-scale magnetic field. Our magneto-gas-dynamic (MHD) model of the accretion discs includes equations of Shakura and Sunyaev, induction equation, equations of thermal and collisional ionization. Induction equation takes into account Ohmic and magnetic ambipolar diffusion, magnetic buoyancy. We also consider the influence of Ohmic and ambipolar heating on thermal structure of the accretion discs. We analyse the influence of considered dissipative MHD effects on the temperature of the accretion discs around classical T Tauri star. The simulations show that Ohmic and ambipolar heating operate near the borders of the region with low ionization fraction (`dead' zone). Temperature grows by $\sim 1000$ K near the inner boundary of the `dead' zone, $r\approx (0.5-1)$ au, and by $\sim 100$ K near its outer boundary, $r\approx (30-50)$ au.

Figures

Equations

$$\frac{\partial \rho }{\partial t}+\text{Div}\left(\rho \text{bvec}v\right)$$

$$\rho \frac{\partial \text{bvec}v}{\partial t}+\left(\text{bvec}v\nabla \right)\text{bvec}v$$

$$\rho T\frac{\partial s}{\partial t}$$

$$\frac{\partial \text{bvec}B}{\partial t}$$

$$\Gamma ={\Gamma }_{\mathrm{cr}}+{\Gamma }_{\mathrm{o}}+{\Gamma }_{\mathrm{a}}.$$

$${\Gamma }_{\mathrm{o}}=\frac{{\nu }_{\mathrm{m}}}{4\pi }{\left(\text{Rot}\text{bvec}B\right)}^2.$$

$${\Gamma }_{\mathrm{a}}=\frac{1}{2}{v}_{\mathrm{ad}}^2{R}_{\mathrm{in}},$$

$$\text{bvec}{v}_{\mathrm{ad}}=\frac{[\text{Rot}\text{bvec}B,\text{bvec}B]}{4\pi R_{\mathrm{in}}}$$

$$R_{\mathrm{in}}={\mu }_{\mathrm{in}}{n}_{\mathrm{i}}{n}_{\mathrm{n}}⟨\sigma v{⟩}_{\mathrm{in}},$$

$$\dot{M}$$

$$\dot{M}{\Omega }_{\mathrm{k}}f$$

$$v_{\phi }$$

$$H$$

$${\sigma }_{\mathrm{sb}}{T}_{\mathrm{eff}}^4$$

$$T^4$$

$$B_r$$

$$B_{\phi }$$

$$B_z$$

$$\eta ={\nu }_{\mathrm{m}}+{\eta }_{\mathrm{a}}+{\eta }_{\mathrm{b}}$$

$${\nu }_{\mathrm{m}}=\frac{c^2}{4\pi \sigma },$$

$${\eta }_{\mathrm{a}}=\frac{B_z^2}{4\pi R_{\mathrm{in}}}.$$

$${\eta }_{\mathrm{b}}=H|v_{\mathrm{b}}|.$$

$$\left(1-x_{\mathrm{s}}\right)\xi ={\alpha }_{\mathrm{r}}{x}_{\mathrm{s}}^2{n}_{\mathrm{n}}+{\alpha }_{\mathrm{g}}{x}_{\mathrm{s}}{n}_{\mathrm{n}},$$

$$x=x_{\mathrm{s}}+\sum _j{\nu }_j{x}_j^{\mathrm{T}},$$

Full text

We also consider the influence of Ohmic and ambipolar heating on thermal structure of the accretion discs. We analyse the influence of considered dissipative MHD effects on the temperature of the accretion discs around classical T Tauri star. The simulations show that Ohmic and ambipolar heating operate near the borders of the region with low ionization fraction (‘dead’ zone). Temperature grows by $\sim 1000$ K near the inner boundary of the ‘dead’ zone, $r\approx(0.5-1)$ au, and by $\sim 100$ K near its outer boundary, $r\approx(30-50)$ au.

Introduction.

Observations in optical, infrared (IR), submillimeter and radio ranges have shown that young stars are surrounded by gas-dust discs [1]. The discs are geometrically thin and optically thick. Typical masses of the discs are $0.001-0.1M_{\odot}$, radii are 100-1000 au111$1\,\mathrm{AU}=1.5\times 10^{13}$ cm. Matter accretes from the disc onto the star with mass rates $\dot{M}=(10^{-9}-10^{-6})\,M_{\odot}\,\mathrm{yr}^{-1}$ over time $10^{6}-10^{7}$ yr.

It has been found that thermal emission of dust grains in the accretion discs of young stars (ADYS) is polarized [2, 3, 4]. The polarization can be caused by paramagnetic dust grain alignment by the magnetic field in the disc. Therefore, the polarimetric studies point to existence of large scale magnetic field in the accretion discs. Measurements of the magnetic field strength in the ADYS are still challenging.

Young stars with accretion discs form as a result of gravitational collapse of magnetized cores of molecular clouds. Numerical simulations have shown that initial magnetic flux of the cloud cores partially conserves during the star formation (see reviews [5, 6]). Therefore it is natural to assume that the large-scale magnetic field of the accretion disc has fossil nature, i.e. it is the remnant of the magnetic field of parent molecular cloud core [7, 6]. Possibility of dynamo generation of the magnetic field in the accretion discs is also considered in a number of works (see, for example, [8], [9], [10]).

Accretion and differential rotation lead to generation and amplification of the magnetic field in the regions with high ionization fraction. Cosmic rays, X-rays and radioactive decay are the main sources of ionization in the ADYS. Attenuation of the cosmic and X-rays in the surface layers of the discs leads to formation of layered ioinization structure of the discs [11]. Region of low ionization fraction, $x\leq 10^{-12}$, so-called ‘dead’ zone, lies typically between 0.5-1 au and 10-30 au near the midplane of the accretion discs of T Tauri stars [12, 13, 14]. Ohmic and ambipolar diffusion hinder generation of the magnetic field inside the ‘dead’ zones [14]. Dudorov and Khaibrakhmanov  [14] have shown that the magnetic field geometry varies over the disc. The magnetic field is frozen into gas and has quasi-azimuthal geometry in the inner region, where thermal ionization operates. The magnetic field preserves its initial poloidal geometry inside the ‘dead’ zones. The magnetic field can have quasi-radial or quasi-azimuthal geometry in the outer regions of the ADYS depending on the intensity of ionization sources and dust grain parameters. We have not considered contribution of the dissipative MHD effects to thermal processes yet. Lizano et al. [15] calculated thermal structure of the accretion discs taking into account turbulent and Ohmic heating, and found that both processes can have comparable effect on the temperature of the discs of T Tauri stars. But they assumed that resistivity is uniform in the disc that does not allow to correctly study role of the Ohmic heating in the discs with ‘dead’ zones.

In this work, our aim is to investigate effect of Ohmic and ambipolar heating on thermal structure of ADYS on the basis of our MHD model of the accretion disc [14, 16].

The paper is organized as follows. In section 1, we describe the problem statement and main equations of the model. Results of the simulations are presented in section 2. We summarize and discuss the results in section 3.

1 Problem statement and main equations.

We consider stationary geometrically thin and optically thick accretion disc around solar mass classical T Tauri star (see schematic picture in Figure 1). We use cylindrical coordinates $(r,\,\varphi,\,z)$. Disc rotation axis is directed along $z$-axis. The disc is considered to be in hydrostatic equilibrium in the $z$-direction, so that the velocity components are $\bvec{v}=(v_{r},\,v_{\varphi},\,0)$. Magnetic field components: $\bvec{B}=(B_{r},\,B_{\varphi},\,B_{z})$. Inner edge of the disc $r_{{\rm in}}$ is determined by the boundary of stellar magnetosphere. Outer edge $r_{{\rm out}}$ is the contact boundary with the interstellar medium.

We investigate the dynamics of the accretion discs using following system of MHD equations taking into account viscous stresses, thermal conductivity, magnetic ambipolar diffusion, Ohmic diffusion, and magnetic buoyancy:

[TABLE]

where $\Phi$ is the gravitational potential of the star, $\sigma^{{}^{\prime}}$ is the viscous stress tensor (see [17]), $s$ is the entropy, $\kappa_{{\rm r}}$ is the radiative heat conductivity (see [18]), $\Gamma$ describes additional heating sources, $\bvec{v}_{{\rm ad}}$ is the velocity of ambipolar diffusion, $\nu_{{\rm m}}$ is the magnetic viscosity, $\bvec{v}_{{\rm b}}$ is the buoyant velocity (see [19]), $\bvec{B}_{{\rm t}}=(0,\,B_{\varphi},\,0)$ is the toroidal component of the magnetic field, $P=\rho R_{{\rm g}}T/\mu$, $s=c_{{\rm v}}{\rm ln}\left(P/\rho^{\gamma}\right)$, $R_{{\rm g}}$ is the universal gas constant, $\mu=2.3$ is the mean molecular weight of the gas, $c_{{\rm v}}$ is the specific heat at constant volume, $\gamma$ is the adiabatic index. Other quantities are used in their usual physical notation.

The first term on right-hand side of Equation (3) describes heating of the gas by viscous friction, $\Gamma_{{\rm v}}$. In addition to the viscous heating, we consider heating by cosmic rays and radioactive decay, Ohmic (Joule) and ambipolar heating,

[TABLE]

Rate of heating by cosmic rays and radioactive decay $\Gamma_{{\rm cr}}$ is calculated according to [20]. Rate of Ohmic heating

[TABLE]

Rate of ambipolar heating [21]

[TABLE]

where ambipolar velocity in stationary approximation

[TABLE]

and coefficient of friction between ions and neutrals

[TABLE]

$\mu_{\rm{in}}$ is the reduced mass of ions and neutrals, $n_{\rm{i}}$ is the number density of ions, $n_{\rm{n}}$ is the number density of neutrals, $\langle\sigma v\rangle_{\rm{in}}=2\times 10^{-9}\,\mbox{cm}^{3}\,\mbox{s}^{-1}$ is the ‘slowing-down’ coefficient [22]. We consider plasma consisting of hydrogen, helium and metals with standard cosmic abundance. Mean mass of neutrals is $2.3\,m_{{\rm p}}$, mean mass of metal ions $m_{{\rm i}}=30\,m_{{\rm p}}$, where $m_{{\rm p}}$ is the proton mass.

We solve equations (1-4) in frame of Shakura and Sunyaev [23] approximation (see [14] and [16] for details). It is considered that angular momentum is transported radially by turbulent stresses. Molecular viscosity in the stress tensor $\sigma^{{}^{\prime}}$ is replaced by turbulent viscosity $\nu_{{\rm t}}=\alpha v_{{\rm s}}H$, where $\alpha<1$ is the non-dimensional parameter characterizing turbulence efficiency, $v_{{\rm s}}$ is the sound speed, $H$ is the accretion disc’s scale height. We assume that the disc is in centrifugal and hydrostatic balance, and that the heat released inside the disc is transported vertically by radiation. Radial derivatives are neglected in comparison with vertical derivatives, and the latter are replaced by finite differences. Then Equations (1-3) transform to

[TABLE]

where $\dot{M}$ is the accretion rate, $\Sigma=2\rho H$ is the surface density of the disc, $\Omega_{{\rm k}}=\sqrt{GM/r^{3}}$ is the Keplerian angular speed, $G$ is the gravitational constant, $f=1-\sqrt{r_{{\rm in}}/r}$, $T_{{\rm eff}}$ is the effective temperature of the disc, $T$ is the temperature inside the disc, $\tau=\kappa\Sigma$ is the optical thickness, and $\kappa$ is the opacity. Solution of the induction Equation (4) in approximations of our model:

[TABLE]

where $B_{z0}$ and $\Sigma_{0}$ are the magnetic field strength and surface density at the outer edge of the disc, $x$ is the ionization fraction, $R_{{\rm m}}$ is the magnetic Reynolds number,

[TABLE]

is the effective diffusivity including effects of Ohmic diffusion (the first term), ambipolar diffusion (the second term), and magnetic buoyancy (the third term). Magnetic viscosity

[TABLE]

where $\sigma$ is the electric conductivity (see book [24]). Ambipolar diffusivity

[TABLE]

Diffusivity due to magnetic buoyancy

[TABLE]

Velocity of magnetic buoyancy $v_{{\rm b}}$ is determined from the balance between buoyant and drag forces according to [25]. Aerodynamic and turbulent drags are considered.

Solutions (16) and (17) determine radial and azimuthal components of the magnetic field from the balance between advection and diffusion in corresponding directions. Solution (18) for $R_{{\rm m}}\gg 1$ describes the magnetic field that is frozen into gas in regions with high conductivity. Solution (18) for $R_{{\rm m}}\leq 1$ determines $B_{z}$ from the equality of advection and ambipolar diffusion velocities in the regions with small conductivity.

Ionization fraction is calculated following [26], $x=(n_{{\rm e}}+n_{{\rm i}})/(n_{{\rm e}}+n_{{\rm i}}+n_{{\rm n}})$, where $n_{{\rm e}}\approx n_{{\rm i}}$ is the electrons concentration. Collisional ionization fraction $x_{{\rm s}}$ is calculated from the equation

[TABLE]

where $\xi$ is the total ionization rate by cosmic rays, X-rays and decay of radionuclides, $\alpha_{{\rm r}}$ is the rate of radiative recombinations, $\alpha_{{\rm g}}$ is the rate of recombinations onto the dust grains. We consider spherical grains with initial mean radius $0.1\,\mu m$ and standard composition. Evaporation of grains is taken into account.

In the regions with high temperature $T>400$ K, we calculate thermal ionization $x_{j}^{{\rm T}}$ of hydrogen ($j={\rm H}$) and alkali metal atoms ($j={\rm m}$) using the Saha equation. ‘Mean’ metal with ionization potential of $5.76$ eV and abundance logarithm of $5.97$ is considered for simplicity. Total ionization fraction

[TABLE]

were $\nu_{j}$ is the abundance of element $j$ with respect to hydrogen.

System of algebraic equations (10-24) is non-linear. We solve it numerically using the iterative method together with the bisection method. Solution without Ohmic and ambipolar heating found by [14] is used as an initial approximation.

2 Results.

In this work we present illustrative results for the following parameters: $M=1\,M_{\odot}$, stellar radius $R=2\,R_{\odot}$, magnetic field strength at stellar surface $B_{\star}=2$ kG, $\dot{M}=10^{-7}\,M_{\odot}\,{\rm yr}^{-1}$, $\alpha=0.01$, cosmic ray ionization rate $\xi_{0}=10^{-17}\,{\rm s}^{-1}$, cosmic ray attenuation length $100\,{\rm g}\,{\rm cm}^{-2}$. For simplicity, we use opacity $\kappa=10\,{\rm cm}^{2}\,{\rm g}^{-1}$ to approximately describe thermal structure of the disc at low temperatures. We set value of $B_{z0}$ at the outer boundary of the disc equal to the magnetic field strength of parent protostellar cloud (see details in [16]).

2.1 Ionization fraction and magnetic field.

The Ohmic and ambipolar heating are determined by ionization fraction and magnetic field strength, as Equations (6) and (7) show. In Figure 2, we plot radial profiles of the ionization fraction (panel (a)) and magnetic field components (panel (b)). It is considered that $B_{\varphi}$ strength is limited by the equipartition value $B_{{\rm eq}}$. The equipartition value corresponds to plasma beta $\beta=1$.

Figure 2(a) shows that the radial profile of the ionization fraction is non-monotonic. The ionization fraction is high in the inner region of the disc, $r<1$ au, due to thermal ionization. The ionization fraction increases with distance due to decrease of gas density and more efficient ionization by cosmic rays at $r>1$ au. There is region of low ionization fraction, $x\leq 10^{-12}$, at $0.8$ au $<r<$ $21$ au. This region is called ‘dead’ zone [11].

Figure 2(b) shows that there is three regions with different magnetic field geometry in the disc. In the region of thermal ionization, $r<1$ au, the magnetic field is frozen into gas, and efficient generation and amplification of all three components take place. Azimuthal component has fastest generation rate, and $B_{\varphi}>(B_{r},\,B_{z})$ in this region. Inside the ‘dead’ zone, the magnetic field is quasi poloidal due to efficient Ohmic diffusion, $B_{z}\gg(B_{r},\,B_{\varphi})$. In the outer region, $r>20$ au, ambipolar diffusion hinders generation of $B_{r}$, and the magnetic field has quasi-azimuthal geometry, $B_{\varphi}\geq B_{z}$.

2.2 Effect of Ohmic and ambipolar heating on thermal structure of the disc.

As it was mentioned in section 1, four heating sources are considered in our model of the disc: viscous heating ($\Gamma_{{\rm v}}$), heating by cosmic rays and radioactive decay ($\Gamma_{{\rm cr}}$), Ohmic ($\Gamma_{{\rm o}}$) and ambipolar heating ($\Gamma_{{\rm a}}$). We carried out two simulations: one with viscous and heating by cosmic rays and radioactive decay ($\Gamma_{{\rm v}}+\Gamma_{{\rm cr}}$), and another with additional effect of Ohmic and ambipolar heating ($\Gamma_{{\rm v}}+\Gamma_{{\rm cr}}+\Gamma_{{\rm o}}+\Gamma_{{\rm a}}$). In Figure 3, we plot radial profiles of disc temperature (panel (a)) and surface density (panel (b)) for this two cases. Figure 3(a) shows that Ohmic and ambipolar heating leads to increase of the temperature near the edges of the ‘dead’ zone. Temperature increases by $1000-3000$ K near the inner edge, $r\approx(0.5-1)$ au, and by $100$ K near the outer edge, $r\approx(30-50)$ au. Figure 3(b) shows that temperature growth leads to decrease in surface density. It happens because the accretion speed $v_{r}$ is higher in regions with high temperature, and therefore matter is accreted more efficiently from these regions.

Temperature growth and density decrease lead to increase of the ionization fraction in the regions of Ohmic and ambipolar heating. As a consequence, size of the ‘dead’ zone reduces. The ‘dead’ zone is situated between $1.3$ au and $20$ au, when Ohmic and ambipolar heating operate.

3 Conclusions.

We investigated thermal structure of the magnetized accretion discs of classical T Tauri stars. Ohmic and ambipolar heating are included in the MHD model of the accretion discs.

The simulations have shown that the dissipative MHD effects operate near the boundary of the ‘dead’ zone. The ‘dead’ zone is situated between $0.8$ au and $21$ au in the simulation without Ohmic and ambipolar heating. Ambipolar heating is determined by $B_{\varphi}$ strength in the regions of maximal ambipolar velocity. In the case of maximum $B_{\varphi}=B_{{\rm eq}}$, temperature increases by value $\leq 3000$ K near the inner boundary of the ‘dead’ zone, $r\approx(0.5-1)$ au, and by 100 K near its outer boundary $r\approx(30-50)$ au. If the magnetic buoyancy operates, then $B_{\varphi}\approx B_{z}$, and ambipolar heating affects only inner boundary of the ‘dead’ zone.

Radial extent of the ‘dead zone’ reduces due to the action of the Ohmic and ambipolar heating. We suppose that size of the ‘dead’ zone in the vertical direction would also reduce. Amplification of the fossil magnetic field is possible outside the ‘dead’ zones, as well as dynamo generation from the seed fossil field. We conclude that the dissipative MHD effects support amplification and generation of the magnetic field in the accretion discs of young stars.

Temperature growth due to magnetic field dissipation can lead to evaporation of dust grains near the borders of ‘dead’ zones in the accretion discs of young stars. These regions with evaporated dust grains can probably be observed in the maps of discs’ thermal radiation.

\Thanks

We thank anonymous referee for useful comments that helped to improve our paper. The work is supported by Russian Foundation for Basic Research (project 18-02-01067) and by the Ministry of Science and High Education (the basic part of the State assignment, RK no. AAAA-A17- 117030310283-7).