2 D MHD Simulations of the State Transitions of X-Ray Binaries Taking into Account Thermal Conduction

Thermal conduction plays an important role in bimodal accretion flows consisting of high-temperature flow and cool flow, especially when the temperature is high and/or has a steep gradient. For example, in hard-to-soft transitions of black hole accretion flows, thermal conduction between the high-temperature region and the low-temperature region is appropriately considered. We conducted two-dimensional magnetohydrodynamic (MHD) numerical simulations considering anisotropic heat conduction to study condensation of geometrically thick hot accretion flows driven by radiative cooling during state transitions. Numerical results show that the intermediate region appears between the hot corona and the cool accretion disk when we consider heat conduction. The typical temperature and number density of the intermediate region of the 10 M black hole at 10Rg(Rg = 3.0× 106cm is the Schwarzschild radius) are 4× 1010 < T [K] < 4× 1012 and 5× 1015 < n [cm−3] < 5× 1017, respectively. The thickness of intermediate region is about half of the radius. By comparing two models with or without thermal conduction, we demonstrate the effects of thermal conduction.


Introduction
Black hole X-ray binary systems consisting of either a high-mass X-ray binaries (HMXB) such as "Cyg X-1" or a low-mass X-ray binaries (LMXB) such as "GX339-4" show two quasi-steady X-ray spectral states-hard state and soft state-and state transitions between them.X-ray satellites and other instruments have revealed the hidden nature of X-ray binaries, such as existence of the intermediate state, jet ejections and quasi-periodic oscillations [1].Theoretically, accretion flows in the hard state are explained by optically thin high-temperature accretion flows, "advection-dominated accretion flows (ADAFs)" [2,3].On the other hand, the soft state is explained by the optically thick cool accretion disks, "Shakura Sunyaev disks (SSDs)" [4].To understand the structural change of accretion flows during state transitions; however, two-dimensional or three-dimensional numerical studies are necessary.
Hard-to-soft state changes has been studied by Machida et al. by performing three-dimensional MHD simulations considering radiative cooling by bremsstrahlung [5].They found that the magnetically supported mildly cool accretion disk is formed during the vertical contraction of the disk driven by the loss of internal energy due to radiative cooling.
Galaxies 2019, 7, 22 2 of 13 Recently hydrodynamical (HD) simulations have been performed using the alpha viscosity.Das and Sharma carried out two-dimensional HD simulations taking into account bremsstrahlung and isotropic thermal conduction [6].By changing the accretion rate, they reproduced transitions from hard state to soft state and vice versa.Wu et al. carried out two-dimensional HD simulations including synchrotron and inverse-Compton scattering in addition to thermal bremsstrahlung [7].Moreover, they studied two-temperature structures of accretion flow, although they assume a simple equation of ion and electron temperatures.Based on their numerical results, they proposed two-phase accretion model for the intermediate state of black hole X-ray binaries.HD simulations with alpha viscosity, however, cannot simulate the formation of low beta (β = p/(B 2 /8π) 1) disks.Magnetohydrodynamical approach is preferable, because the origin of viscosity in accretion disk is believed to be the angular momentum transport driven by the magneto-rotational instability [8,9].
Thermal conduction is regarded as important in solar astronomy.Yokoyama and Shibata performed MHD simulations to investigate solar flares by including anisotropic thermal conduction [10][11][12].The heat released by magnetic reconnection above the chromosphere is transported downward to the chromosphere along the magnetic field line and evaporates them.Protostellar flares are also studied taking into account thermal conduction [13].Meyer and Meyer-Hofmeister and their collaborators intensively studied the role of thermal conduction in accretion disk system, such as dwarf novae and black hole accretion systems from stellar size to the size of active galactic nuclei [14,15].Once we consider the coexistence of hot corona and cool accretion disk, inevitably evaporation occurs.We must include the effect of thermal conduction, when we study the accretion disk surrounded by hot corona.
We perform two-dimensional MHD simulations including thermal conduction to study the transition from hard state to soft state, in other words, to vertical contraction of the hot accretion flow by extracting the energy by radiative cooling.Our purpose is to study the effect of thermal conduction for the coexistence of hot and cool region.In Section 2 we introduce the basic equations and initial and boundary conditions.In Section 3 we present numerical results.In Section 4 we discuss and compare the results with previous works.

Numerical Methods
We adopt the cylindrical coordinate (r, ϕ, z).The rotational axis is z axis.We assume axisymmetry, therefore partial derivatives in the azimuthal direction are set equal to zero.We use the "Coordinated Astronomical Numerical Software (CANS)" to solve the equations of resistive MHD including thermal conduction [16].The CANS integrates differential equations by the operator-splitting method dividing calculations into the MHD part and thermal conduction part [12].We adopt the modified Lax-Wendroff method with artificial viscosity for the MHD part.For thermal conduction part, we choose the Biconjugate gradient stabilized method (Bi-CGSTAB method) which is an implicit method to solve the large, sparse, unsymmetrical and linear systems of equations [17,18].The implicit method is necessary to calculate our model because the time scale of thermal conduction is much shorter than one of the other time scales in high-temperature regions.

Basic Equations
The basic equations are as follows: ∂ ∂t where ρ is the density, v is the velocity, p is the pressure, I is the unit tensor, B is the magnetic field, g = −∇ψ is the gravitational acceleration of pseudo-Newtonian potential , G is the gravitational constant, M = 10M is the black hole mass, R g = 2GM/c 2 is the Schwarzschild radius, c is the speed of light, γ = 5/3 is specific heat ratio, E is the electric field, Q hc is the heating rate due to thermal conduction, Q br is the radiative cooling rate by bremsstrahlung, k B is the Boltzmann constant, T is the temperature, µ = 0.5 is the mean molecular weight of pure hydrogen medium, m H is the hydrogen mass, α D is the dynamo alpha, η is the resistivity, j is the current density.
Let us explain Q hc , Q br , α D and η in more detail.We employ the formula derived by Spizter for the coefficient of thermal conduction [20].The heating rate of thermal conduction is written as Here κ is the coefficient of thermal conduction parallel to the magnetic field, b is the unit vector of magnetic field.We ignore κ ⊥ which is the tangential coefficient of thermal conductivity, since it is much smaller than κ .We set the upper limit of temperature to evaluate the coefficient of thermal conduction, T < 0.03T 0 .The Radiative cooling rate by bremsstrahlung is written as The parameter of dynamo alpha α D is derived from the mean field theory of galactic dynamo [21,22].
where, V S = γp/ρ is the sound speed, V A = B 2 /(4πρ) is the Alfvén velocity.According to our assumption of axisymmetry, our models cannot regenerate the poloidal components of magnetic field enough.Therefore, we refer the mean field theory and mimic the dynamo action for poloidal magnetic field.We assume anomalous resistivity [11,12].The resistivity is written as follows: where V drift = j/(en) is the drift velocity, e is the elementary electric charge, n = ρ/(µm H ) is the number density, V crit = 0.9c is the threshold of this anomalous resistivity.We assume that there is the upper limit of resistivity, η max = 4π × 10 −7 .

Initial Conditions
The initial state is a torus surrounded by low-density halo.The initial gas halo is assumed to be isothermal and non-rotating in hydrostatic equilibrium.The density of the halo is written as follows: where 3 is the halo density at (R 0 , 0), T halo = 15T 0 = 5.1 × 10 12 K is the uniform halo temperature, T 0 is the temperature of initial torus at (R 0 , 0).Initial gas torus which is embedded in the halo satisfy the polytropic relation p torus = Kρ torus γ = Kρ torus 1+1/n , K being constant, n = 1/(γ − 1) = 3/2 being polytropic index [6,7].We assume that the gas torus has the constant angular momentum.Moreover, it is in dynamical equilibrium.From the conditions of rotational equilibrium, we obtain where d is the distortion parameter [23].We employ d = 1.5.The torus has the maximum density ρ 0 at (R 0 , 0).We obtain K by substituting above values into the Equation (13).
Initial magnetic field is obtained according to the following manner.At first, we assume the azimuthal component of vector potential is proportional to the torus density.A φ = Cρ torus , C being constant.Therefore, magnetic field is calculated as, B r = −C∂A φ /∂z and B z = C(A φ /r + ∂A φ /∂r).Next, we rearrange the strength of magnetic field to satisfy plasma beta β = p torus /(B 2 /8π) = 1000 in the torus.

Grids and Boundary Conditions
We take a rectangular computation domain in the r-z plane.The region is r min < r < r max and z min < z < z max , where r min = −0.15rg , r max = 100r g , z min = −100r g and z max = 100r g .We assume that r = 0 is the rotational axis.A logarithmic grid spacing is adopted along the both r and z directions.Here 256 × 512 grids are used.At the cell-surface boundary at r = 0 between grids which are next to each other across the rotational axis, symmetric-boundary conditions are imposed.For ρ, p, v z and B z , we do not change signs at the boundary.On the other hand, for v r , v φ , B r and B φ , we change signs at the boundary.At the boundaries at r = r max , z = z min and z = z max , free-boundary conditions are imposed.Physical quantities are extrapolated from the values inside the boundary.We impose damping condition in the shell, where q is the original value, q is the value modified by the damping treatment, q init is the initial value at t = 0, f damp is the coefficient of damping which depends on the distance from the origin, r in = r g is the inner radius of the shell, r out = 2r g is the outer radius of the shell.

Results
We carried out numerical simulations of two models.One is for a model which does not include thermal conduction (Model 1).The other is for a model which includes thermal conduction (Model 2).To include the heating by dissipation of magnetic energy, we carry out simulations without including radiative cooling until t/t 0 = 10.At t/t 0 > 10, we incorporate radiative cooling in both cases.The variables are normalized by the units shown in Table 1 in calculations.Results also are shown in the units of Table 1.

Variable
Quantity Unit Value

Results of Model 1
Figure 1 shows the results of Model 1. Upper panel shows the density distribution.Middle panel is the distribution of temperature.Lower panel shows the distribution of magnetic energy density.Left column is distribution at t/t 0 = 10.The initial torus and magnetic field evolve and mass accretion proceeds.After 10 rotations at the fiducial radius r/R 0 = 1, the hot accretion flow formed is geometrically thick and has the low density.It is regarded as ADAF.Middle column is the distribution of physical quantities at t/t 0 = 30.After radiative cooling is included, the accretion flow begins to condense and settle down to the midplane (z/R 0 = 0).Hot component of ADAF disappears by t/t 0 = 15.After that, the hot corona and the cool accretion disk co-exist.The temperature and density of the cool disk is about T ∼ 10 −6 T 0 ∼ 3 × 10 5 K and n ∼ 10n 0 ∼ 5 × 10 19 cm −3 respectively.The interface between the cool disk and the hot corona has steep gradients of density and temperature.Right column shows the vertical distributions of density, temperature, and magnetic energy density at t/t 0 = 30 and of the radius r/R 0 = 0.8.Black lines show the initial distributions.Blue lines show the distribution at t/t 0 = 10.We can identify the geometrically thick ADAF whose density is low, and temperature is high.Red lines show the results at t/t 0 = 30.The density of cool accretion disk around the midplane increases.The cool accretion disk decreases the geometrical thickness.On the other hand, the temperature of the cool accretion disk decreases.The quasi-steady state is achieved in the cool accretion disk, where radiative cooling is balanced with heating by Joule heating and compression.There is the steep gradient of temperature.We expect that thermal conduction plays an important role.The region which has the large magnetic energy was confined into the cool accretion disk.In this region plasma beta become smaller than unity.The magnetic pressure is dominant in the cool accretion disk.
other hand, the temperature of the cool accretion disk decreases.The quasi-steady state is achieved in the cool accretion disk, where radiative cooling is balanced with heating by Joule heating and compression.There is the steep gradient of temperature.We expect that thermal conduction plays an important role.The region which has the large magnetic energy was confined into the cool accretion disk.In this region plasma beta become smaller than unity.The magnetic pressure is dominant in the cool accretion disk.

Results of Model 2
Figure 2 shows the results of Model 2 carried out by including thermal conduction.Upper, middle, and lower panels, and left, center, and right columns have same meaning as Figure 1.Although we include thermal conduction, almost the same hot accretion flow appears as Model 1.After the time passes t/t 0 = 10, the hot accretion flow condenses and produces a cool accretion disk inside.However, the hot accretion flow does not disappear perfectly in Model 2. The hot accretion flow changes its appearance to the intermediate region between the hot corona and the cool accretion disk.We see the sandwich structure of three components.The intermediate region becomes denser and has lower temperature than the corona surrounding the intermediate region.At this stage, the heat is transported from the corona to the cool accretion flow.Temporarily the cool accretion disk is formed; however, according to thermal conduction the disk partially evaporates gradually.In the intermediate region, magnetic field at t/t 0 = 30 is weaker than that at t/t 0 = 10; however it remains stronger than Model 1 at same time.Figure 3 shows the vertical distributions of nine physical quantities averaged between r/R 0 = 0.75 and r/R 0 = 0.85 at t/t 0 = 30.Dashed and Solid lines represent the results of Model 1 and Model 2 respectively.Figure 3a-c 3g-i show the vertical distributions of three components of magnetic field.Figure 3h shows that the azimuthal component of magnetic field is dominant.This has important meaning.If the magnetic field lines are wound tightly, it reduces the coefficient of thermal conduction in vertical direction.Effective optical depth integrated vertically is much smaller than unity, τ eff < 0.052 at r/R 0 = 0.8.The integration was done as follows: where κ es = 0.34 cm 2 g −1 is the electron scattering opacity, κ f f = 6.4 × 10 22 ρT −7/2 cm 2 g −1 is the free-free absorption opacity.The system that include intermediate region and cool accretion disk is optically thin enough.
Galaxies 2019, 6, x FOR PEER REVIEW 8 of 14 Figure 3g-i show the vertical distributions of three components of magnetic field.Figure 3h shows that the azimuthal component of magnetic field is dominant.This has important meaning.If the magnetic field lines are wound tightly, it reduces the coefficient of thermal conduction in vertical

Discussion
We solved the 2D MHD equations including the thermal conductivity.Previous works done by Das and Sharma, and Wu et al. are not MHD simulations but HD simulations using alpha viscous parameter [6,7].We did not assume alpha viscosity.We assume that the angular momentum is Galaxies 2019, 7, 22 9 of 13 transported mainly by the rϕ element of the magnetic stress in the rotating gas.We define viscous alpha parameter α MHD as where "Domain" is the region of integration, 0.4 < r/R 0 < 1.2, 0 < ϕ < 2π and −0.3 < z/R 0 < 0.3.Figure 4a shows the time evolution of viscous parameter α MHD .Black and red lines show the results of Model 1 and Model 2, respectively.Before radiative cooling is included (t/t 0 < 10), both calculations exhibit the magnitude is around α MHD ∼ 0.01.This is a typical value used in hot accretion flows [6,7].
There is no large difference between our alpha parameter and ones of previous works.Hot accretion flows do not have large difference between models with thermal conduction or not.After radiative cooling is taken into account, the cool and dense accretion disk appears and α MHD rises and drops in a short time.There is the maximum α MHD ∼ 0.1.In the early stage of transition, the magnetic field changes the structure drastically and azimuthal magnetic field becomes strong by condensation [5].
That is the reason parameter α MHD shows rapid variation.The evolution of alpha parameter is not considered in previous works.Further analysis of the angular momentum transportation will be done in future.Figure 4b shows the evolution of the luminosity normalized by the Eddington luminosity L Edd = 4πGMm p c/σ T = 1.26 × 10 39 erg/s.Here m p is the proton mass and σ T is the Thomson cross section.We calculate the luminosity as Where the "Domain" means the region as 0 < r/R 0 < 10, 0 < ϕ < 2π and −10 < z/R 0 < 10.Black and red lines show the results of Model 1 and Model 2, respectively, as well as Figure 4a.Despite the appearance of the intermediate region, there is no difference between Model 1 and Model 2.
Since mainly the cool accretion disk plays a role of the radiation of bremsstrahlung, we cannot distinguish the two models.
Figure 5a shows the radiative cooling rate of Model 2 at t/t 0 = 30.Region where radiative cooling is strong is localized in the cool accretion disk whose density is high.There is the temperature minimum at the midplane.Therefore, the radiative cooling rate at midplane become smaller than accretion flows around midplane.Figure 5b shows the thermal conduction rate of Model 2 at t/t 0 = 30.Thermal conduction works at the large temperature gradient.In this region z/R 0 ∼ 0.05, the vertical magnetic field is greater than the radial magnetic field (Figure 3d,f).Heat is transported from the intermediate region to the cool accretion disk effectively.
To compare the efficiency of thermal conduction and radiative cooling, we computed the Field length λ F without including the heating term [24,25].When the Field length is longer than the scale of the medium, thermal conduction dominates radiative cooling.On the other hand, when the Field length is shorter than the scale of the medium, radiative cooling dominates thermal conduction.Our numerical results show the magnetic field changes temporally and spatially in the r − z plane.Since magnetic turbulence develops in the disk, we simply assume that heat conduction isotropic in the r − z plane.The Field length computed without including the heating term is Figure 6a,b show the distribution of the Field length λ F /R 0 of Model 2 at t/t 0 = 10 and t/t 0 = 30, respectively.In the coronal region at t/t 0 = 10, since the Field length is longer than the size of the simulation region, thermal conduction dominates radiative cooling.However, in the region r > z, since radiative cooling becomes dominant over heat conduction, hot accretion flows can condensate by the cooling instability.On the other hand, thermal conduction suppresses the cooling instability in the intermediate region at t/t 0 = 30.Figure 7 shows the evolutions of mass and volume whose medium is in three temperature ranges in the cylinder defined as 0 < r/R 0 < 1.6, 0 < ϕ < 2π and −1.6 < z/R 0 < 1.6 except for the sphere √ r 2 + z 2 /R 0 < 0.3.Three temperature range are T/T 0 > 10 (denoted by red line), 0.1 < T/T 0 < 10 (denoted by green line) and T/T 0 < 0.1 (denoted by blue line), these temperature ranges are supposed to be the hot corona, the intermediate region and the cool accretion disk, respectively.Figure 7a,b show that the intermediate region remains by thermal conduction against radiative cooling.Although the mass of intermediate region is small compared to cool accretion disk, the volume of intermediate region is large compared to other temperature region (Figure 7c,d).These situations and sandwich structure of three components is helpful to produce the hard X-ray spectral component by inverse-Compton processes.Post process calculations by a Monte Carlo simulation for X-ray spectrum will be presented in subsequent papers.
In this study, we considered radiative cooling by bremsstrahlung.Synchrotron radiation and inverse-Compton scattering can change the structure of accretion disks in the equatorial region and the intermediate region.Das and Sharma pointed out the possibility that synchrotron radiation and inverse-Compton scattering have same dependency on the density as bremsstrahlung [6].We think that increase of the density may compensate the lack of radiative cooling.When we consider radiative cooling by synchrotron radiation and inverse-Compton scattering, the formed cool accretion disk has lower temperature and higher number density than the results of our present study.Inverse-Compton scattering reduces the thickness of the intermediate region.On the other hand, thermal conduction suppresses the shrinks of the intermediate region.We will perform the parameter survey of density in future.
We set the upper limit of temperature to evaluate the coefficient of thermal conductivity, since we solve one-temperature MHD equations.We suppose that the upper limit of temperature is electron temperature in two-temperature accretion flows [2].It is also useful to avoid the unusual thermal conductivity in high-temperature and low-density corona.As is well known, ADAFs has the two-temperature structures.The electron temperature is suppressed to be lower than the ion temperature since the energy exchange rate between ions and electrons by the Coulomb coupling is low due to the low density and the short free-fall time scale.It is preferable to simulate the two-temperature MHD to consider the state transition when the intermediate region exists.Wu et al. consider the two-temperature HD models; however, they do not calculate the two energy equations for ions and electrons independently [7].We are seeking to update our codes which solve two-temperature MHD equations.
that increase of the density may compensate the lack of radiative cooling.When we consider radiative cooling by synchrotron radiation and inverse-Compton scattering, the formed cool accretion disk has lower temperature and higher number density than the results of our present study.Inverse-Compton scattering reduces the thickness of the intermediate region.On the other hand, thermal conduction suppresses the shrinks of the intermediate region.We will perform the parameter survey of density in future.We set the upper limit of temperature to evaluate the coefficient of thermal conductivity, since we solve one-temperature MHD equations.We suppose that the upper limit of temperature is electron temperature in two-temperature accretion flows [2].It is also useful to avoid the unusual thermal conductivity in high-temperature and low-density corona.As is well known, ADAFs has two-temperature structures.The electron temperature is suppressed to be lower than the ion temperature since the energy exchange rate between ions and electrons by the Coulomb coupling is low due to the low density and the short free-fall time scale.It is preferable to simulate the twotemperature MHD to consider the state transition when the intermediate region exists.Wu et al.The state transition is driven by the variation of the mass accretion rate.However, we do not include the mass supply calculation in our codes.We studied the various initial condition of density, 1.0 × 10 −8 < ρ 0 [g/cm 3 ] < 1.0 × 10 −4 .Instead of the mass supply, we choose the critical initial condition of the density which can make the hot accretion flow condense.It is easy to obtain the results of the MHD simulation of the state transition.Machida et al. and Wu et al. adopt this concept [5,7].However, it is preferable to implement the subroutine of the mass supply to achieve the cycle of the state transition in simulations by controlling the mass accretion rate [6].This difficulty should be solved in future.

Figure 1 ..
Figure 1.Results of Model 1.(a) Distribution of number density at

Figure 1 .
Figure 1.Results of Model 1.(a) Distribution of number density at t/t 0 = 10 before cooling term is included.(b) Distribution of number density at t/t 0 = 30.(c) Vertical distribution of number density at three different time and at r/R 0 = 0.8.Black line shows initial condition at t/t 0 = 0. Blue line shows result at t/t 0 = 10.Red line shows result at t/t 0 = 30.(d) Distribution of temperature at t/t 0 = 10.(e) Distribution of temperature at t/t 0 = 30.(f) Vertical distribution of temperature.Color scale is the same as (c).(g) Distribution of magnetic energy density (B/B 0 ) 2 /8π at t/t 0 = 10.(h) Distribution of magnetic energy density at t/t 0 = 30.(i) Vertical distribution of magnetic energy density.Colors have the same meaning as (c) and (f).
disk.We see the sandwich structure of three components.The intermediate region becomes denser and has lower temperature than the corona surrounding the intermediate region.At this stage, the heat is transported from the corona to the cool accretion flow.Temporarily the cool accretion disk is formed; however, according to thermal conduction the disk partially evaporates gradually.In the intermediate region, magnetic field at  remains stronger than Model 1 at same time.

Figure 2 .
Figure 2. Results of Model 2. Same quantities and conditions of visualization as Figure 1.

Figure 3
Figure 3 shows the vertical distributions of nine physical quantities averaged between 0 0.75 rR and

Figure 2 .
Figure 2. Results of Model 2. Same quantities and conditions of visualization as Figure 1.
Figure3shows the vertical distributions of nine physical quantities averaged between r/R 0 = 0.75 and r/R 0 = 0.85 at t/t 0 = 30.Dashed and Solid lines represent the results of Model 1 and Model 2 respectively.Figure3a-cshow the number density, the temperature, and the pressure, respectively.The pressure does not vary very much in the intermediate region.Therefore, the density decreases with height (∂ ln n/∂ ln z ∼ −2) in Model 2. Meanwhile, the temperature increases with height (∂ ln T/∂ ln z ∼ 2).In Model 1, the density and temperature show steep gradient.The difference between Model 1 and Model 2 indicate the existence of the intermediate region.The number density and temperature of the intermediate region are 5 × 10 15 < n [cm −3 ] < 5 × 10 17 and 4 × 10 10 < T [K] < 4 × 10 12 , respectively.The thickness of intermediate region is about 0.5R 0 .Figure 3d-f show the components of velocity.The azimuthal component of velocity simply shows the Keplerian rotation from the cool accretion disk to the outer corona.The radial and vertical components show the exchange of their signs.It means that there is inflow and outflow component of the intermediate region and the lower corona.Figure3g-ishow the vertical distributions of three components of magnetic field.Figure3hshows that the azimuthal component of magnetic field is dominant.This has important meaning.If the magnetic field lines are wound tightly, it reduces the coefficient of thermal conduction in vertical direction.Effective optical depth integrated vertically is much smaller than unity, τ eff < 0.052 at r/R 0 = 0.8.The integration was done as follows:

Figure
Figure3shows the vertical distributions of nine physical quantities averaged between r/R 0 = 0.75 and r/R 0 = 0.85 at t/t 0 = 30.Dashed and Solid lines represent the results of Model 1 and Model 2 respectively.Figure3a-cshow the number density, the temperature, and the pressure, respectively.The pressure does not vary very much in the intermediate region.Therefore, the density decreases with height (∂ ln n/∂ ln z ∼ −2) in Model 2. Meanwhile, the temperature increases with height (∂ ln T/∂ ln z ∼ 2).In Model 1, the density and temperature show steep gradient.The difference between Model 1 and Model 2 indicate the existence of the intermediate region.The number density and temperature of the intermediate region are 5 × 10 15 < n [cm −3 ] < 5 × 10 17 and 4 × 10 10 < T [K] < 4 × 10 12 , respectively.The thickness of intermediate region is about 0.5R 0 .Figure 3d-f show the components of velocity.The azimuthal component of velocity simply shows the Keplerian rotation from the cool accretion disk to the outer corona.The radial and vertical components show the exchange of their signs.It means that there is inflow and outflow component of the intermediate region and the lower corona.Figure3g-ishow the vertical distributions of three components of magnetic field.Figure3hshows that the azimuthal component of magnetic field is dominant.This has important meaning.If the magnetic field lines are wound tightly, it reduces the coefficient of thermal conduction in vertical direction.Effective optical depth integrated vertically is much smaller than unity, τ eff < 0.052 at r/R 0 = 0.8.The integration was done as follows:

Figure 3 ...
Figure 3. Results of Model 1 and Model 2 at

Figure 3 .
Figure 3. Results of Model 1 and Model 2 at t/t 0 = 30.Dashed lines represent the result of Model 1. Solid Lines represent the results of Model 2. These panels show the vertical distributions of the nine physical quantities averaged between r/r 0 = 0.75 and r/r 0 = 0.85.Black lines indicate positive values.Red lines indicate negative values, that is, mean the absolute values of quantities.(a) Number density.(b) Temperature.(c) Pressure.(d) Radial component of velocity.(e) Azimuthal component of velocity.(f) Vertical component of velocity.(g) Radial component of magnetic field.(h) Azimuthal component of magnetic field.(i) Vertical component of magnetic field.

Figure 4 .
Figure 4. Evolutions of the averaged magnetic stress and integrated luminosity.Black lines show the results of Model 1 without thermal conduction.Red lines show the results of Model 2 including thermal conduction.(a) Averaged viscous alpha parameter

Figure 4 .Figure 5 .
Figure 4. Evolutions of the averaged magnetic stress and integrated luminosity.Black lines show the results of Model 1 without thermal conduction.Red lines show the results of Model 2 including thermal conduction.(a) Averaged viscous alpha parameter

Figure 5 .
Figure 5. (a) Distributions of radiative cooling rate −Q br /Q 0 of Model 2 at t/t 0 = 30.Q 0 = ρ 0 V 0 3 /R 0 = 5.1 × 10 16 erg/ cm 3 • s is normalization unit. (b) Distributions of thermal conduction rate the Field length is longer than the size of the simulation region, thermal conduction dominates radiative cooling.However, in the region rz  , since radiative cooling becomes dominant over heat conduction, hot accretion flows can condensate by the cooling instability.On the other hand, thermal conduction suppresses the cooling instability in the intermediate region at 0 / 30 tt .

Figure 6 .
Figure 6.Distributions of the Field length F 0 / R  of Model 2. (a) Result at

Figure 7 Figure 6 .
Figure 7 shows the evolutions of mass and volume whose medium is in three temperature ranges in the cylinder defined as  0 0 / 1.6 rR , 02   and    0 1.6 / 1.6 zR except for the sphere

Figure 7 ..
Figure 7. Time evolution of mass and volume of Model 1 and Model 2. Results are obtained by the integration in the domain of  0 0 / 2 rR , 02  

Figure 7 .
Figure 7. Time evolution of mass and volume of Model 1 and Model 2. Results are obtained by the integration in the domain of 0 < r/R 0 < 2, 0 < ϕ < 2π and −1 < z/R 0 < 1 except for the region of √ r 2 + z 2 /R 0 < 0.3.These are normalized by initial total mass and volume of the above domain, respectively.(a) Model 1. Blue line shows the mass of gas whose temperature range is T/T 0 < 0.1.Green line shows the mass of 0.1 < T/T 0 < 10.Red line shows the mass of 10 < T/T 0 .Black line shows sum of three components.(b) Model 2. Same as (a).(c) Model 1.Time evolution of volume of each temperature range same as (a).Colors mean same as (a).(d) Model 2. Same as (c).