Local models of two-temperature accretion disc coronae. II. Ion thermal conduction and the absence of disc evaporation

We use local stratified shearing-box simulations with magnetic field-aligned thermal conduction to study an idealized model of the coupling between a cold, radiatively efficient accretion disc, and an overlying, hot, two-temperature corona. Evaporation of a cold disc by conduction from the hot corona has been proposed as a means of mediating the soft-to-hard state transitions observed in X-ray binary systems. We model the coronal plasma in our local disc patch as an MHD fluid subject to both free-streaming ion conduction and a parameterized cooling function that captures the collisional transfer of energy from hot ions to colder, rapidly cooling leptons. In all of our models, independent of the initial net vertical magnetic flux (NF) threading the disc, we find no evidence of disc evaporation. The ion heat flux into the disc is radiated away before conduction can heat the disc's surface layers. When an initial NF is present, steady-state temperature, density, and outflow velocities in our model coronae are unaffected by conduction. Instead of facilitating disc evaporation, thermal conduction is more likely to feed the disc with plasma condensing out of the corona, particularly in flows without NF. Our work indicates that uncertainties in the amount of NF threading the disc hold far greater influence over whether or not the disc will evaporate into a radiatively inefficient accretion flow compared to thermal conduction. We speculate that a change in net flux mediates disc truncation/evaporation.


INTRODUCTION
Thermal conduction from hot to cold gas can act to evaporate the cold medium, driving an evaporative flow (Cowie & McKee 1977;McKee & Cowie 1977;Doroshkevich & Zel'dovich 1981;Giuliani 1984;Draine & Giuliani 1984;Balbus 1986).In the absence of cooling and gravity, the thermal energy flux provided by conduction into the cold medium must be balanced by the enthalpy flux of the ensuing evaporative outflow (Draine 2011).These ideas, dating back to work on the solar corona by Pikel'ner (1948Pikel'ner ( , 1950)), have found widespread application to a diverse range of astrophysical phenomena, from the destruction of cold clouds embedded in hot supernova remnants to accretion discs around white dwarfs and black holes.
In this paper, our focus is on accretion discs in X-ray binary systems (XRBs) and active galactic nuclei (AGN).XRBs exhibit a variety of spectral states, most notably, the 'high/soft' and 'low/hard' states, where soft and hard refer to the spectral 'hardness' of power-law emission in the X-ray band (see Remillard & McClintock (2006) and Done, Gierliński & Kubota (2007) for reviews).
Transitions from the soft to the hard state are often described via the disc truncation model of Esin, McClintock & Narayan (1997).In this model, the luminous soft state is associated with a standard thin accretion disc (Shakura & Sunyaev 1973) extending down to the innermost stable circular orbit (ISCO) of a stellar mass black hole, ★ E-mail: cbambic@princeton.edu while the low-luminosity hard state is associated with a truncated disc feeding a hot, optically thin, radiatively inefficient accretion flow (RIAF; Ichimaru 1977;Rees et al. 1982;Narayan & Yi 1994).A hot (electron temperature   ∼ 10 9 K) corona overlying the thin disc in the soft state, as would be expected in the observationally motivated two-phase model of a corona 'sandwiching' an accretion disc (Haardt & Maraschi 1991), provides a reservoir of thermal energy that can be tapped to evaporate the thin disc and form a RIAF (Liu et al. 1999).The evaporation itself occurs through thermal conduction from the corona into the disc (Meyer, Liu & Meyer-Hofmeister 2000), an idea first invoked to describe salient features in the phenomenology of dwarf novae in accreting white dwarf stars, i.e., cataclysmic variable stars (Meyer & Meyer-Hofmeister 1994, hereafter, MM94).If the heat flux from the corona into the disc is large enough to overcome cooling, disc evaporation and truncation can ensue.Thus, thermal conduction offers a physical mechanism for mediating disc truncation, and subsequently, state transitions, in XRBs.
Conductive evaporation has been extended to XRB and AGN accretion flows by Meyer-Hofmeister & Meyer (1999) and Liu et al. (1999); however, even with the recognition that coronae in accretion flows around black holes should be two-temperature, with the ion temperature much larger than that of the electrons (Di Matteo, Blackman & Fabian 1997), the vast majority of works on XRBs and AGN have restricted their attention to electron thermal conduction (Różańska & Czerny 2000;Liu et al. 2002), rather than considering the heat flux carried by ions.This assumption is justified in the colder transition region just above the disc; however, ion thermal conduction is certainly critical in the hot, weakly collisional corona (Spruit & Deufel 2002;Dullemond & Spruit 2005), and even the temperature profiles of Liu et al. (2002) imply that ion thermal conduction would dominate over that of the electrons.
Electron thermal conduction has been shown to be inefficient at inducing evaporation within a few tens of gravitational radii   of a central accreting compact object, where   ≡  / 2 and  is the mass of the compact object.At distances of ≳10 3   , evaporation is far more efficient, and the mass-outflow rate provided through conductive evaporation can outpace accretion.Equivalently, since the surface of a white dwarf with a mass  ≈ 1.2M ⊙ is at ≈10 3   , conduction can induce evaporation in the inner regions of white dwarf discs.Cho & Narayan (2022, hereafter CN22) recently revisted this problem, analyzing the interplay of electron thermal conduction and bremmstrahlung cooling in accretion flows around black holes.Similar to previous works, these authors find that evaporation is inefficient at small radii near the black hole.Instead, conduction works to cool coronae, forming a condensing flow onto the disc.This process of conductive condensation could re-form a disc from a RIAF near the ISCO and influence the hard-to-soft transition in XRBs (Liu et al. 2007;Meyer, Liu & Meyer-Hofmeister 2007).
It is unclear whether a mechanism of disc truncation that operates at large radii alone can explain the observed properties of XRBs.Observations clearly indicate that the disc truncates at a few tens of   in the hard state of Cygnus X-1, the best studied of the XRBs (Zdziarski et al. 1999;Churazov et al. 2001); however, a number of other XRBs may possess a cold disc near the ISCO, even in the hard state (Miller et al. 2006;Kara et al. 2019).If discs are in fact truncated on scales of ∼ few × 10 R g , then a physical mechanism beyond electron thermal conduction is necessary to induce disc evaporation and truncation in XRB and AGN accretion flows.This paper explores the conductive coupling between a cold disc and hot corona to determine if evaporation can occur in the innermost regions of XRB and AGN accretion discs.Because our work is focused on the innermost few gravitational radii of these accretion flows, the thermodynamic regime explored here is that of a two-temperature ion-lepton plasma.In Paper I of this series, we demonstrated that, in this regime, the relevant cooling mechanism for the ions is Coulomb collisions with rapidly Compton-cooled leptons, and ion, rather than electron, thermal conduction in the saturated, free-streaming limit dominates the transport of heat (Bambic, Quataert & Kunz 2024, hereafter, BQK24).These thermodynamic processes, distinct from those studied by CN22, offer a possible means for truncating XRB and AGN discs close to the central object.
While conductive evaporation and disc truncation are the main motivations for this work, conduction may produce a variety of other observable effects in black hole accretion flows that can be explored by our models.Spruit & Deufel (2002) contend that the hard X-ray continuum emission characteristic of coronae results from the bombardment of electrons by MeV protons in a two-temperature RIAF.Ion conduction is a natural way to channel these hot ions into the X-ray emitting electrons.Conduction may affect the rate of energy release in coronae.Goodman & Uzdensky (2008) argue that the rate of magnetic reconnection may be regulated by conductive evaporation, which acts to modify the plasma density in the corona.Similar processes may occur on the Sun (Uzdensky 2007), where nearly half of the thermal energy released in a flare leaves the reconnecting X-point via conduction.Finally, there have been claims that a truncated disc can still produce a detectable iron line at the ISCO (Tanaka et al., 1995;Reis, Fabian & Miller 2010) if colder, denser clumps can form within a RIAF (Liska et al. 2022).However, just as conduction may act to destroy cold clouds within supernova remnants (Cowie & McKee 1977;McKee & Cowie 1977), so too may conduction erase this signal.
This paper is organized as follows.In §2, we introduce models for two-temperature cooling and free-streaming thermal conduction relevant for the innermost regions of XRB and AGN accretion flows.We then apply these models to a suite of vertically stratified magnetohydrodynamic (MHD) shearing-box simulations (Stone et al. 1996) that, for the first time, include field-aligned free-streaming ion thermal conduction.The set-up of these simulations is described in §3, while the results are presented in §4.In §5, we explore a series of toy models in one dimension to better understand our results, namely, that the inclusion of saturated ion thermal conduction cannot evaporate our model thin accretion discs.We discuss the implications of our results to state transitions in §6, and we summarize and conclude in §7.

CONDUCTION IN A TWO-TEMPERATURE CORONA
We begin by introducing models for free-streaming ion thermal conduction and Coulomb collisional cooling relevant for the thermodynamic properties of plasma in the innermost regions of XRB and AGN accretion flows.Using these models, we derive an order-of-magnitude estimate for the conditions under which conductive heating can win out over two-temperature cooling to induce evaporation of an accretion disc.These estimates motivate three-dimensional MHD simulations that can evaluate if conductive evaporation can ever be realized.
Ions exchange energy with leptons on the temperature equilibration timescale,  eq , defined in BQK24 as where the lepton temperature Θ  ≡  B   /   2 is defined in terms of the electron rest mass    2 , the lepton number density is   , and we have introduced the electron Coulomb logarithm Λ  .A simple estimate based on observational implications that the optical depth to electron scattering  es is ∼O (1) in coronae implies that the lepton number density is   ≈ 10 17 ( BH /10 M ⊙ ) −1 cm −3 .
We describe the cooling of ions through a simple, optically thin cooling function Λ of the form where E int = (3/2)   B  is the internal energy of the plasma,   is the ion number density, and  is the ion temperature, with  ≫   (Di Matteo, Blackman & Fabian 1997).Observations constrain the electron temperature and the optical depth in the corona, such that the product of the equilibration time and the Keplerian orbital frequency, Ω ≈ 6.4 × 10 2 (/10 M ⊙ ) −1 /10   −3/2 s −1 , is given by (3) Here,  is the radial size of the corona, while  is an O (1) (possibly O (/), where  is the disc scale height) factor set by the unknown geometry of the corona.Based on the optical depths implied by observations, Coulomb collisions should be weak enough to allow ions to heat up to an order unity fraction of the virial temperature , where   is the proton mass.At 10  g , this ion temperature is   ≈ 31 MeV ≈ 4 × 10 11 K. Magnetic-field-aligned ion thermal conduction in the free-streaming limit can be captured via a heat flux of the form (Malone, McCrory & Morse 1975) Here, the mass density of the coronal ions is , the speed of sound   is defined via  2  =  B /  , the thermal gradient length scale is defined as   ≡ −(d ln /d) −1 , the unit vector aligned with the magnetic field  is b ≡ /||, and we work in a local coordinate system in our disc patch such that  is the vertical coordinate with unit direction ẑ.
The free-streaming heat flux is appropriate when ion collisions are sufficiently infrequent that the Spitzer (1962) collisional heat flux saturates at the free-streaming limit.We can write the Spitzer heat flux as where  =    B  is the ion pressure and the ion-ion collision time   is longer than the equilibration time (1) by a factor of ∼(  /  ) 3/2 (  /  ) 1/2 .For our estimate of   ≃  virial , we find that  ii / eq ≈ 130.The Spitzer heat flux saturates for any thermal gradient length scale satisfying where  ,c ≡ √ 2  /Ω is the vertical scale height of the corona.Thus, even small temperature gradients with large |  | will still saturate the heat flux at the free-streaming value.Similarly, the large ion-to-lepton temperature ratio implies that, unlike for the scenarios studied in MM94, Liu et al. (2002), Liu et al. (2007), Meyer, Liu & Meyer-Hofmeister (2007), and CN22, electron thermal conduction is negligible.Even at the saturated limit, ion thermal conduction is larger than that of electrons by the same factor of (  /  ) 3/2 (  /  ) 1/2 ≈ 130, such that ion conduction carries the bulk of the conductive heat flux within a two-temperature corona.
For conduction to evaporate an accretion disc, the energy released in the disc must overcome two-temperature cooling and heat the disc, driving a substantial enthalpy flux into the corona.Conduction overcomes two-temperature cooling when Here, we have introduced a geometric suppression factor   to account for the fact that conduction across the magnetic field is strongly suppressed relative to conduction along the magnetic field (see Paper I).In principle, this condition could be satisfied in the corona where Ω eq ≳ 1, so long as the suppression factor is not ≪1.However, in the disc's surface layers, where densities are large enough that Ω eq ≪ 1, Coulomb collisions will rapidly cool the ions before evaporation can occur.To determine the temperature gradient length scales   , suppression factors   , and Coulomb cooling timescales Ω eq in an accretion disc, we turn to numerical simulations.

METHODS
The simulation set-up, initial conditions, and boundary conditions imposed on the MHD variables are identical to those presented in BQK24.To implement conduction, we must include the conductive heat flux  con in the MHD energy equation and evaluate an equation for the evolution of  con .Here, we describe a two-moment method, as presented in Jiang & Oh (2018) and originally applied to the problems of radiation transport and cosmic-ray MHD, which we use to evolve the free-streaming heat flux in our simulations.We can then solve the full system of equations for MHD, thermal conduction, and cooling within a local approximation for an accretion disc, the stratified shearing box (Hawley, Gammie & Balbus 1995), using the Athena++ MHD code (Stone et al. 2020).The shearing box 'zooms in' on a local patch of plasma orbiting a black hole at a radius  0 and the Keplerian angular velocity  = Ωẑ through mapping the global coordinates (, , ) to Cartesian coordinates  =  −  0 ,  =  0 ( − Ω), and  = , where the simulation evolves with time .

The Coulomb cooling function
Rather than evolve separate energy equations for the ions and leptons, we treat the ions as a single MHD fluid subject to cooling via Coulomb collisions.Our goal is to model properly the approximately virialized ions in the surface layers of the disc while keeping the mid-plane of the simulation at a fixed temperature and thus, scale height.We have already introduced a cooling function in Equation 2. By ignoring the Coulomb logarithm, which depends only weakly on   and   , we arrive at a cooling function that can be implemented in our simulations, Here,  B =   = 1, density  and temperature  are measured relative to their mid-plane values,  0 and  0 respectively, and the uncertain lepton physics is absorbed into a constant free parameter A, which we vary in our model (BQK24).For weak Coulomb coupling, i.e., A ∼ 10 − 10 2 , this form for the cooling function ensures that the cooling time will be very short in the bulk of the disc, where  ≃  0 ; however, the cooling time can become long relative to the orbital timescale in the diffuse surface layers of the disc where  ≪  0 , consistent with the properties of ions in the corona.

Equations solved
Traditional methods for thermal conduction suppose a closure relation for  con and treat −∇ •  con as a source term in the total energy equation.In the free-streaming regime relevant to ion conduction in accretion disc coronae, standard explicit methods for computing this source term introduce spurious, large-amplitude oscillations in the heat flux, which can only be controlled through a rather stringent condition on the timestep, Δ ∼ (Δ) 3 (Sharma, Colella & Martin 2010).Such an onerous constraint is computationally prohibitive.The two-moment method that we employ (Jiang & Oh 2018;Tan et al. 2021) circumvents this timestep constraint by introducing a maximum signal velocity  max , which controls the timestep via the standard Courant condition.Conduction acts to modify the internal energy of the MHD fluid through This is the first of the moment equations.Assuming a closure for the equilibrium, time-steady heat flux along field lines, 1965), where  is a conductivity that depends on local plasma conditions, the heat flux can be evaluated via solution of the second moment equation, 1 This moment equation reduces to the desired form of the heat flux,  con = − b • ∇ b, in the limit Crucially, Equation (11) describes the saturated, free-streaming limit of thermal conduction when  is chosen such that where  min = 10 −10  0  2  Ω is a negligibly small conductivity included for numerical stability.All simulations with free-streaming thermal conduction use this expression for , and the two-moment equations (10 and 11) are solved to evaluate the heat flux.We demonstrate the accuracy of this method in Appendix A.

Initial Conditions and Boundary Conditions
We adopt units in which the Keplerian orbital frequency in our local disc patch is Ω = 1, the mid-plane density  0 = 1, and the mid-plane temperature  0 = 1/2.The scale height of the disc at the mid-plane,   , is determined through vertical force balance: while hydrostatic equilibrium with the local potential Φ determines the initial density distribution.Pressure is determined by assuming the atmosphere is initially at a single temperature  0 .Time is measured relative to the orbital timescale  orb = 2/Ω = 2.The conductive heat flux is initially set to  con = 0 throughout the domain, corresponding to the fact that we have assumed an initially uniform-temperature atmosphere.To seed the magnetorotational instability (MRI; Balbus & Hawley 1991) in our simulations, we initialize the domain in the region || ≤ 0.5  with randomly distributed, spatially uncorrelated velocity and adiabatic pressure fluctuations with maximum amplitude  = 5 × 10 −3  s and / = 0.025 (Hawley et al. 1995).These small perturbations introduce a small heat flux initially in the domain; however, cooling and conduction rapidly smooth out these initial temperature fluctuations, and our results are independent of the initial values for the heat flux.
While the boundary conditions of the density , pressure , velocities , and magnetic field  are the same as those in BQK24, we must impose an additional boundary condition on  con .We choose a simple outflow boundary condition, where the conductive heat flux is copied into the ghost zones with zero gradient.Note that this boundary condition is inconsistent with the condition on temperature, which is simply copied into the ghost zones from the last live zone such that  (, , ± ghost ) =  (, , ± max ), where ± ghost are the -coordinates in the ghost zones.We find, though, that the choice of a simple zero-gradient boundary condition ensures the correct behaviour of the heat flux near the upper/ lower boundaries.
The influence of magnetic fields is described in terms of the net vertical magnetic flux threading the disc, For the zero-net-flux (ZNF) simulations, we initialize the field of Miller & Stone (2000), and for the net flux (NF) simulations, we initialize a constant vertical magnetic field,   .The magnetic-field strength is determined by the initial mid-plane plasma  parameter at time  = 0:  0 ≡ 2 0 ( = 0)/| 0 ( = 0)| 2 .For the ZNF simulations, we study cases where  0 = 10 2 (weak ZNF) and 10 (moderate ZNF); for NF fields, we analyze cases with  0 = 10 4 (weak NF) and 10 3 (moderate NF).

Numerical Details
Because the shearing box loses mass via winds, we inject material to maintain an approximately constant steady-state disc mass.This mass injection ensures that the the density below the corona remains high enough to keep cooling and heating roughly balanced throughout the disc.Cooling and mass injection are handled in the same way as was done in Paper I. Our cooling function ( 9) is evaluated through the exact integration method presented in Townsend (2009).Mass is injected at a rate  src according to the initial density profile to compensate for the mass-outflow rate of the wind,  wind , When we inject mass into a given cell, we maintain the velocity  and temperature  already in that cell such that we are injecting mass, momentum, and energy into the domain and maintaining the constant mass of the box,  box .These source terms as well as the cooling loss term are implemented in operator-split fashion.Our simulations are advanced in time using second-order-accurate van Leer time integration (vl2; van Leer 1979) and spatial integration is performed using the Harten-Lax-van Leer Discontinuity (HLLD) Riemann solver with second-order-accurate piecewise-linear-method (PLM) reconstruction.As in Paper I, we do not use orbital advection, and we use the smoothed vertical gravitational potential from Davis et al. (2010) (see BQK24 Equation 23).
Simulations with conduction are more expensive than those without conduction because a large value of  max must be chosen to ensure that the equilibrium heat flux is realized, i.e., to satisfy equation 12.This large signal speed must be resolved by a short time step.To cut down computational costs, we use  floor / 0 = 10 −4 for the NF simulations and  floor / 0 = 10 −6 for the ZNF simulations, with the pressure floor set as  floor =  floor  0 .Further, while we use the same domain size as BQK24, i.e., (  ,   ,   ) = (4  , 8  , 12  ), which allows for a comparison of these conduction runs with our earlier results, we run the more expensive conduction simulations at half BQK24's resolution: 16 grid cells per   .We perform higher resolution, 32 cells/  , simulations of three of our models: the moderate NF A = 10, weak NF A = 10, and moderate ZNF A = 10 models, to verify that our conclusions do not depend on the chosen resolution.All simulations are evolved for 100 orbits until time Ω ≈ 628.

Choice of 𝑣 max
The two-moment method used to evolve the conductive heat flux requires a choice of the maximum signal speed  max , which sets the size of the time-derivative term in the evaluation of the heat flux.This term is negligible so long as Equation 12is satisfied.To determine a sensible choice for  max , we leveraged intuition gained from the simulations in BQK24 to realize that, in those simulations, the Courant condition is set by the large Alfvén speed in the corona: Using the simulation data from BQK24, we found Conduction NF A = 10  0 = 10 4 1.6 14 ± 4 2.5 +0.8 −0.9 6.1 +1.7 −1.8 4.2 +0.8 −0.9 0.08 ± 0.07 Conduction NF A = 10  0 = 10 4 HiRes 1.5 10 ± 2 1.9 ± 0.6 5.1 ± 0.5 4.7 7.1 ± 0.7 3.3 ± 0.4 0.17 ± 0.05 No Conduction NF A = 10 3  0 = 10 4 3.5 -6.5 +1.8 −1.9 5.4 +1.2 −1.1 3.9 ± 0.7 0.17 Table 1.Summary of key diagnostics for all 19 simulations.Quantities are quoted with 1 error bars, where the errors are computed from the 16th and 84th percentiles of the respective time series from 30 to 100 orbits in net flux (NF) simulations, from 30 to 50 orbits in 'No Conduction' zero net flux (ZNF) simulations, and from 70 to 100 orbits in the 'Conduction' ZNF runs.Horizontally averaged temperature first rises above the mid-plane temperature  0 at height   .The box-averaged turbulent injection power  turb is measured in code units, while the cooling within the corona (  cor cool ) is measured relative to  turb .
that the horizontally averaged Alfvén velocity peaks at around 5 ,0 in the corona, i.e., ≈3.5   Ω.We choose a maximum velocity that greatly exceeds this averaged Alfvén velocity,  max = 30   Ω, to ensure accurate computation of the heat flux.We further explore the effect of  max on our conclusions in Appendix B.

Simulations and Diagnostics
We run a series of 19 simulations as detailed in Table 1.These simulations scan over two values of the Coulomb coupling parameter A ∈ {10, 10 3 } and four different field configurations: weak NF ( 0 = 10 4 ), moderate NF ( 0 = 10 3 ), weak ZNF ( 0 = 10 2 ), and moderate ZNF ( 0 = 10).For each combination of A and field configuration, we provide an analysis of runs with and without field-aligned, free-streaming thermal conduction.
For any simulations with NF, we run full simulations from 0 to 100 orbits with thermal conduction either always on or always off.However, in simulations initialized with ZNF, we found that early on at times  < 5 orbits, conduction acted to collapse the upper atmosphere onto the disc before a corona could properly form.This collapse evacuated the atmosphere above || = 2  , resulting in the continual addition of floored material at large ||.Subsequently, the timestep dropped dramatically.To avoid this initial collapse, we run all ZNF simulations from 0 to 50 orbits with conduction effectively turned off at the level of  =  min (denoted as 'No Conduction' in Table 1).Then, at  = 50 orbits, conduction is turned on.As we show in §4.5, a condensing inflow onto the disc forms, and a steady state is reached by 70 orbits.Thus, a single ZNF simulation provides a 'No Conduction' run between 0 and 50 orbits and a 'Conduction' run from 50 to 100 orbits.All of our simulations are run at a resolution of 16 grid cells/  with the exception of the 3 runs denoted as 'HiRes' in Table 1, which are run at a resolution of 32 cells/  .
To assess if conductive evaporation can mediate the destruction of the accretion disc, we define an evaporation time  evap as where the surface integral evaluates the peak conductive flux into the disc measured at the point where the flux into the disc is maximal, a height that we denote   .Note that |  | is not necessarily the same above the below the mid-plane due to symmetry breaking in the shearing box.This definition of the evaporation time  evap measures the evaporative outflow rate  evap for the maximum time-averaged heat flux into the disc-a 'best case scenario' for the timescale over which the disc can evaporate in the absence of cooling.Similarly, the wind depletion time is given by Energy is injected into the domain at the rate (Hawley et al. 1995) where  ≡ −d ln Ω/d ln  = 3/2 is the shear parameter for a Keplerian flow,   is the radial extent of the domain, and the surface integral is taken over the shearing boundary.Energy injection is modulated by the - component of the turbulent stress tensor, where the turbulent velocity  =  + 3 2 Ω ŷ does not include the background Keplerian shear.
Thermal equilibrium, i.e., a balance of heating and cooling, is established in the disc on the thermal timescale, where we have introduced the Shakura & Sunyaev (1973) '' parameter, which we define in terms of the mid-plane pressure  0 such that  ≡ ⟨T   / 0 ⟩   .We define  mid ≡ (|| < 2  ).Note that this definition of  mid is different than that used in BQK24.For the lower resolution 16/  simulations explored in this work, we find that there is increased numerical dissipation at the mid-plane, resulting in a substantial drop in  within || <   /2, primarily for the ZNF simulations.To compensate for this resolution-induced drop in , we choose to measure the thermal time based on the range || < 2  .Note that at higher 32/  resolution, even in runs with conduction, the thermal times are similar to the comparable simulations in Paper I for the same combination of A and field configuration.
Finally, the cooling rate in the corona is measured by evaluating the following integral, where the volume integral is taken over the entire corona.Following BQK24, we define the corona as the region || > 2  above the disc mid-plane.In a simulation with radiation transport, the corona is self-consistently defined by the surface above which the optical depth drops to  es ≤ 1.Since our simulations do not include radiation transport, a more appropriate definition for the corona may be the region above   , where   is the height at which the temperature begins to rise rapidly above  0 .Table 1 lists   for all simulations.Since   > 2.5  in strongly cooled (A = 10 3 ) simulations, the cooling rates in the corona (|| > 2  ) may be over-estimated.All simulations display temperature inversions, with a hot corona surrounding a cold disc, though runs without conduction (dashed lines) generally show higher temperatures and larger density enhancements relative to those with conduction (solid and dash-dot lines).Decreased densities and temperatures in the ZNF simulations with conduction relative to runs without conduction are the result of the coronae condensing onto the discs once conduction is activated at 50 orbits.
Ultimately, radiation transport simulations are necessary to evaluate the amount of cooling in the optically thin corona vs. the optically thick disc.However, evaporation or survival of a cold  ≤  0 accretion disc is independent of the definition of the corona.

The absence of disc evaporation
The accretion disc does not evaporate in any of our simulations, independent of the presence of net vertical magnetic flux and the strength of Coulomb coupling in the corona, as parameterized by A.
In what follows, we show quantitatively that the coronal structure and outflow rates do not change in the presence of conduction, consistent with no significant evaporation.Figure 1 illustrates the temperature, density, and magnetic-field structures in our simulations.All of our simulations form temperature inversions, with a hot, diffuse corona 'sandwiching' a colder, thin disc (see panels a and b of Fig. 1).These inversions establish a temperature gradient that channels the conductive heat flux on average toward the disc mid-plane at heights above   .Conductive heat fluxes are directed along magnetic-field lines (panels c and d in Fig. 1), and the vertical heat flux into the disc is suppressed by the toroidal fields formed via the strong Keplerian shear.Simulations with NF, particularly when the field is of moderate strength ( 0 = 10 3 ), exhibit well-ordered magnetic fields, including many nearly vertical field lines threading the corona yet anchored into the disc.Heat fluxes along these field lines are largely uninhibited.
Yet, in spite of the heat flux from the corona, discs with mid-plane temperatures  ≤  0 persist effectively unchanged in all simulations.This feature of our simulations is demonstrated by the horizontally and time averaged temperature and density profiles in Figure 2.
Examining the depletion times in Table 1, we find that conduction has no statistically significant effect on the depletion times, i.e., the mass-outflow rates are not enhanced in any meaningful way by an evaporative outflow.In NF runs,  dep values from simulations with and without conduction agree within the 1 error bars when comparing simulations with the same initial field configuration and the same A.Even though simulations with NF launch strong magnetocentrifugally driven outflows (see BQK24), these outflows are in no way fed by a conduction-induced evaporative flow.
A key goal of this paper is to understand which terms in Equation 8 are responsible for the lack of evaporation in our models, and if we should expect these results to carry over to a realistic, global system.In the subsections that follow, we explore (1) the steepness of the temperature gradient, as measured by |  | ( §4.2), (2) the rate of Coulomb cooling relative to conductive heating and advection by an outflow in the corona ( §4.3), and (3) the role of the magnetic field in geometrically suppressing the heat flux into the disc ( §4.4).

Temperature inversions and conductive cooling
Figure 2 displays the horizontally and temporally averaged temperature and density profiles for all simulations with and without conduction.All ZNF simulations form temperature inversions, with a hot corona forming above || = 3  .Conduction has no effect on the location of   , the height at which the horizontally averaged temperature profile first rises above  0 .This feature indicates a lack of substantial conductive heating in the disc's surface layers.We find  22; orange, solid) and conductive cooling timescale (Equation 22; orange, dashed) to the wind outflow time  wind (red; Equation 24), the measured cooling time  cool (blue; Equation 23), and the orbital timescale  orb (black dashed).For zero-net-flux cases (left column), conductive cooling is important relative to two-temperature cooling above  = 3  , leading to modest condensation of the corona onto the disc.Conductive heating, on the other hand, is always subdominant to two-temperature cooling.In the net flux simulations (right column), the conductive heating and cooling timescales are always longer than the wind-outflow timescales.minimum values of |  |/ ,c in the range of 0.1 < |  |/ ,c < 0.5 for the ZNF simulations and 0.9 < |  |/ ,c < 1.5 for the NF simulations.The effect of conduction is to increase |  |/ ,c , but only by a few tens of per cent relative to runs without conduction.
Rather than heating the disc and moving   inward, conduction in ZNF simulations acts primarily to cool the coronae.The result is a 25 per cent drop in temperature in the weakly cooled, A = 10 ZNF runs.Drops in temperature precipitate decreases in the coronal density as material condenses onto the disc; however, these changes in density structure likely have a negligible observable effect.At || = 2  , all ZNF simulations have approximately the same horizontally averaged densities.Above || = 4  , the average coronal density can drop by as much as a factor of 2 once conduction is turned on.This drop is more substantial for the A = 10 runs, indicating that the density drop is simply due to a decrease in thermal pressure support caused by conductive cooling.
In NF simulations, temperature and density profiles with and without conduction are nearly indistinguishable from one another.The only notable feature, present in both the ZNF and NF simulations, is that in runs without conduction, there is a characteristic peak feature in the temperature profiles above || = 4  and an outward decrease in temperature above this peak.This feature is a boundary effect; the box-height scan in Appendix B of Paper I demonstrates that the peak moves outward with increasing box height.With conduction, the temperature continues to rise monotonically into the boundary and no peak is evident.At these large heights above the disc, the plasma is quite diffuse, and conduction dominates over cooling enough to modify the thermal structure of the corona.

Thermodynamic timescales
Conduction redistributes heat amongst different heights in the accretion disc and corona.We define the conductive heating/cooling timescale as where  con, =  con • ẑ is the conductive heat flux in the -direction.
When / ⟨ con, ⟩   < 0, the above expression refers to the conductive heating timescale, and when the gradient in the heat flux is positive,  cond is the conductive cooling timescale.The conductive heating/ cooling timescale can then be compared to the cooling time, and the wind-outflow timescale, i.e. the timescale over which the outflow can transport plasma over a height 2  , Figure 3 displays a comparison of these different timescales.In the ZNF simulations, conductive heating (solid orange line in Figure 3) acts only over a very narrow range of , less than half a mid-plane scale height above   .Within this range, conductive heating is slow compared to cooling (blue curves), with the conductive heating timescale exceeding the cooling timescale by an order of magnitude in the strongly cooled, A = 10 3 cases, and by a factor of a few in the weakly cooled, A = 10 runs.
Conductive heating occurs over a much broader range of heights in the NF simulations than in the ZNF simulations.The range of conductive heating extends from   to ≈(  +   ) for all but the NF, A = 10 3 ,  0 = 10 3 run.As in the ZNF simulations, cooling outpaces conductive heating at all .Unlike the ZNF simulations, though, the wind-outflow timescale is shorter than the conductive heating timescale in all NF runs, with the exception of a very narrow range near  =   = 1.5  in the NF, A = 10,  0 = 10 4 simulation.Hot plasma in the corona is transported upwards more rapidly than conduction can channel thermal energy down towards the disc.
We can extrapolate these results to a global system by comparing the best-case scenario evaporation timescale,  evap , to the inflow (viscous) timescale of the disc,  visc ≡ ( 0 /  ) 2  thm .Thinner discs slow the inward flow of plasma.In the absence of two-temperature cooling, we can compute a critical disc thickness for conduction to evaporate the disc within the inflow time.The critical disc thicknesses to allow  evap <  visc are   / 0 ≈ (1 − 7) × 10 −2 for the NF simulations, and   / 0 ≈ 0.05 − 0.3 for the ZNF simulations.Particularly for the NF simulations, these discs are quite thin; however, such a thin disc may be realizable.This simple argument demonstrates that the lack of evaporation is not due to a negligible heat flux into the disc.Instead, Coulomb cooling radiates away the heat channeled toward the disc, preventing evaporation.
Notably, while the conductive heating/cooling timescales are long compared to the two-temperature cooling time in all NF simulations, for weakly cooled (A = 10) ZNF simulations, the conductive cooling timescale is actually shorter than the cooling timescale imposed by Coulomb collisions above a height  ≥ 3  .Efficient conductive cooling in the diffuse coronae of the ZNF simulations thus results in a collapse of the corona onto the disc in the form of a condensing inflow.We study the behaviour of this inflow further in §4.5.

Geometric suppression of the heat flux
The vertical conductive heat fluxes are suppressed by the geometry of the predominantly toroidal magnetic field.To quantify this geometric suppression, we introduce the suppression factor   , and we define the horizontal average of the suppression factor as Profiles of the conductive suppression factor ⟨  s ⟩   for all simulations, averaged from 30-100 orbits for the NF simulations, and averaged from 70-100 orbits in the ZNF cases.The suppression factor in regions where conduction is heating the plasma (solid) and cooling the plasma (dash-dotted) are shown separately, and a comparison to a naïve estimate for the suppression factor (viz., the horizontal average of | ẑ • b|), is shown by the dotted lines.Where available, we use the high-resolution run for the profiles, but show the same result for the lower-resolution analogs using more transparent curves.We only plot the suppression factor in regions where  ≥   , as there is no consistent temperature gradient below this height.lower in the ZNF simulations than those in the NF simulations.This trend is a result of the fact that fields with larger initial NF are more rigid, allowing them to remain more inclined with respect to the disc mid-plane and thus maintain a stronger coupling between the upper, hotter regions of the corona and cooler regions just above   .Generally, the suppression factor is in the range 0.1 < ⟨   ⟩   < 0.2 for the ZNF simulations, while a much larger range is spanned by the NF simulations: 0.07 < ⟨   ⟩   < 0.6.Stronger cooling (larger A) results in less suppression of the heat flux.This is because stronger cooling inhibits winds and evaporation of material into the corona.Thus, less mass is loaded on field lines, the field can maintain larger inclination angles, and heat fluxes can more easily pass among layers of the stratified flow.
Interestingly, the magnitude of the suppression factor is resolution-dependent, and simulations performed at higher resolution exhibit less suppression of the heat flux.At high resolution, we find suppression factors similar to those used by CN22.Lower resolution runs, particularly the NF, A = 10,  0 = 10 4 simulation at 16 cells/  resolution, display suppression factors well below these values.
Naïvely using the field geometry as a proxy for the suppression factor by horizontally averaging the quantity | ẑ • b|, i.e. the pitch angle cosine, shown as dotted lines in Figure 4, also underestimates the heat flux by an order unity factor of ≈2−5.The pitch angle cosine accounts for the geometry of the field, but not the direction of the heat flux, which is set by the temperature gradient.Because the flow is highly turbulent, the heat flux is not always directed toward the disc.This feature of the solutions is visible in panel (d) of Fig. 1.There, we see that the same field line can have neighbouring red regions (where the heat flux is toward the disc) and blue regions (where the heat flux is away from the disc), with a sharp transition between them, caused by a change in the direction of the temperature gradient along the field line.Failing to account for the direction of the gradient results in an underestimate of the suppression factor, and subsequently, an underestimate of the field-aligned heat flux into the disc.

Condensing inflows with zero net flux
By forming a condensing inflow, conductive cooling allows plasma cooling out of the corona to feed the thin disc.Figure 5 demonstrates this property of our ZNF solutions.We define the hot phase as plasma with a temperature  ≥ 5 0 1 and the mass in the corona as the total mass above || =   .Without conduction turned on, the hot phase comprises ≈10 to 25 per cent of the coronal plasma by mass for the A = 10 runs and ≈5 per cent in the A = 10 3 cases.Activating conduction at  = 50 orbits results in a sudden decrease in the 1 Our results are largely independent of our choice for the minimum temperature of the 'hot phase,' at least for the runs with A = 10.Using a minimum hot phase temperature of 3 0 or 1.5 0 results in a 60 per cent decrease and 50 per cent decrease in the hot phase mass, respectively, once conduction is turned on in the A = 10 ZNF simulations.For A = 10 3 , turning on conduction does not affect the hot phase mass when the minimum hot phase temperature is 3 0 or 1.5 0 .and / 0 for all simulations, normalized such that the integrals under the curves equal 1, with black lines corresponding to runs without conduction, and red lines corresponding to simulations with conduction.Conduction reduces the spread in ion temperature within these multiphase coronae; however, the spread in density is largely unaffected, implying that conduction should not affect the optical depth and 'clumpiness' of two-temperature coronae.
hot phase mass as the corona cools via thermal conduction and a condensing inflow forms.The hot phase mass drops by a factor of ≈3 for weak cooling and ≈2 for strong cooling, and coronal plasma can rain onto the disc.This 'coronal rain' carries very little mass-2 to 3 per cent of the disc's mass in the weakly cooled, A = 10 simulations, and only 0.1 per cent of the disc's mass in the more strongly cooled runs.Thus, mass accretion rates should be unaffected by this condensing inflow.However, such an inflow may rain optically thick material into the optically thin surface layers of the disc, if this material can cool and condense during its infall (see §6.4).

Smoothing of multiphase temperature structure
Figure 6 shows log-density () and temperature () distributions measured at ||/  = 3 for the ZNF  0 = 10, A = 10 simulations without conduction (i.e., from 30 to 50 orbits; top left) and with conduction active (from 70 to 100 orbits, long enough after conduction is turned on at 50 orbits for the solution to reach a steady state; bottom left).Beside these distributions, we show similar distributions for the NF,  0 = 10 4 , A = 10 simulations without (top right) and with (bottom right) conduction, measured at ||/  = 5 from 30 to 100 orbits.As in BQK24, we find that the solutions  26).(a) Cold plasma mass (defined as the mass of material with temperature  ≤ 1.5 0 ) as a fraction of the initial cold plasma mass,  cold,0 , as a function of  and the conductive heat flux averaged above  = 5  (panel b).Temperature profiles (c) and density profiles (d) for all toy model simulations, with darker colours denoting lower values of  and the black, dashed lines indicating the initial profiles.(e) Conductive heat flux profiles for all toy model simulations, measured in terms of the local saturated heat flux,  sat ≡ 0.6 3 s .We define the constant  sat,0 ≡ 0.6 0  3 s,0 .Models do not evaporate, i.e., the cold plasma mass remains  cold > 0 for the duration of the simulation, unless  > 10, or equivalently, the heat flux into the disc is larger than 10 3  sat,0 .These heat fluxes are unphysically large, far above the free-streaming limit (see panel e).
without conduction exhibit multiphase structure, with broad density and temperature distributions that form along curves of constant pressure (shown as light blue dotted curves in Fig. 6).
In the ZNF sims, this characteristic structure vanishes once conduction is activated at 50 orbits.Temperature distributions become sharply peaked, transitioning from a spread of ≈11 0 at || = 3  in the ZNF simulations to only ≈3-4 0 at the same height.The spreads in density remain virtually unchanged in both ZNF and NF simulations, although there may be a slight broadening in log-density space in addition to a shift of the distributions with conduction toward densities lower than those characteristic of simulations without conduction.While a condensing inflow in the ZNF simulations with conduction results in a substantial drop in density, the NF simulations remain outflowing at all || > 2  , and the mean density is unchanged by the inclusion of conduction.

A TOY MODEL FOR EVAPORATION
Conductive evaporation does not occur in any of our simulations.The geometric suppression of the heat flux plays some role in this outcome, decreasing the heat flux into the disc by an order unity factor.However, as we assessed in §4.3 and Figure 3, the dominant process opposing evaporation is the strong two-temperature cooling in the surface layers of the disc.This conclusion leads us to the simple question: 'How large must the heat flux into the disc be to overcome two-temperature cooling?' To answer this question, we design a toy model for the corona-disc system.Our goal is to set up an atmosphere in hydrostatic equilibrium (to control for the effects of the outflow) with a purely vertical magnetic field (to ensure that the heat flux is unsuppressed, with   = 1), where a hot corona sandwiches a cold disc (to ensure that the heat flux is toward the disc).We then turn on conduction and cooling from this initial set up, and observe whether or not an evaporative flow forms.Instead of using the free-streaming form of the heat flux (4), we use a Spitzer-esque heat flux of the form where  is a free parameter.Note that we intentionally use a constant thermal diffusivity set by the mid-plane density and temperature, which can lead to heat fluxes well in excess of the free-streaming value.By artificially varying  and therefore the magnitude of the heat flux into the disc, we can assess how large the heat flux into the disc must be to form an evaporative flow.The choice of the heat flux is discussed further in §5.1.

Problem set-up
We restrict our calculation to one dimension, the -direction, which extends from −8 < /  < 8, although the results do not change for a domain with a maximum height of 6  .Such a tall domain minimizes the effect of the boundary condition on the thermal structure of the solutions near the mid-plane.The resolution chosen for this domain is the same as that used in the high-resolution three-dimensional simulations, 32 grid cells/  .Just as in the shearing-box simulations, we measure distances in terms of the initial mid-plane scale height   ; however, these toy model simulations do not include the effects of rotation.Thus, we measure time in units of the free-fall time,  ff ≡   / s,0 .
Motivated by the peak temperature of ≈9 0 achieved in the ZNF A = 10 simulations without conduction, we initialize an atmosphere with a maximum temperature of  max = 10 0 above || ≈ 2  (in the corona) and a minimum temperature of  0 below this height (in the disc).We initialize the atmosphere with a temperature profile, This profile is shown by the dashed black line in panel (c) of Fig. 7.The corresponding density profile is computed by solving the equation of hydrostatic equilibrium for the disc's local gravitational potential, Φ = ( s,0 /  ) 2  2 .As a result, we achieve a similar initial density profile within the disc to that exhibited in the steady state of the ZNF and weak NF simulations within the same region of || < 2  (Fig. 2).Above || = 2  , the density in the toy-model initial condition drops off much less steeply compared to the density profiles in the three-dimensional simulations.Because our atmosphere is initially hydrostatic and not outflowing, the initial density in the corona is higher than that of the coronae of the ZNF and weak NF simulations.The upper and lower boundaries are handled in the same way as was done for the fiducial simulations, i.e., the conductive heat flux  con is copied into ghost zones with zero gradient.Unlike the fiducial simulations, however, we pin the outer boundary temperature to 10 0 .This condition ensures that there is a continuous heat flux from the outer boundary into the domain and toward the disc.We choose a high value of  max = 10 3   Ω to ensure that, even with artificially high heat fluxes well in excess of the free-streaming limit, the solutions remain accurate.The heat flux into the disc is removed through our cooling function ( 9), where we choose A = 10 to match our most weakly cooled three-dimensional models.We choose an extremely low floor of  floor / 0 = 10 −10 , although the solutions never come within an order of magnitude of this floor density.
While we tried to implement the free-streaming heat flux in these toy model simulations, we found that this form of the heat flux does not achieve a steady state sufficient to allow for a measurement of evaporation.The atmosphere collapses onto the disc, causing a shock wave, which reverses the direction of the heat flux such that conduction channels heat out of the disc.The reason for this behaviour is that, unlike the heat flux used in our toy model (Equation 26), the free-streaming heat flux depends on density .Because of the large density gradient (relative to the temperature gradient), in the corona-disc transition zone, the heat flux nearest the disc is much higher than that higher up in the corona.As a result, the bottom of the corona cools first and begins to collapse onto the disc before sufficient thermal energy from higher in the corona can arrive to offset this cooling.Thus, no steady state with the free-streaming heat flux is possible, and using the free-streaming heat flux would prevent us from accurately measuring the magnitude of the heat flux required to induce evaporation.
We run the simulation until  ≈ 89 ff .By this time, the solutions have relaxed into a steady state, and the fate of the disc, whether the disc evaporates or remains intact, is sealed for all time.

Toy model results
Figure 7 shows the results of these calculations for 16 values of  in the range 10 −2 ≤  ≤ 10 3 .We measure the cold gas mass  cold , defined as the mass of material with temperature  ≤ 1.5 0 , relative to the initial cold gas mass  cold,0 measured at  = 0.As is clear from panels (a) and (b) of Fig. 7, the cold disc is completely destroyed by  ≈ 89 ff for any heat flux corresponding to  > 10.
Panel (b) of Fig. 7 shows the final cold gas mass at  ≈ 89 ff as a function of the heat flux  con averaged above  > 5  .We measure the heat flux relative to the peak saturated value achievable for a plasma with temperature  0 and density  0 , i.e.,  sat,0 ≡ 0.6 0  3 s,0 .Profiles of heat flux relative to the local saturated value,  sat ≡ 0.6 3 s , are shown in panel (f).The critical heat flux for evaporation is  crit ≈ 10 3  sat,0 -far larger than any heat flux that can be physically attained in the system.Similarly, models that evaporate display heat fluxes of | con, |/ sat ≳ 10 above || = 2  .Two-temperature cooling is so strong in the disc's surface layers that only a completely unphysical heat flux can induce evaporation.This analysis underscores the primary conclusion of this paper: ion heat fluxes from the corona into the disc are insufficient to evaporate the disc, even when the heat flux achieves its maximum value at the free-streaming limit.

Comparison to previous work
Our conclusions are similar to CN22, in that we find that conduction in the inner regions of XRB accretion discs is more likely to lead to a condensing inflow from the corona onto the disc than to result in disc evaporation.Yet, we arrive at our conclusions by examining very different physical mechanisms: free-streaming ion thermal conduction vs. Spitzer electron conduction, and two-temperature ion cooling via Coulomb collisions rather than Bremsstrahlung cooling of a one-temperature plasma.

Limitations of our model
Physical effects not captured in our simple two-temperature model may be important for understanding the role of conduction in disc evaporation.In the bulk of the disc where the plasma becomes one-temperature, electron thermal conduction dominates the heat flux.Electron conduction is likely to operate in the Spitzer regime (as in CN22), rather than in the free-streaming regime.Further, radiation transport may substantially modify the temperature structure of the disc's surface layers, potentially forming a steep temperature gradient approaching the mid-plane.
In addition, the form of our two-temperature cooling function (9) implies that cooling becomes stronger with increasing temperature.Bremsstrahlung cooling behaves similarly, with the cooling rate scaling ∝ 1/2  .However, if the disc is radiation pressure dominated, gas temperatures may be low enough that substantial line cooling takes place near the mid-plane.At  ≈ 10 5−6 K, the cooling rate decreases with increasing temperature, and conductive heating would lead to even weaker cooling.Such thermally unstable heating may play some role in the bulk of the disc.Indeed, this is the regime relevant for the solar corona and for white dwarf coronae, where conductive evaporation is traditionally studied.We note that such low temperatures are difficult to achieve in XRB discs, whose inner edges are hot enough to emit soft, thermal X-rays.The 10 6 K branch of the cooling curve is more likely to be approached at larger radii from the black hole in XRB discs, or in much cooler, much more radiation dominated flows, such as those fueling AGN.

Soft-to-hard state transitions in XRBs
Our primary motivation for exploring ion thermal conduction was to find the long-sought physical mechanism that causes disc evaporation, truncation, and therefore, state transitions in XRB accretion flows.Based on our calculations, ion thermal conduction cannot serve this role.Heat fluxes orders of magnitude larger than what can be mustered by free-streaming ion conduction are necessary to force evaporation in the presence of two-temperature cooling.
The presence and magnitude of NF, rather than the inclusion of thermal conduction, has the most substantial effect on our two-temperature disc-corona systems.Simulations with zero net vertical magnetic flux (ZNF) do not evaporate, but instead, their coronae cool via conduction and form condensing inflows onto the disc.Similarly, models with NF launch magnetocentrifugal outflows.Such outflows may act to deplete the disc, reducing its surface density until Coulomb collisions become sufficiently infrequent that the flow becomes two-temperature.The runaway ion heating that ensues might evaporate the disc.Similarly, magnetic pressure support afforded by NF fields allows for lower disc densities and optical depths compared to those in gas-pressure-dominated and radiation-pressure-dominated flows (Mishra et al. 2022;Huang et al. 2023).These low densities, combined with inefficient cooling at these reduced optical depths, encourage runaway ion heating and evaporation.Thus, the introduction of NF, through providing magnetic pressure support and through launching winds that deplete the disc, might be the agent of disc evaporation and truncation in XRBs (Ferreira 1997;Ferreira et al. 2006;Begelman & Armitage 2014;Liska et al. 2023).

Hard-to-soft transitions and iron lines
Conduction acts as an efficient cooling mechanism in coronae threaded by ZNF, producing a condensing inflow onto the disc.In this way, just as the introduction of NF from larger scales may prompt the evaporation of the disc through wind depletion, the removal of NF and the subsequent cessation of such a wind may lead to sudden condensation of the inner hot RIAF into a disc.Indeed, our simulations show that the introduction of conduction into the two-temperature ZNF simulations can decrease the mass of plasma in the hot phase of the corona by as much as 80 per cent in the weakly cooled, A = 10 simulations (Fig. 5).If this physics of condensation of the corona onto the disc operates in real systems, it could contribute to a hard-to-soft state transition and the re-formation of Fe lines in the inner regions of XRBs.Thus, the impact of ion thermal conduction on the disc-corona system is intimately tied to whether or not NF fields thread the accretion flow.
While the presence and strength of the Fe line signal may be affected by conductive cooling and condensation in the corona, the hard X-ray power-law emission produced by coronae is likely unaffected by conduction.Even though conduction significantly decreases the spread in temperature in the coronae of our models (Fig. 6), the spread in density is virtually unchanged.Since we are evolving the ion temperature rather than the lepton temperature, the temperature spread would be unlikely to affect the cut-off in the X-ray spectrum.Similarly, the negligible effect on the density spread implies that the optical depth to electron scattering, and therefore the power-law index measured from observations, is unlikely to be changed by free-streaming ion conduction.

SUMMARY AND CONCLUSION
We have implemented field-aligned thermal conduction at the saturated, free-streaming limit, into local simulations of vertically stratified accretion flows.By treating the ions in the corona as a single MHD fluid subject to cooling via Coulomb collisions with radiatively efficient electrons, we capture temperature inversions reminiscent of a hot corona surrounding a colder disc, which then direct a conductive heat flux from the corona into the disc.Through the application of our two-temperature model with conduction to a large suite of three-dimensional stratified shearing-box simulations and a one-dimensional, highly idealized, toy model, we have shown that ion thermal conduction is unable to evaporate thin accretion discs into RIAFs, as has previously been proposed (Spruit & Deufel 2002).Instead, two-temperature cooling in the surface layers of discs removes energy from the ions, radiating this thermal energy away before any substantial conductive heating can impact the bulk of the dense disc.
Our main results are as follows: (i) All of our two-temperature models with thermal conduction form temperature inversions (Figure 2, left panels), with a hotter corona surrounding a colder disc.In NF simulations, conduction has no effect on the ion temperatures in the corona.For ZNF simulations, conduction acts to cool coronae, decreasing their maximum horizontally and temporally averaged temperatures by about 25 per cent in the weakly cooled, A = 10 simulations, independent of the initial strength of the ZNF fields.
(ii) The temperature inversions allow for a net conductive heat flux from the hot corona into the cold disc.This heat flux is not negligible.In discs thinner than   / 0 ≈ a few ×10 −2 , conductive evaporation would destroy the discs more rapidly than inflow could empty them, if cooling were ignored.
(iii) Cooling through Coulomb collisions between hot ions and rapidly cooling leptons far outpaces the heating delivered by thermal conduction into the disc's surface layers.The two-temperature cooling timescale is shorter than the conductive heating timescale throughout the coronae of all of our simulations ( §4.3 and Fig. 3), independent of the presence and strength of NF and the strength of Coulomb cooling (as parameterized through the Coulomb coupling parameter, A) in the corona.
(iv) The dominantly toroidal magnetic fields formed through MRI turbulence and Keplerian shear in the stratified shearing box act to geometrically suppress the heat flux from the corona into the disc.We find suppression factors of ⟨   ⟩   ≈ 0.1 − 0.6 (Fig. 4).The heat flux is less suppressed in moderate NF simulations, which exhibit field lines that are more rigid and more vertically aligned, on average compared to ZNF simulations.Higher-resolution simulations show less geometric suppression of the heat flux, and estimates based on post-processing generally under-estimate the heat flux into the disc.
(v) Using a one-dimensional toy model of the disc-corona system, where the 'corona' above two mid-plane scale heights is 10× hotter than the disc and the system is initially in hydrostatic equilibrium, we showed that the heat flux from the corona into the disc must be unphysically large to overcome two-temperature cooling in the disc and induce evaporation ( §5 and Fig. 7).
(vi) Rather than forming outflows, ZNF simulations with conduction form condensing inflows, which can remove as much as 80 per cent of the hot phase plasma mass from the coronae ( §4.5 and Fig. 5).However, because of how little mass is contained in the corona, this condensing inflow is unlikely to affect the accretion rate or luminosity of the system.
(vii) Thermal conduction smooths out multiphase temperature structure, shrinking the spread in temperature distributions and eliminating the characteristic nearly isobaric density-temperature structure at fixed height || that was observed in the two-temperature simulations of Paper I. The spread in densities, and consequently, observational signatures that can be probed by the power-law index of the X-rays released by the corona, are unaffected by the inclusion of thermal conduction.
Our work indicates that ion thermal conduction cannot provide the physical mechanism that enables soft-to-hard state transitions in XRBs.The introduction of NF changes the vertical structure of two-temperature coronae from a condensing inflow cooled via thermal conduction into a magnetocentrifugal outflow.Similarly, removal of NF results in condensation of a hot corona and the feeding of a thin, optically thick accretion disc near the ISCO.The presence of significant magnetic support also decreases disc densities, leading to less efficient cooling and promoting the formation of a RIAF (e.g., Mishra et al. 2022 andHuang et al. 2023).We thus speculate that a change in the amount of NF threading an accretion flow, either through in-situ generation of poloidal flux via a dynamo or transport of this flux from larger scales, is responsible for the state transitions observed ubiquitously in XRBs, consistent with works by Begelman & Armitage (2014) and recently, Liska et al. (2023).velocities, are more likely to introduce advective errors that increase the size of the time-derivative term in Equation 11.
Figure B1 shows profiles of the horizontally averaged temperature, density, turbulent  parameter, and vertical conductive heat flux for these three simulations.The solid lines denote the results computed with  max = 30   Ω (the fiducial  max used throughout the paper) while dashed lines show the results after increasing  max to 10 2   Ω at 50 orbits.The temperature, density,  and vertical heat flux profiles are relatively similar between simulation runs presented in the main paper and runs with increased  max , indicating that our chosen  max is sufficiently large that the heat flux is reasonably close to the saturated value.Just as in the profiles with  max = 30   Ω, increasing  max does not result in evaporation of the cold disc, in agreement with the conclusions presented throughout the main paper.
Similarly, in Figure B2, we demonstrate that our conclusion that thermal conduction causes a condensing inflow onto the disc in simulations without NF (see Fig. 5) does not depend on the chosen value of  max .Finally, we note that if  max is not large enough, such that Equation 12 is not satisfied, the two-moment method is more likely to under-estimate the total heat flux rather than over-estimate this flux.Thus, we checked that the profiles of the horizontally averaged conductive heat flux produced by the simulations, i.e. ⟨ con, ⟩   were in agreement with the expectations from the free-streaming value: ⟨ sat ( con, /| con |)⟩   .To compute these profiles, we restricted our computations to hot zones in the corona where  > 1.5 0 , since in the disc midplane, a uniform temperature would imply no heat flux while  sat would be 0.6 0  3 ,0 .We find that the post-processed and simulated heat fluxes agree to within 5 per cent error for all simulations, with smaller errors corresponding to higher resolution and higher  max runs.
The constraint on  max implied by Equation 12 is more difficult to satisfy if | b • ∇ ln  | is small.The reason is that the conductivity  in Equation 13 is then larger and so a larger  max is necessary for the heat flux to equal the equilibrium (saturated) value in the two-moment method.The dominantly toroidal nature of the field lines in our local disc patches implies that temperature will be relatively uniform horizontally (e.g., the narrow temperature distributions at fixed height in the ZNF simulations in Fig. 6).Thus with our chosen values of  max we somewhat underestimate the horizontal heat fluxes relative to the correct saturated value.However, when field lines possess a significant vertical component, the field-aligned temperature gradient samples a large variation in temperature,  is small enough to be controlled by the chosen  max , and the heat fluxes are accurate such that | con | ≈  sat .These zones with substantial vertical components to the field lines dominate the vertical transport of heat and thus mediate conductive heating, conductive cooling, and evaporation (or lack thereof) in our simulations.Restarting the simulation at  = 50 orbits with a higher value of  max has a negligible effect on the time-averaged amount of hot gas in the corona (and thus the amount of material rained onto the disc).The solutions diverge from one another over time due to chaos intrinsic to the turbulent system.

Figure 2 .
Figure2.Left: Horizontally averaged temperature () profiles from all simulations.Profiles for net-flux simulations with and without conduction are averaged over 30-100 orbits, while zero-net-flux simulations are averaged over two different time frames: 30-50 orbits for 'No Conduction' runs and 70-100 orbits for 'Conduction' runs.Right: Horizontally averaged density () profiles for all simulations over the same time frames.Profiles from both the standard-resolution (16 grid zones per   ) and high-resolution (32 grid zones per   ) simulations are shown.All simulations display temperature inversions, with a hot corona surrounding a cold disc, though runs without conduction (dashed lines) generally show higher temperatures and larger density enhancements relative to those with conduction (solid and dash-dot lines).Decreased densities and temperatures in the ZNF simulations with conduction relative to runs without conduction are the result of the coronae condensing onto the discs once conduction is activated at 50 orbits.

Figure 3 .
Figure3.Timescale profiles for all simulations with thermal conduction.We compare the conductive heating timescale (Equation22; orange, solid) and conductive cooling timescale (Equation22; orange, dashed) to the wind outflow time  wind (red; Equation24), the measured cooling time  cool (blue; Equation23), and the orbital timescale  orb (black dashed).For zero-net-flux cases (left column), conductive cooling is important relative to two-temperature cooling above  = 3  , leading to modest condensation of the corona onto the disc.Conductive heating, on the other hand, is always subdominant to two-temperature cooling.In the net flux simulations (right column), the conductive heating and cooling timescales are always longer than the wind-outflow timescales.

)
Figure 4. Profiles of the conductive suppression factor ⟨  s ⟩   for all simulations, averaged from 30-100 orbits for the NF simulations, and averaged from 70-100 orbits in the ZNF cases.The suppression factor in regions where conduction is heating the plasma (solid) and cooling the plasma (dash-dotted) are shown separately, and a comparison to a naïve estimate for the suppression factor (viz., the horizontal average of | ẑ • b|), is shown by the dotted lines.Where available, we use the high-resolution run for the profiles, but show the same result for the lower-resolution analogs using more transparent curves.We only plot the suppression factor in regions where  ≥   , as there is no consistent temperature gradient below this height.

Figure 5 .
Figure5.Evolution of hot phase ( ≥ 5 0 ) in the corona as a function of time for all ZNF simulations.The period when thermal conduction is turned off at the minimum value of  =  min is indicated by the dashed portions of the curves while the period during which conduction is active at the free-streaming value is shown by the solid portion of the curves.We define the initial coronal mass M cor,0 as the total mass in the region |  | ≥   from  = 30-50 orbits, when conduction is not active.For the A = 10  0 = 10 simulation, both the high resolution run (darker) and the fiducial resolution run (lighter/ more transparent) are shown.Conductive cooling produces a condensing flow from the corona onto the disc, reducing the mass of plasma in the hot phase by as much as 80 per cent when cooling is weak (A = 10) and 50 per cent for much stronger two-temperature cooling (A = 10 3 ), independent of the initial strength of the ZNF field.

Figure 6 .
Figure 6.Density-temperature distributions for all plasma in high resolution  0 = 10 A = 10 ZNF simulations at times 30-50 orbits ('No Conduction'; top left) and 70-100 orbits ('Conduction'; bottom left) at height |  | = 3  (left column).We show similar density-temperature distributions for all of the plasma in high resolution NF  0 = 10 4 A = 10 simulations at times 30-100 orbits for the simulation without conduction ('No Conduction'; top right) and the run with conduction ('Conduction'; bottom right) at height |  | = 5  (right column).White curves over the distributions represent lines of constant entropy, while light blue curves represent lines of constant pressure.We show one-dimensional probability density functions (PDFs) in log (/ 0 )and / 0 for all simulations, normalized such that the integrals under the curves equal 1, with black lines corresponding to runs without conduction, and red lines corresponding to simulations with conduction.Conduction reduces the spread in ion temperature within these multiphase coronae; however, the spread in density is largely unaffected, implying that conduction should not affect the optical depth and 'clumpiness' of two-temperature coronae.

Figure 7 .
Figure 7. Toy results at time ≈ 89  / s,0 into the simulation.The colours denote the value of the conduction amplification factor  used in the toy model heat flux (see Equation26).(a) Cold plasma mass (defined as the mass of material with temperature  ≤ 1.5 0 ) as a fraction of the initial cold plasma mass,  cold,0 , as a function of  and the conductive heat flux averaged above  = 5  (panel b).Temperature profiles (c) and density profiles (d) for all toy model simulations, with darker colours denoting lower values of  and the black, dashed lines indicating the initial profiles.(e) Conductive heat flux profiles for all toy model simulations, measured in terms of the local saturated heat flux,  sat ≡ 0.6 3 s .We define the constant  sat,0 ≡ 0.6 0  3 s,0 .Models do not evaporate, i.e., the cold plasma mass remains  cold > 0 for the duration of the simulation, unless  > 10, or equivalently, the heat flux into the disc is larger than 10 3  sat,0 .These heat fluxes are unphysically large, far above the free-streaming limit (see panel e).

Figure A1 .Figure B1 .
Figure A1.Dispersion relation for a sound wave subject to thermal conduction in the Spitzer regime (A1).Solid lines show the numerical solution to the dispersion relation (A3), while points indicate measurements based on a series of linear sound wave tests (see text).The real frequency   is indicated with the color blue, while the imaginary frequency/ damping rate −   (the minus sign indicates damping of waves) is indicated by the color orange.

Figure B2 .
Figure B2.Evolution of hot phase ( > 5 0 ) in the corona as a function of time for the high resolution ZNF A = 10  0 = 10 simulation, similar to Figure5, for different values of  max .Restarting the simulation at  = 50 orbits with a higher value of  max has a negligible effect on the time-averaged amount of hot gas in the corona (and thus the amount of material rained onto the disc).The solutions diverge from one another over time due to chaos intrinsic to the turbulent system.