Influence of Ohmic and ambipolar heating on thermal structure of accretion discs
Sergey A. Khaibrakhmanov Alexander E. Dudorov



Sparse Nonnegative CANDECOMP/PARAFAC Decomposition in Block Coordinate Descent Framework: A Comparison Study
    Deqing Wang, Student Member, IEEE,
        Fengyu Cong, Senior Member, IEEE,
        and Tapani Ristaniemi, Senior Member, IEEE
This work was supported in part by the National Natural Science Foundation of China (Grant No. 81471742), in part by the Fundamental Research Funds for the Central Universities [DUT16JJ(G)03] in Dalian University of Technology in China, and in part by the scholarship from China Scholarship Council (No. 201600090043). (Corresponding author: Fengyu Cong)
D. Wang and F. Cong are with the School of Biomedical Engineering, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, Dalian 116024, China, and also with the Faculty of Information Technology, University of Jyväskylä, Jyväskylä 40100, Finland (e-mail: deqing.wang@foxmail.com; cong@dlut.edu.cn).
T. Ristaniemi is with the Faculty of Information Technology, University of Jyväskylä, Jyväskylä 40100, Finland (e-mail: tapani.e.ristaniemi@jyu.fi).
Manuscript received Month Day, Year; revised Month Day, Year.
    December 30, 2023
=============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



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 ∼ 1000 K near the inner boundary of the `dead' zone, r≈ (0.5-1) au, and by ∼ 100 K near its outer boundary, r≈ (30-50) au.







§ INTRODUCTION.

Observations in optical, infrared (IR), submillimeter and radio ranges have shown that young stars are surrounded by gas-dust discs <cit.>. The discs are geometrically thin and optically thick. Typical masses of the discs are 0.001-0.1M_⊙, radii are 100-1000 au[1 AU = 1.5× 10^13 cm]. Matter accretes from the disc onto the star with mass rates Ṁ=(10^-9-10^-6) M_⊙ 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 <cit.>. 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 <cit.>). 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 <cit.>. Possibility of dynamo generation of the magnetic field in the accretion discs is also considered in a number of works (see, for example, <cit.>, <cit.>, <cit.>).

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 <cit.>. Region of low ionization fraction, x≤ 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 <cit.>. Ohmic and ambipolar diffusion hinder generation of the magnetic field inside the `dead' zones <cit.>. Dudorov and Khaibrakhmanov  <cit.> 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. <cit.> 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 <cit.>.

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



§ 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 <ref>). We use cylindrical coordinates (r, φ, 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 v=(v_r, v_φ, 0). Magnetic field components: B=(B_r, B_φ, B_z). Inner edge of the disc r_ in is determined by the boundary of stellar magnetosphere. Outer edge r_ 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:

    ∂ρ/∂ t + (ρv)    =    0,
    ρ∂v/∂ t + (v∇)v   =   - P - ρΦ +1/4π[B, B] + σ^',
    ρ T∂ s/∂ t   =   σ_ik^'∂ v_i/∂ x_k + (κ_ r T) + Γ,
    ∂B/∂ t   =   [v + v_ ad, B] + ν_ m∇^2B + [v_ b, B_ t],

where Φ is the gravitational potential of the star, σ^' is the viscous stress tensor (see <cit.>), s is the entropy, κ_ r is the radiative heat conductivity (see <cit.>), Γ describes additional heating sources, v_ ad is the velocity of ambipolar diffusion, ν_ m is the magnetic viscosity, v_ b is the buoyant velocity (see <cit.>), B_ t=(0, B_φ, 0) is the toroidal component of the magnetic field,  P=ρ R_ gT/μ, s=c_ v ln(P/ρ^γ), R_ g is the universal gas constant, μ=2.3 is the mean molecular weight of the gas, c_ v is the specific heat at constant volume, γ is the adiabatic index. Other quantities are used in their usual physical notation.



The first term on right-hand side of Equation (<ref>) describes heating of the gas by viscous friction, Γ_ v. In addition to the viscous heating, we consider heating by cosmic rays and radioactive decay, Ohmic (Joule) and ambipolar heating,

    Γ = Γ_ cr + Γ_ o + Γ_ a.

Rate of heating by cosmic rays and radioactive decay Γ_ cr is calculated according to <cit.>. Rate of Ohmic heating

    Γ_ o = ν_ m/4π(B)^2.

Rate of ambipolar heating <cit.>

    Γ_ a = 1/2v_ ad^2R_ in,

where ambipolar velocity in stationary approximation

    v_ ad = [B, B]/4π R_ in
 
and coefficient of friction between ions and neutrals

    R_ in = μ_inn_in_n⟨σ v⟩_in,

μ_in is the reduced mass of ions and neutrals, n_i is the number density of ions, n_n is the number density of neutrals, ⟨σ v⟩_in=2× 10^-9 ^3 ^-1 is the `slowing-down' coefficient <cit.>. We consider plasma consisting of hydrogen, helium and metals with standard cosmic abundance. Mean mass of neutrals is 2.3 m_ p, mean mass of metal ions m_ i=30 m_ p, where m_ p is the proton mass.

We solve equations (<ref>-<ref>) in frame of Shakura and Sunyaev <cit.> approximation (see <cit.> and <cit.> for details). It is considered that angular momentum is transported radially by turbulent stresses. Molecular viscosity in the stress tensor σ^' is replaced by turbulent viscosity ν_ t = α v_ s H, where α<1 is the non-dimensional parameter characterizing turbulence efficiency, v_ 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 (<ref>-<ref>) transform to

    Ṁ   =    - 2 π r v_r Σ, 
    ṀΩ_ k f    =    2 πα v_ s^2 Σ, 
    
      v_φ   =   √(GM/r)(1 + z^2/r^2)^-3/4,
    
      H    =   v_ s/Ω_ k,
    σ_ sb T_ eff^4    =   3/8πṀΩ_ k^2 f + 2HΓ_ cr + 
          ν_ m/4πB_r^2 + B_φ^2/H + (B_zB_r)^2 + (B_zB_φ)^2 + (B_r^2 + B_φ^2)^2/32π^2 R_ inH,
    
      T^4    =    T_ eff^4 (1/2 + 3/4τ) ,

where Ṁ is the accretion rate, Σ=2ρ H is the surface density of the disc, Ω_ k=√(GM/r^3) is the Keplerian angular speed, G is the gravitational constant, f = 1 - √(r_ in / r),  T_ eff is the effective temperature of the disc, T is the temperature inside the disc, τ=κΣ is the optical thickness, and κ is the opacity. Solution of the induction Equation (<ref>) in approximations of our model:
  
    B_r    =    -v_r H/ηB_z, 
    
      B_φ   =    -3/2(H/r)^2v_φH/ηB_z - 1/2(H/r)v_φH/ηB_r, 
    
     B_z    =   
      B_z0Σ/Σ_0    R_ m≫ 1,
    √(4π xρ^2 r|v_r|)    R_ m≤ 1,

where B_z0 and Σ_0 are the magnetic field strength and surface density at the outer edge of the disc, x is the ionization fraction, R_ m is the magnetic Reynolds number,

    η = ν_ m + η_ a + η_ b

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 

    ν_ m = c^2/4πσ,

where σ is the electric conductivity (see book <cit.>). Ambipolar diffusivity

    η_ a =  B_z^2/4π R_ in.

Diffusivity due to magnetic buoyancy

    η_ b = H|v_ b|.

Velocity of magnetic buoyancy v_ b is determined from the balance
between buoyant and drag forces according to <cit.>. Aerodynamic and turbulent drags are considered.

Solutions (<ref>) and (<ref>) determine radial and azimuthal components of the magnetic field from the balance between advection and diffusion in corresponding directions. Solution (<ref>) for R_ m≫ 1 describes the magnetic field that is frozen into gas in regions with high conductivity. Solution (<ref>) for R_ m≤ 1 determines B_z from the equality of advection and ambipolar diffusion velocities in the regions with small conductivity.

Ionization fraction is calculated following <cit.>, x=(n_ e + n_ i) / (n_ e + n_ i + n_ n), where n_ e≈ n_ i is the electrons concentration. Collisional ionization fraction x_ s is calculated from the equation

    (1 - x_ s)ξ = α_ rx_ s^2n_n + α_ gx_ sn_n,

where ξ is the total ionization rate by cosmic rays, X-rays and decay of radionuclides, α_ r is the rate of radiative recombinations, α_ g is the rate of recombinations onto the dust grains. We consider spherical grains with initial mean radius 0.1 μ 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^ T of hydrogen (j= H) and alkali metal atoms (j= 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

    x = x_ s + ∑_jν_jx_j^ T,

were ν_j is the abundance of element j with respect to hydrogen.

System of algebraic equations (<ref>-<ref>) is non-linear. We solve it numerically using the iterative method together with the bisection method. Solution without Ohmic and ambipolar heating found by <cit.> is used as an initial approximation.



§ RESULTS.


In this work we present illustrative results for the following parameters: M=1 M_⊙, stellar radius R=2 R_⊙, magnetic field strength at stellar surface B_⋆=2 kG, Ṁ=10^-7 M_⊙  yr^-1, α=0.01, 
cosmic ray ionization rate ξ_0=10^-17  s^-1, cosmic ray attenuation length 100  g  cm^-2. For simplicity, we use opacity κ=10  cm^2  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 <cit.>).



 §.§ Ionization fraction and magnetic field.


The Ohmic and ambipolar heating are determined by ionization fraction and magnetic field strength, as Equations (<ref>) and (<ref>) show. In Figure <ref>, we plot radial profiles of the ionization fraction (panel (a)) and magnetic field components (panel (b)). It is considered that B_φ strength  is limited by the equipartition value B_ eq. The equipartition value corresponds to plasma beta β=1.



Figure <ref>(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≤ 10^-12, at 0.8 au <r< 21 au. This region is called `dead' zone <cit.>. 

Figure <ref>(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_φ>(B_r, B_z) in this region. Inside the `dead' zone, the magnetic field is quasi poloidal due to efficient Ohmic diffusion, B_z≫ (B_r, B_φ). In the outer region, r>20 au, ambipolar diffusion hinders generation of B_r, and the magnetic field has quasi-azimuthal geometry, B_φ≥ B_z.



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


As it was mentioned in section <ref>, four heating sources are considered in our model of the disc: viscous heating (Γ_ v), heating by cosmic rays and radioactive decay (Γ_ cr), Ohmic (Γ_ o) and ambipolar heating (Γ_ a). 
We carried out two simulations: one with viscous and heating by cosmic rays and radioactive decay (Γ_ v + Γ_ cr), and another with additional effect of Ohmic and ambipolar heating (Γ_ v + Γ_ cr + Γ_ o+Γ_ a). In Figure <ref>, we plot radial profiles of disc temperature (panel (a)) and surface density (panel (b)) for this two cases. Figure <ref>(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≈ (0.5-1) au, and by 100 K near the outer edge, r≈ (30-50) au. Figure <ref>(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.



§ 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_φ strength in the regions of maximal ambipolar velocity. In the case of maximum B_φ=B_ eq, temperature increases by value ≤ 3000 K near the inner boundary of the `dead' zone, r≈ (0.5-1) au, and by 100 K near its outer boundary r≈ (30-50) au. If the magnetic buoyancy operates, then B_φ≈ 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.

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). 





 mhd
 




		
		





