A 3D radiation-hydrodynamic AGB binary model
DDraft version March 2, 2020
Typeset using L A TEX twocolumn style in AASTeX63
A 3D radiation-hydrodynamic AGB binary model
Zhuo Chen,
1, 2
Natalia Ivanova, and Jonathan Carroll-Nellenback Department of Physics, University of Alberta, Edmonton, AB T6G 2E1, Canada CITA National Fellow Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627, USA
ABSTRACTThe origin of chemically peculiar stars and non-zero eccentricity in evolved close binaries have beenlong-standing problems in binary stellar evolution. Answers to these questions may trace back toan intense mass transfer during AGB binary phase. In this work, we use
AstroBEAR to solve the3D radiation-hydrodynamic equations and calculate the mass transfer rate in asymptotic-giant-branch(AGB) binaries that undergo the wind-Roche-lobe-overflow or Bondi-Hoyle-Lyttleton (BHL) accretion.
MESA produces the density and temperature of the boundary condition of the AGB star. To improvethe resolution of the dynamics of a circumbinary disk, we implement an azimuthal angle-dependent 3Dradiation transfer. We consider optically thin cooling and obtain the number density of the coolantsby solving Saha equations. One of the goals of this work is to illustrate the transition from the wind-Roche-lobe-overflow to BHL accretion. Both circumbinary disks and spiral structure outflows canappear in the simulations. Circumbinary disks may form when the optical thickness in the equatorialregion increases. The increase of the optical thickness is due to the deflected wind. The resultingmass transfer efficiency in our models is up to a factor of 8 times higher than what the standard BHLaccretion scenario predicts, and the outflow gains up to 91% of its initial angular momentum whenit reaches 1.3 binary separations. Consequently, some AGB binaries may undergo orbit shrinkage,and some will expand. The high mass transfer efficiency is closely related to the presence of thecircumbinary disks.
Keywords: binaries: symbiotic – methods: numerical – stars: AGB and post-AGB – stars: winds,outflows INTRODUCTIONAsymptotic-giant-branch (AGB) stars have a signif-icantly larger size ( ∼ Corresponding author: Zhuo [email protected] author: Natalia [email protected] −
20 km · s − (H¨ofner & Olofsson 2018), and a compan-ion star may capture the wind with its gravity. In thecase that there is a MS star close to an AGB star, a sub-stantial fraction of the mass-loss may be accreted ontothe MS companion (Chen et al. 2017; Saladino et al.2018, 2019). As a result, the metallicity of the com-panion may change. Such early-stage low-mass starsbecome chemically peculiar, and their future evolutionwill be strongly affected. Carbon-enhanced-metal-poor(CEMP) stars (Beers & Christlieb 2005; Abate et al.2013, 2015), Barium stars (Bidelman & Keenan 1951;Escorza et al. 2019), CH stars (Keenan 1942; McClure &Woodsworth 1990) and dwarf carbon stars (Dahn et al.1977; Roulston et al. 2019) are common examples ofthe chemically peculiar stars. Their existence could bethe evidence of the mass transfer during the previousAGB binary phase. The binarity of CH stars and CEMPstars has been studied (McClure & Woodsworth 1990; a r X i v : . [ a s t r o - ph . S R ] F e b Starkenburg et al. 2014; Jorissen et al. 2016), confirm-ing that many of them have companions. A number ofrecent studies show that the eccentricity of some of theaforementioned chemically peculiar stars may be large(Hansen et al. 2016; Jorissen et al. 2016; Van der Swael-men et al. 2017; Oomen et al. 2018; Jorissen et al. 2019),and their orbital periods ranging from hundreds to thou-sands of days. The non-zero eccentricity in these closebinary stars indicates some intense interactions that canpump the eccentricity may happen during their AGB bi-nary phases. A strong correlation between a circumstel-lar disk and binarity has also been established in galac-tic RV Tauri stars (Manick et al. 2017). Furthermore,many RV Tauri stars show a lack of refractory elements,which is called ’depletion’ (Giridhar et al. 1994; VanWinckel et al. 1998). Some researches suggest that thereaccretion of gas from a circumstellar disk around thepost-AGB star (Gezer et al. 2019; Oomen et al. 2019)may induce the ’depletion’. Besides the ’smoking gun’evidences, observations also reveal that dusty circumbi-nary disks exist in binary systems with evolved stars(Kervella et al. 2015; Hillen et al. 2016; Homan et al.2017; Ertel et al. 2019). UV excess of some AGB starsalso suggests that there could be accreting MS compan-ions near them (Sahai et al. 2008; Ortiz & Guerrero2016).In close AGB binaries, the mass transfer efficiency al-most determines the intensity of the interaction betweenthe two stars. However, it has been a difficult problemto quantify the mass transfer efficiency observationallybecause of the uncertainty in the measurement of themass-loss rate and accretion rate. Fortunately, with therapidly increasing computational power, it is now possi-ble to quantify the mass transfer efficiency in AGB bi-naries by carrying out 3D global numerical calculations.Figure 1 shows a schematic and global picture of an AGBbinary undergoing wind-Roche-lobe-overflow (WRLOF)as described in Podsiadlowski & Mohamed (2007). Themain difference between the WRLOF and Roche-lobeoverflow (RLOF) is that the AGB star is not fillingits Roche-lobe to the full. Therefore, the binary sep-aration of an AGB binary undergoing WRLOF shouldbe larger than the binary separation if it is undergoingRLOF. On the other hand, if the binary separation isso significant that not only the physical size of the AGBstar but also its dust formation region is much smallerthan the Roche lobe of the AGB star, the mass trans-fer mechanism approaches the Bondi-Hoyle-Lyttleton(BHL) accretion limit (Hoyle & Lyttleton 1939; Bondi& Hoyle 1944; Edgar 2004). In that scenario, the ra-diation driven AGB wind becomes supersonic before itcrosses the Roche lobe of the AGB star. The two stars y / d L1-potentialL2-potentialL3-potentialdust formationionization
Figure 1.
A schematic picture of the WRLOF. The coor-dinate is scaled to the binary separation d . The contours ofthe L1, L2, and L3 Lagrangian potential are drawn in red,green, and blue dashed lines. The AGB star is in yellowcolor, which is filling a fraction of its Roche lobe. Its dustforming region is represented by the gray ring. The compan-ion star may have an ionization region, and we indicate thatregion by a purple line. become not causally connected by hydrodynamics. Thesecondary may accrete mass from the AGB wind with-out imposing much effect on the AGB star’s atmosphere.Such a scenario has been carefully investigated in sim-ulations that focus on the hydrodynamics around thesecondary (Huarte-Espinosa et al. 2013). Although wequalitatively know that WRLOF is different from theBHL accretion (Mohamed & Podsiadlowski 2012), thetransition from the WRLOF to the BHL accretion isstill poorly understood. One goal of this research is toillustrate such a transition with incremental change ofthe binary separation.Dust formation is the key to model the mass transferprocess in AGB binaries. When dust forms, the opac-ity of the fluid may increase by a factor of 10 , and theoutflow may become supersonic in a short time underthe radiation pressure. The location of the dust for-mation region almost determines the domain of depen-dence (circumbinary disk scenario is an exception) ofboth stars and the mass transfer mechanism. The dustformation problem has been studied from a kinetic per-spective by modeling the dust growth and sublimationfrom clusters of small particles (Gail et al. 1984; Gail& Sedlmayr 1985, 1986, 1988, 1999). The method hasbeen used by Helling & Woitke (2006) to study the dustin brown dwarfs and applied to the AGB wind problemin 1D (Jeong et al. 2003), 2D (Woitke 2006a,b), and3D (Freytag & H¨ofner 2008). More recently, H¨ofner &Freytag (2019) simulated M-type AGB stars in 3D withtheir new kinetic dust formation model (H¨ofner et al.2016). The first principle approaches are invaluable be-cause the dust formation and destruction are susceptibleto the radiation and hydrodynamical processes. How-ever, an accurate (and sensitive) dust formation modelwould intrinsically require other parts of physics such aschemistry, radiation transfer, and shock physics to be asaccurate.The chemistry plays a vital role in the formation ofdust and the dynamics of AGB winds. Boulangier et al.(2019a) modeled the AGB atmospheres and the AGBwinds in a dust-free and radiation-free environment.They emphasize the chemical non-equilibrium naturearound the AGB stars and conclude that the chemicalprocesses may regulate the thermodynamics of the AGBstars’ atmosphere, thus strongly affect the properties ofthe winds. Recently, Boulangier et al. (2019b) endeav-ored to model the onset of nucleation process from thenon-equilibrium chemical reaction. The nucleation pro-cess may explain the origin of the small particles in thekinetic model. However, the biggest obstacle of such anapproach is the lack of quantitative information aboutsome reaction rates.The radiation in the atmosphere of an AGB star maynot be in local thermal equilibrium (LTE) with thegas, and the scattering of the radiation may be non-negligible. Woitke (2006a) solved the radiative trans-fer problem in their 2D model with the Monte Carlomethod. Although Woitke (2006a) reported a negativeresult on the formation of the dusty wind, they pointedout that the composition-dependent dust formation andnon-gray radiative transfer might play a crucial role ifone wants to model the dusty wind correctly. Freytag &H¨ofner (2008) performed a 3D calculation of the AGBatmosphere with non-local radiation transfer in an ex-plicit scheme. Their results not only exhibited the dustformation in the post-shock region but also showed theAGB winds with reasonable mass-loss rates. At thesame time, radiation may affect the dynamics of thecircumbinary disks in AGB binaries but is still poorlyunderstood. A critical physics in the circumbinary disksis the evolution of the dust and the scattering of the ra-diation by the dust.The hydrodynamic problem with realistic equation ofstate (EoS) has not been studied systematically in thecontext of the WRLOF. In AGB stars’ atmospheres, theinternal energy change in phase transitions, e.g., H (cid:10) H and H (cid:10) H + + e − is significant compared to thekinetic part of thermal energy (Chen et al. 2019) andradiation energy. Freytag & H¨ofner (2008) used a Roe-type Riemann solver with a tabulated equation of stateto resolve the ionization. Phase transitions such as thedissociation of H and the ionization of H take place inthe accreting flow, too (Omukai & Nishi 1998). These phase transitions are very similar to the ones in starformation, and many similar physical phenomena mayhappen around the accreting MS stars in AGB binaries.The morphology of the outflow and the formation ofaccretion disk in evolved binary system has been studiedby many colleagues (Theuns & Jorissen 1993; Mastrode-mos & Morris 1998, 1999; de Val-Borro et al. 2009; Chenet al. 2017; de Val-Borro et al. 2017; Liu et al. 2017; Sal-adino et al. 2018, 2019; Kim et al. 2019), to name buta few. All the above mentioned works were successfulin predicting the formation of spiral structure outflow.Mastrodemos & Morris (1999) was among the first tomodel the hydrodynamics of the circumstellar materialwhich is caused by the focused wind in evolved binary.Chen et al. (2017, hereafter C17) performed a non-localradiation transfer calculation in 2D, and they observedthe formation of circumbinary disks in some AGB bina-ries. Circumbinary disk does not appear in some recentworks (de Val-Borro et al. 2017; Liu et al. 2017; Saladinoet al. 2018, 2019) where no such calculation is carriedout. In this work, we resolve the radiation transfer indifferent azimuthal angles to improve the modelling ofthe dynamics of the circumbinary disks.Numerical accuracy is of paramount importance ifsomeone tries to tackle a system that is highly non-linear. Grid-based codes in an Eulerian picture andSmoothed-Particle Hydrodynamics (SPH) (Lucy 1977;Gingold & Monaghan 1977) codes in a Lagrangian pic-ture are two popular kinds of numerical codes that aresuitable for radiation-hydrodynamic astrophysical prob-lems. With mesh-refinement techniques and high orderalgorithms, Eulerian grid-based codes are effective in re-solving shocks but are not Galilean invariant (Springel2010). The angular momentum is in general not strictlyconserved in grid-based codes. On the other hand, SPHcodes need more effort in resolving strong shocks (Taskeret al. 2008) and instabilities (Agertz et al. 2007) whilegood SPH codes can evolve angular momentum conser-vatively. Both methods have been used in the study ofthe orbital evolution in AGB binaries. Chen et al. (2018)studied the change of the binary separation by extract-ing the information of the mass transfer rate, mass-lossrate, and angular momentum loss rate in their Carte-sian grid-based radiation hydrodynamic model (C17).The error in angular momentum conservation was testedand showed that it would not affect their conclusions inChen et al. (2018). SPH method has been used to studythe evolution of the binary separation (Saladino et al.2019) and eccentricity (Saladino & Pols 2019). We findthat many previous numerical studies of AGB binarieswith SPH codes do not model the pulsating AGB winds.In our opinion, shocks are very important in AGB bi-naries. In the atmosphere of the AGB star, shocks canlevitate the atmosphere to the radius where dust couldcondensate (Lamers & Cassinelli 1999; Freytag et al.2017). Some post-shock material may fall back to theAGB star. Without shocks, the mass-loss rate of AGBwinds is usually too low. Another important shock in anAGB binary is the bow shock near the secondary. Theregion within the bow shock is a domain of dependenceof the accretion disk. Material that crosses this bowshock may be eventually accreted by the secondary.In this research, we present a 3D radiation hydrody-namic model for the AGB binary systems to illustratethe transition from the BHL accretion to the WRLOF.We carried out Four simulations on the Cedar cluster ofCompute Canada . Each simulation costs about 20 coreyears, or equivalently, about two months on 128 cores.The rest of this paper is organized as follows: we dis-cuss the adopted physics models in Section 2. Section3 outlines the governing equations of our 3D radiation-hydrodynamic binary model. Section 4 discusses thesetup of the simulation domain. Section 5 describes theboundary conditions for the AGB star and secondary.The results of the single AGB star, mass transfer ef-ficiency, outflow morphology, outflow angular distribu-tion, and orbital stability can be found in Section 6. Weconclude our work and discuss the implications of ourresults in Section 7. PHYSICS2.1.
Overall setup of the problem
We consider binary systems that consist of an AGBstar with a fixed mass of 1.02 M (cid:12) and a MS star witha fixed mass of 0.51 M (cid:12) . The separation of the binaryvaries from 5.4 AU to 6.6 AU. This range of orbital sepa-rations was found to encapsulate the transition from theWRLOF to BHL accretion, as will be discussed in Sec-tion 6.2 and 6.3. The orbit of the binary is consideredto be circular.The stellar model of the primary star is obtained withdetailed stellar evolution code and then the stellar modelis used as a boundary condition (see § MESA (Modules of Experiments in Stellar Astro-physics), release 10398 (Paxton et al. 2011, 2013, 2015,2018). We evolve a 1.5 M (cid:12) zero-age main-sequence(ZAMS) star with an initial Z = 0 .
02 and X = 0 . to 1.02 M (cid:12) . During the mass-loss, the star goes throughlengthily episodes of expansions and contractions. Thespecific model we chose for this work has an effectivetemperature T eff = 2874 K, luminosity L AGB = 4384L (cid:12) and radius of the photosphere r photo = 267 R (cid:12) .2.2. Dust formation and destruction
Dust forms at about 2-3 stellar radii of the AGB star(H¨ofner & Olofsson 2018) and may be destroyed nearthe secondary (see below). We assume that destructionand formation of dust are dependent on the gas density ρ and the gas temperature T , and dust does not exist inthe following regions:1. T > K.2. T eq > T sub , where T sub is the dust sublimationtemperature.3. ρ > − g · cm − and T > K thus destroy the dust.The physical meaning of criterion 2 is that the dustis removed if its ”radiation equilibrium” temperature ishigher than the sublimation temperature. By ”radia-tion equilibrium” temperature we mean that a sphericaldust particle is heated to the temperature that its black-body radiation can balance the radiation it absorbs if theemission opacity and the absorption opacity are equal,i.e., (cid:90) ∞ κ λ I λ πa dλ = (cid:90) ∞ κ λ B λ ( T eq )4 π a dλ, (1)where λ, κ λ , I λ , B λ are the wavelength, Planck opac-ity, incoming radiation intensity, and Planck function,respectively. a dust is the radius of the spherical dustgrain. If κ λ is a constant, Equation 1 can be simplifiedto (cid:90) ∞ I λ dλ = (cid:90) ∞ πB λ ( T eq ) dλ. (2)If only short wavelength radiation can be absorbed bythe newly formed dust, the fraction of that absorbedradiation can be expressed by α = (cid:82) λ λ πB λ ( T eff ) dλσ sb T . (3)Here σ sb is the Stefan-Boltzmann constant. For λ =0 . µ m and λ = 0 . µ m, α = 0 . composition r pyro /r photo r pyro /r photo Mg Fe SiO Fe SiO Fe SiO Fe SiO Table 1.
The second and third column show the formationradius of the test dust grain whose radius is a dust = 0 . µ m and a dust = 0 . µ m, respectively. Figure 2.
Examples of the wavelength dependence of theabsorption opacity of dust. The left y-axis is the opacityand the right y-axis is the radiation power. The black linesshow the opacity of different dust species. The blue dashedline shows the black-body spectrum for T eff = 2874 K. Grayband covers the wavelengths 0 . µ m ≤ λ ≤ . µ m andthe yellow region covers the wavelengths 2 . µ m ≤ λ ≤ . µ m. The shown opacities are found assuming the dust’snumber density distribution of dn dust ≈ a − . for dust radius0 . µ m < a dust < . µ m (Mathis et al. 1977). that could be effectively absorbed by the newly formedsmall ( a dust < . µ m) dust grains.If the material between the AGB star and the dustformation region is transparent, I λ = B λ ( T eff ). Makeuse of Equation 1 and 3 and further assume that κ λ isa constant, we can find T eq = (cid:18) αL p πσ sb r (cid:19) . , (4)where r p is the distance from the dust to the center ofthe star. A dust formation radius can be calculated if asublimation temperature T sub is given r p = √ αT r photo T . (5)When T sub = 1300 K, r dust ≈ . r photo .Indeed, κ λ for dust is a function that has a greatervalue in the short wavelength, see Figure 2 where weshow κ λ for pyroxene with different concentration ofiron. We take pyroxene as an example because it has a high sublimation temperature of about 1500 K(Kobayashi et al. 2011) and exists near M-type AGBstars (Bladh et al. 2013, 2015). We use databases ofDust Optical Properties for the optical properties ofthe specific pyroxene composition (Jaeger et al. 1994;Dorschner et al. 1995). The dust opacity is calculatedusing the Mie model (M¨atzler 2002). We use the variable κ λ and the whole spectrum to find the formation radiifor several types of dust grains, see Table 1. We can seethat the radius predicted by our simplified dust forma-tion model falls in the range of the formation radius ofpyroxene.In our simulations, r dust is found using the actual in-coming intensity. The incoming intensity may be atten-uated by the matter between the AGB star and the dustformation region, i.e., it depends on the dust-free opti-cal depth τ , which itself is the function of the directionfrom the AGB star. As a result, the location of dustformation or destruction depends on the direction fromthe AGB star and is r dust ( θ, φ ) = √ αT r photo e − τ ( θ,φ ) T , (6)where θ is the polar angle and φ is the azimuthal anglein the spherical coordinate system centered at the AGBstar’s center (see Appendix A for detail).The purpose of criterion 3 is to describe dust destruc-tion around the secondary star. As the density getshigher near the secondary, the dust and gas may comeinto LTE again and the dust sublimates.Dust formation is a complicated subject because manyspecies form at different T s, and the formation sequenceis of great importance (Jeong et al. 2003; Woitke 2006a).We anticipate that our approach is simplified. We notehowever that our goal is not to study the dust formationbut to use a relatively reasonable dust formation modelto drive a phenomenologically pulsating AGB wind.2.3. Dust and its ability to capture radiation
We integrate the blackbody spectrum for T eff = 2874K and find that within the band 0 . µ m ≤ λ ≤ . µ mit covers 75% of the total radiation power; the region2 . µ m ≤ λ ≤ . µ m covers an additional 5%. In thiswork, we focus on the absorption of small dust grains. Inthe case that the dust may grow to 1 µ m, the reddeningof the spectrum may cover the whole K band (Bladhet al. 2013, 2015). In summary, the overall momentumcontribution of the blackbody radiation beyond the Kband is relatively small, as a result, we only take 75% of the incoming radiation energy as the momentum budgetthat can be transferred to the dust; we hence choose touse L eff = f L AGB = 0 . L AGB (7)as the effective luminosity for the momentum transfer.2.4.
Dust opacity
When r dust is known, we use an opacity profile thatmimics the dust growth from the dust-free region to thedusty region κ = κ dust (cid:110) − r p − r dust − r scale r scale (cid:111) + κ mol . (8)Here κ mol = 2 . × − cm · g − is the dust-free opacityand κ dust = 6 cm · g − . r scale = 0 .
16 AU is the adoptedlength-scale which mimics the dust growth length-scale.We introduce the corresponding radius where the opac-ity of the dust-gas mixture reaches the half of its maxi-mum r half = r dust + r scale ≈ .
75 AU . (9)The choice of r scale is empirical but an order of mag-nitude estimation could be made to justify the choice.In our model, the characteristic speed of a pulsation atthe dust formation radius is v ch ≈ · s − . We canestimate the implied timescale for the dust growth as t growth ≈ r scale /v ch ≈ . . (10)The multiple 2 comes from the adopted opacity profilewhose opacity grows to half of its maximum in about r scale distance. This value is smaller than the typicaldynamical timescale and pulsation period ( ∼ Optically thin cooling
We use the following two groups of processes for thecooling process.1. The rotational and vibrational cooling of H , H O,and CO (Neufeld & Kaufman 1993).2. Cooling due to the collisional ionization, recombi-nation, collisional excitation, and Bremsstralungof H and H + (Cen 1992).We set (Neufeld & Kaufman 1993; Lacy et al. 1994) n H O /n H = 2 × − , (11) n CO /n H = 2 × − . (12) In principle, the abundance of molecules cannot beknown without carrying out chemical reaction calcula-tions. In the AGB wind, the C/O ratio is vital to thespecies of the molecules and dust (H¨ofner & Olofsson2018). The ratios we adopted in our simulations makesense only for some oxygen-rich AGB stars.We calculate n H + , n H , and n H by solving the thermalequilibrium problem, i.e., the Saha equations with thelocal ρ and T (Chen et al. 2019)H (cid:10) H + + e − , (13)H (cid:10) . (14)In summary, the cooling strength is˙Λ = ˙Λ H + ˙Λ H O + ˙Λ CO + ˙Λ H,H + . (15)Ideally, optically thin cooling should only be appliedto the optically thin region. In this research, we use ρ alone to identify whether an individual cell is opticallythick. If ρ > − g · cm − , we assume such a cell is opti-cally thick and the optically thin cooling will be turnedoff. However, with such a criterion, the vicinity of theaccretion disk may become optically thick. The cool-ing will be halted and the temperature may rise rapidlyto more than 10 K. To make the temperature profilesmoother and keep the accretion process efficient, wekeep the cooling mechanism around the secondary (in-side the accretion disk). In a real accretion disk, radia-tion (thermal and non-thermal) should be able to carryaway energy. We currently do not have such a consis-tent method to let the energy goes away, therefore, wekeep the simple cooling method. This may lead to anoverestimate of the cooling strength inside the accretiondisk. GOVERNING EQUATIONSThe governing equations are the compressible Eulerequations in a rotating frame with relevant sink andsource terms. We solve the radiation hydrodynamicequations by
AstroBEAR with static mesh refinement(Carroll-Nellenback et al. 2013). ∂ρ∂t + ∇ · ( ρ(cid:126) v ) = 0 , (16) ∂ρ(cid:126) v ∂t + ∇ · ( ρ(cid:126) v (cid:126) v ) = −∇ p + ρ ( (cid:126) a rad + (cid:126) a g + (cid:126) a i ) , (17) ∂E∂t + ∇ · [ (cid:126) v ( E + p )] = ˙Λ + ρ(cid:126) v · ( (cid:126) a rad + (cid:126) a g + (cid:126) a i ) , (18)where t, ρ, p , and E are the time, density, pressure, andtotal energy, respectively. (cid:126) v the is the velocity in therotating frame. (cid:126) a rad , (cid:126) a g , (cid:126) a i , and ˙Λ are the accelerationdue to the radiative momentum transfer, gravity of thetwo stars, and inertial force, and the energy change ratedue to the optically thin cooling. We will discuss theconsidered physics in detail in the following subsections.Also, we adopt a perfect gas EoS, p = ρk b Tµm amu , (19) (cid:15) = p ( γ − ρ , (20) E = ρ ( (cid:15) + v , (21)where k b and m amu are the Boltzmann constant andatomic mass unit. The mean atomic weight is µ =1 . γ = 5 /
3. We solvethe hydrodynamics by the Harten-Lax-van Leer-Contact(HLLC) Riemann solver with the corner transport up-wind method (Colella 1990) and piece-wise linear recon-struction.3.1.
Gravitational force and inertial force (cid:126) a g is the gravitational acceleration (cid:126) a g = − Gm p ˆ r p r − Gm s ˆ r s r , (22)where G is the gravitational constant and (cid:126) r p = (cid:126) r − (cid:126) l p , (23) (cid:126) r s = (cid:126) r − (cid:126) l s . (24)Subscript p and s stand for the primary star and the sec-ondary star, respectively. (cid:126) l p , (cid:126) l s , and (cid:126) r are the Cartesiancoordinates of the primary, the secondary, and the cell.We choose the origin to be at the center of mass of thebinary system. Self-gravity is ignored because the totalmass of the ejected material in all of our simulations isless than 5 × − M (cid:12) which is small compared to themass of the binary stars. The outflows also have somesymmetry which may cancel out most of the net effectof self-gravity. (cid:126) a i is the acceleration due to the inertial force in therotating frame. (cid:126) a i = − (cid:126) Ω b × (cid:126) v − (cid:126) Ω b × ( (cid:126) Ω b × (cid:126) r ) , (25)where (cid:126) Ω b = (0 , , ω ) is the pseudo-vector of the rotatingframe, and ω = (cid:113) G ( m p + m s ) /d (26)is the orbital frequency. d is the binary separation.Throughout our simulations, (cid:126) l p , (cid:126) l s , and (cid:126) Ω b are kept con-stant for computational simplicity. That means we donot track the orbital evolution of the binaries in our sim-ulations. We examine this assumption in Section 6.5. 3.2. Non-local momentum transfer by radiation
In this work, we only consider a single interaction be-tween the photon and dust, i.e., we adopt that the pho-tons are destroyed by the dust once they interact withthe dust. This is a reasonable approach as the dusttemperature is low, the thermal radiation from a dustgrain is less likely to be absorbed by another dust grain.On the other hand, scattering by the dust usually onlychanges the photon’s direction by a small angle (Henyey& Greenstein 1941) thus the anisotropy of the radiationfrom the AGB star may not be strong. Therefore, weonly consider the radiation in the radial direction of theAGB star. (cid:126) a rad is the acceleration due to the radiative pressurefrom the primary star. The radiation force is given by, (cid:126) a rad = κL eff e − τ ˆ r p πcr . (27)Here c is the speed of light. τ ( x, y, z ) is obtained bylinearly interpolating τ ( r, θ, φ ). κ is given by Equation8. We treat the AGB star as a point radiation source.3.3. Time-step and control of the integration
Some of the physics we used, e.g., the cooling, may op-erate on the timescale shorter than the hydro time-stepin a few cells. This makes the energy equation stiff. Wesolve the stiff ODE problem by comparing the 4th and5th order Runge-Kutta solution and adaptively evolvethe time-step (Press & Teukolsky 1992). For each inte-gration step, we only allow a maximum of 20% changein the internal energy, i.e., (cid:15) n+1 > . (cid:15) n , (28)here n is the step number of sub-cycle of the integration. THE SIMULATION DOMAINWe fix the initial mass ( m p = 1 .
02 M (cid:12) ) of theAGB star and the secondary ( m s = 0 .
51 M (cid:12) ) in allof the models. We consider four binary separations, d = 5 . , . , . , . ×
48 AU ×
24 AU in the x, y, and z directions,respectively. The base resolution of the domain is 60 × ×
30 cells and we use five levels of mesh refinement.Each level of mesh refinement doubles the resolution ineach dimension. The physical size of the finest cell is(2 . × − AU) .The actual simulation domain is a cylinder whose axiscoincides with the z-axis. The radius and the height ofthe cylinder are r domain = 23 . h domain = 24AU, respectively. This cylinder is inside the simulationbox. Radially outward supersonic flow in the lab frameis set outside the cylinder. The temperature of the su-personic flow is 1000 K and the speed is 15 km · s − .Such a supersonic flow ensures that there is no informa-tion propagating into the cylinder from its side surface.The boundary condition of the top and bottom surfacesis set diode, i.e, only outgoing flow is allowed.The refined region around the secondary is two con-axis cylinders, one with five levels of refinement has aheight of 0 . .
35 AU. The secondcylinder has four levels of refinement and has a height of0 . . .
85 AU. The sphere is aimed to resolvethe dust forming region. We resolve the equatorial zoneby using a cylinder with three levels of mesh refinement.The height and radius of the equator’s cylinder are 1 . . INNER BOUNDARY CONDITIONSThe primary star and the secondary star are two im-portant inner boundaries in our simulations. We discusstheir implementations separately.5.1.
Primary star
The primary star is the AGB star. We use the pistonmodel (Bowen 1988) to approximate the pulsating AGBstar. The boundary condition of the primary star de-termines the temperature, density, radius, velocity, andperiod of the piston. Figure 3 shows the density andtemperature profiles of the AGB star as obtained from
MESA .Similar to some other researches, we place the pistonbelow the photosphere at r piston = 247 R (cid:12) (Bowen 1988;Liljegren et al. 2016). At this radius, MESA model has T = 9484 K and ρ = 4 . × − g · cm − . The profilesand location of the piston determine the temperatureand the density of the piston, which are kept constantthroughout the simulations. The sound speed of thepiston is roughly 10 km · s − . To create a pulsating at-mosphere, we adopt a subsonic piston model. The radialvelocity of the fluid inside the AGB star is described by v piston ( r p , t ) = r p v amp r piston sin (cid:18) πtP piston (cid:19) r p ≤ r piston , (29)where v amp = 9 km · s − is the amplitude of the piston(similar amplitude has also been used in Jeong et al.(2003)) and P piston = 1 year is the adopted period of the
220 240 260 r [R ]8.48.38.28.18.07.9 l o g [ g c m ] r piston = 247 R piston = 4.77 × 10 g/cm
220 230 240 250 260 r [R ]3.53.63.73.83.94.04.1 l o g T [ K ] r piston = 247 R T piston = 9484.0 K Figure 3.
The zoomed in density and temperature profilesof the primary star produced by
MESA . The piston position isplaced at r piston = 247 R (cid:12) from the center of the AGB star. piston. The v amp we choose corresponds to 0 . Mach .We have tested smaller amplitude and find that the sub-sequent mass-loss of the AGB star is lower. In this work,we adopt the AGB star with high mass-loss rate there-fore we set the amplitude as 9 km · s − . We summarizethe parameters of the AGB star model in Table 2The AGB stars in all four binary models are not spin-ning in the lab frame. That means the AGB star iscounter-rotating in the co-rotating frame. The reasonfor setting non-spinning AGB star is twofold:1. Complete tidal spin-orbit coupling state of anAGB star with its orbital rotation is not war-ranted, and the degree of that coupling is un-known. Saladino et al. (2019) analyzed the spin-orbit coupling in AGB binaries with high mass parameter name physical quantity r photo
267 R (cid:12) T eff L AGB (cid:12) L eff (cid:12) r piston
247 R (cid:12) ρ piston . × − g · cm − T piston P piston v amp · s − Table 2.
Parameters of the AGB star model. loss rates. From Figure D.1a of their paper, wecan see that for d ∈ [5 − q = m p /m s = 2, the AGB star is spinning at80% −
90% of the orbital angular frequency. TheAGB star they analyzed has a radius of 330 R (cid:12) andis 33% larger than the boundary of the AGB star(247 R (cid:12) ) we use. If the tidal spin-orbit couplingtimescale scales with r p /d , we anticipate the spin-ning rate of our AGB star should be similar totheir d = 8 AU model which is about 60%.2. The single AGB star model we tested is a non-spinning one. We would have to test and verifymany more single star models to justify the binarysimulations with a spinning AGB star.By setting the spin of the AGB star to be zero, we areunderestimating the initial angular momentum in theoutflow. We now show that the underestimation is notbig compared to the orbital angular momentum. The z-component of the specific angular momentum in a thinshell of a sphere can be calculated by j z = 23 ωr , (30)where r piston = 247 R (cid:12) is the radius of the piston layer.The orbital specific angular momentum at d/ j p = Ω b d / . (31)In our four binary simulations, the maximum ratio of j z /j p is 0 .
27. In a not fully co-rotating case, the ratio issmaller. 5.2.
Secondary star
The secondary star is modeled as a point gravitationalsource. While it does not interact with the gas hydrody-namically, it can remove material in its vicinity that hasbecome bound to it. The mass of the removed materialis added to the point particle representing the secondary (Krumholz et al. 2004). The point particle does not havethe energy or spin angular momentum, therefore, theenergy and angular momentum in the removed materialis simply removed from the simulation. We define the’vicinity’ of the secondary to be r acc = 0 . r soft = 0 . RESULTS6.1.
Single star
We first present the results of the single AGB star withthe adopted physics and boundary conditions. Figure 4shows the time-dependent mass-loss rate and piston’sradial velocity, the wind’s radial velocity profile and es-cape velocity and the time-dependent ρ and T measuredat a distance r = 2 .
75 AU from the center of the singleAGB star. The phase difference between the mass-lossrate and the piston’s radial velocity is irrelevant becauseit simply depends on the sampling shell of the mass flux.We find that the mass-loss rate varies along with thephase of pulsation. We evaluate the mass flux at r p = 10AU. The average mass-loss rate is ˙ M single = 3 . × − M (cid:12) · year − . Due to falling back material and continuingpulsations, the gas between the dust formation regionand the piston boundary is constantly shocked, whichcan be seen in the middle panel. The wind’s velocityprofile approaches a constant value at large distances.This terminal velocity of the AGB wind is v ∞ ≈ . · s − . A good summary of the relation of the luminos-ity, effective temperature, and mass-loss rate versus thewind terminal speed of M-type AGB stars can be foundin Figure 18 in H¨ofner & Olofsson (2018). The mass-lossrate in our model is within the range of the observedmass-loss rate of 10 − − − M (cid:12) · year − . v ∞ maybea bit high. According to H¨ofner & Olofsson (2018), amore reasonable wind terminal speed should be between7 km · s − and 13 km · s − . However, this fact would notadversely affect the conclusions of this work as an AGBbinary with a slow AGB wind is more likely to form acircumbinary disk (see Section 6.2).0
15 16 17 18 19 20 t [year]2345 × [ M y e a r ] M single v piston ( r piston ) 105051015 [ k m s ] r [AU]10010 [ k m s ] v r v esc
15 16 17 18 19 20 t [year]15.515.014.514.0 l o g [ g c m ] (2.75 AU) T (2.75 AU) 2.753.003.253.503.754.00 l o g [ K ] Figure 4.
Top: the mass-loss rate of the single AGB starmodel and the radial velocity of the piston at its edge. Mid-dle: the radial velocity profile of the wind of the x-axis andthe escape speed. The AGB star is covered with the yellowband and the dust forming region is highlighted with thegray band. Bottom: the time-dependent gas density and gastemperature at r half = 2 .
75 AU on the x-axis.
Phase differences of different quantities may play im-portant roles in the mass-loss rate. Liljegren et al. (2016)studied the impact of the phase difference between thepulsation and luminosity. They find that different phaseshifts could lead to drastically different mass-loss rates.In our model, there is no phase difference in ρ and T in the piston. However, the phase between ρ (2 .
75 AU)and T (2 .
75 AU) differs by about 25% of the pulsationperiod. Because we do not consider the energy trans-fer by radiation, such a phase shift is mainly caused bycooling and shocks. According to criterion 3 of dust for-mation in Section 2.2, we can see from the third panelof Figure 4 that there are moments when the dust is de-stroyed at the radius we probe. However, we need to becareful about the temperature variation here because, in reality, H disassociation may absorb a large amount ofenergy, and the temperature may not rise to the pointthat sublimates the dust. A full consideration of the re-alistic EoS is necessary to an accurate model of the dustformation and AGB wind.6.2. Morphology of the outflow
We present the morphology of the outflows in four bi-nary simulations in Figure 5. The AGB star is in thelower half of each plot. From column 1, we can distin-guish two large-scale patterns of outflow structure: thecircumbinary disk and the spiral outflow. Specifically,circumbinary disks are found in the two simulations withthe smallest separations, while spiral outflows are foundin the two widest AGB binaries. We also find these twomorphologies in C17.At a smaller scale, we can distinguish between theappearance or not of the accretion disk around the sec-ondary. Specifically, there is no accretion disk in thesimulation with d = 6 . d = 6 . d = 5 . d = 6 . ◦ of the azimuthalangle, and its optical depth is 20. If we average its opti-cal depth over all angles, the optical depth increase willbe 1.667, which is big enough to affect the radiation-hydrodynamics in the equatorial plane.A former study of AGB binaries shows that the EoScan affect the formation of the accretion disk signifi-cantly (Theuns & Jorissen 1993). In the adiabatic case(no cooling nor heating), gas with low γ EoS can settleonto a thin disk near the accreting object, while the gaswith a high γ EoS may not be able to form an accretiondisk. In this work, as we consider cooling, the effective γ is lower than the adiabatic case in the accreting flow,1 Figure 5.
The snapshots for four binaries simulations in XY plane at Z = 0 for ρ [g · cm − ], τ , and T [K] (from the left toright). From the top to the bottom we show the simulations with d = 5 . , . , . . L potential with white or black lines and show L potential with red lines. Figure 6.
The snapshots for four binaries simulations in XZ plane at Y = 0 for ρ [g · cm − ]. The figure shows the results fromsimulations with d = 5 . , . , . . L potential with white or black lines and show L potential with red lines. a rad /a AGB where a AGB is the gravitational acceleration from theAGB star. The outflow can escape the binary system if a rad (or the luminosity of the AGB star) is large enough.If we hold the luminosity of the AGB star constant andincrease the density in the spiral outflow, the opticalthickness increases, and the value of a rad acting on theoutflow material decreases. When the optical thicknessincreases so fast that a rad near the secondary orbit be-comes too small to keep the total acceleration positive,the material in the outflow is less likely to escape fromthe binary system and a circumbinary disk forms. Thecloser the binary, the more material will be gravitation-ally focused towards the equatorial plane, and hencea circumbinary disk is more likely to form. In addi-tion, the accretion disk blocks radiation from the pri-mary star. When the material passes behind the disk’sshadow, the amount of acceleration it can receive fromradiation drops to zero, and the material is more likelyto experience a fallback towards binary. In Section 6.1,we find that the AGB wind model we adopted has a highwind speed compared to observations and some theoret-ical models. A high-speed AGB wind is less likely to becaptured or deflected by the secondary. Nevertheless,our simulations still show that circumbinary disks formin the two closest binary models. We anticipate that ifa more realistic (i.e., slower) AGB wind is being usedin our model, circumbinary disks may form in when thebinary separation wider than 5.7 AU.In this work, we do not consider the eccentricitypumping and any non-circular motion because the (un-determined) numerical viscosity in our model preventsus from modeling resonances accurately. The interac-tion between a circumbinary disk and the binary maybe one of the sources of the eccentricity (Artymowicz &Lubow 1994).6.3. Mass transfer efficiency × M y e a r M acc × M y e a r M acc × M y e a r M acc × M y e a r M acc Figure 7. ˙ M acc in the simulations with d = 5 . d = 5 . d = 6 . d = 6 . .
125 year.
We define the mass transfer efficiency as β = ˙ M acc / ˙ M p , (32)where ˙ M acc and ˙ M p are the accretion rate of the sec-ondary and the mass-loss rate of the primary. We showthe ˙ M acc during the final 40 years of simulations in allmodels in Figure 7.We find that the accretion rate varies, and the smallerthe distance between the stars, the larger and more ir-regular is the variation - ˙ M acc changes by a factor of 2for d = 5 . d = 6 . r soft = 0 . d ˙ M p ˙ M acc ˙ M flux β r half /r L1 ˙ M BHL β BHL v ∞ v orbit morpho mechanism[AU] [M (cid:12) · year − ] [M (cid:12) · year − ] [M (cid:12) · year − ] [%] [M (cid:12) · year − ] [%]5.4 3 . × − . × − . × −
31 0.921 1 . × − . × − . × − . × −
26 0.857 1 . × − . × − . × − . × − . × − . × − . × − . × − . × − Table 3.
Average mass-loss rate, accretion rate, and mass flux through the sampling shell of the four simulations. β is themass transfer efficiency. The results of the BHL accretion scenario are also calculated. v ∞ /v orbit is calculated as a reference.Morpho refers to the outflow morphology. CB and SO are short for circumbinary disk and spiral outflow which can be seen inFigure 5. Mechanism refers to the mass transfer mechanism, in our simulations, it is either the WRLOF or the BHL accretion. not smooth numerical viscosity. This numericalviscosity in multidimensional hydrodynamic simu-lations could be larger than the physical viscosity.In an accretion disk, a higher viscosity could leadto a higher accretion rate (Frank et al. 2002). Also,the small number of cells of the accretion zonearound the secondary may contribute to the ir-regularity in the accretion (Krumholz et al. 2004).2. The accretion algorithm is well tested in the BHLaccretion scenario (Krumholz et al. 2004) but isnot examined for the case of accreting from an ac-cretion disk. The direct interaction between theaccreting particle and the disk around it may beimpulsive. In a sub-Keplerian disk, a sudden re-moval of some material at the center of the diskmay lead to unbalanced forces and create artificialpressure waves (Li et al. 2003). Additionally, weuse cooling function to remove the internal energy(and the pressure) inside the disk. The cooled gasis more likely to be accreted according to the accre-tion algorithm. However, the cooling rate dependson the density and velocity and is highly nonlinear.We anticipate that sometimes the cooling makes achunk of gas to be less energetic and accreted.As we have mentioned in Section 2.5, the use of theoptically thin cooling inside the accretion disks may leadto an overestimate of the cooling rate. The impact of”over cool” depends on how much the accretion diskis pressure-supported. In a disk that is partially sup-ported by the pressure, cooling processes may decreasethe pressure gradient of the disk that could lead to acontraction of the disk and a variation in the accretionrate. The amplitude of the variation in the accretionrate depends on the strength of the cooling rate and thepressure gradient of the disk. Saladino et al. (2018);Saladino & Pols (2019) observed variations in the accre-tion rate in their SPH model with a prescribed coolingrate (their cooling rate may be much less stiff than theone used in this work). We anticipate that the accretiondisks in our simulations are partially supported by the pressure, thus, cooling contributes to the accretion ratein our simulation. An over estimation of the coolingrate may eventually lead to an over estimation of theaccretion rate and the mass transfer efficiency.The accretion rate of d = 6 . M p is the sum of the mass flux ˙ M flux thatleaves the binary system and ˙ M acc . We find average˙ M flux by summing the mass flux through a sphericalshell with r shell = 10 AU and centered at the origin.We then average ˙ M flux over 40 years, see Table 3. Atthe same time we can find what accretion rate can bepredicted by the analytic BHL accretion model˙ M BHL = G m ˙ M single v ∞ ( v + v ∞ ) / d , (33)Here v orbit = Ω b d is the relative speed between the twoorbiting objects. and v ∞ is velocity of wind at infinity.Saladino et al. (2018) multiplied the RHS of Equation 33by an efficiency parameter α BHL which may vary from0.8 to 1.5. This parameter is mainly empirical so wedo not study which value is more appropriate for whichAGB binary model. Readers can multiply β BHL in Table3 by a factor to get the corresponding range.To find what the analytic BHL accretion would pre-dict for our system, we adopt that the AGB star losesmass at the constant rate of ˙ M single and use a constantradial velocity of v ∞ = 15 . · s − which is the termi-nal velocity of the AGB wind. Equation 33 is usually theupper limit of the BHL accretion rate because its modelignores all kinds of feedback such as pressure and radi-ation of the accreting flow. The mass transfer efficiencyof the BHL accretion can be found as β BHL = ˙ M BHL ˙ M single = G m v ∞ ( v + v ∞ ) / d . (34)We can see from Table 3 that only the simulation with d = 6 . M acc and β smaller than the BHL model.5The simulation with d = 6 AU has ˙ M acc and β largerthan the BHL model. At the same time we find thatthis simulation also exhibits a spiral outflow structure.The high accretion rate suggests that WRLOF is takingplace. The spiral structure suggests that the radiationpressure is still large enough to push away the material.The two wide AGB binaries exhibit spiral outflow andhave a higher mass-loss rate as compared to the twocloser ones with circumbinary disks. It is because thecircumbinary disk confines some of the gas in the equa-torial region and the gas may fall back to the AGB star.It has been argued that v ∞ /v orbit is a key parameterof the mass transfer efficiency, and the orbital evolutionof a binary with an AGB star can be predicted by thisparameter (Saladino et al. 2018, 2019). In our four AGBbinary simulations, v ∞ /v orbit does not change much.However, the mass transfer efficiency changes drasti-cally. The change in mass transfer efficiency seems tobe closely related to the formation of the circumbinarydisk. The formation of the circumbinary disk is relatedto the increase of the optical depth in the equatorial re-gion (see Section 6.2). The formation of a circumbinarydisk in close AGB binaries questions the applicability ofthe simplified mass transfer models.6.4. Angular distribution of the outflow
The secondary deflects the outflow from the primaryto the equator, making the shape of the outflow morebipolar (see also C17). To quantitatively demonstratethat the secondary focuses the outflow onto the equato-rial plane, we split the simulation domain into six binsin the polar angle, where each bin extends for 30 ◦ . Weagain evaluate the mass flux at 10 AU from the origin.We find the time average mass of the outflow througheach of the bins. We define the fraction of mass of theoutflow in each bin as f i = ˙ M i / ˙ M p , (35)where ˙ M i is the mass of the outflow through the i th bin.Values of f i are provided in Table 4.Since the secondary is located at the equator, we dis-tribute the accreted mass of the secondary among thetwo bins in the equatorial region evenly. We find thatthe outflow of our single star is slightly not isotropic,albeit it is quite symmetric with respect to the polarangle. Its outflow is exceeding a fully isotropic case byabout 20% in the 30 − ◦ region. The Cartesian meshprobably causes anisotropy.From Table 4, one can see that the outflow is stronglyconcentrated in the equatorial region in all four binarysimulations. Let us compare the equatorial bin f forall four binary simulations. If we compare the two close binaries with the circumbinary disk, the closer one (5.4AU) has a greater f . Similarly, if we compare the twowide binaries with spiral outflow, the closer one (6.0 AU)has a greater f . If a circumbinary disk is present, somematerial in the equatorial region may fall back to theAGB star, which can be inferred in the total mass-lossrate ˙ M p in Table 3.6.5. Orbital stability
In Section 3.1, we adopted that the orbit of the bi-nary is circular and does not change. In nature, thebinary orbit would change due to the loss of the angularmomentum. Here we evaluate the potential rate of theorbital change and check posteriorly if our assumptionin Section 3.1 is justified.The total angular momentum J of the binary systemand the specific angular momentum of the binary system j are J = m p m s (cid:115) Gdm p + m s , (36) j = m p m s (cid:115) Gd ( m p + m s ) . (37)It does not include any stellar spin as in our simulationsstars were considered to not rotate. We define γ am tobe the multiple of the specific angular momentum inthe outflow in terms of the specific angular momentumof the binary system j outflow = γ am j. (38)It is a measure of the efficiency of angular momentumloss. Making use of γ am , β , and the mass ratio q = m p /m s , the rate of change of the binary separation canbe expressed by˙ dd = 2 ˙ M p m p (cid:18) − βq − (1 − β )( γ am + 12 ) q q (cid:19) . (39)We have calculated ˙ M p and β of each binary simulationin Section 6.3. The only unknown is γ am . We get γ am from the binary simulations by evaluating the average γ am through a spherical shell whose radius is 1 . d andcenter is at the center of mass of the binary. By settingthe sampling radius at 1 . d , we probably underestimatethe angular momentum loss from the binary because theescaping gas can still gain angular momentum as it goesbeyond 1 . d (Lin 1977; Saladino et al. 2018). However,the angular momentum conservation becomes worse inour code as the radius goes large (Chen et al. 2018). Weprefer not to incur too much uncertainty, so we set the6 model name f f f f f f [0 , ◦ ] [30 ◦ , ◦ ] [60 ◦ , ◦ ] [90 ◦ , ◦ ] [120 ◦ , ◦ ] [150 ◦ , ◦ ]isotropic 6 . × − . × − . × − . × − . × − . × − single 5 . × − . × − . × − . × − . × − . × − . . × − . × − . × − . × − . × − . × − . . × − . × − . × − . × − . × − . × − . . × − . × − . × − . × − . × − . × − . . × − . × − . × − . × − . × − . × − Table 4.
The fraction of mass of outflow through binned polar angular range. The second row list the polar angle range of eachbin in degree. 0 ◦ corresponds to the positive z direction. The third row shows the angular distribution if the flux from a sphereis isotropic. The fourth row shows the result from the single star from Section 6.1. The results of the four binary simulationsare listed from row 4 to row 8. sampling shell small. The specific angular momentumof each star are, j p , s = m , p (cid:115) Gd ( m p + m s ) (40)Therefore, by taking Equation 37 into consideration, wecan find two limiting cases for γ am . When q = 2 and γ am = 0 .
5, the gas emitted from the AGB star does notgain any angular momentum, when γ am ≈
2, the escap-ing gas has a specific angular momentum that similar tothe secondary. In Table 5, we list values of γ am for eachof our simulations. The obtained values are in a betweenof the two extreme cases. As the distance increases, γ am decreases towards the specific angular momentumof AGB star.At the risk of introducing too many variables (readerswill appreciate this variable later), we convert the γ am derived from our simulations to a new quantity definedby η = qγ am (1 + q ) = Jm p m s ( m p + m s ) . (41)The Equation 7 of Saladino et al. (2019) is an empiricalformula that uses q and v ∞ /v orbit to predict the value of η em . Table 5 lists γ am , η , and η em of each of our binarysimulation and the potential rate of change of the orbit.It turns out that η derived from each of our simulationis indeed close to η em . The match suggests that eventhough v ∞ /v orbit may not be effective in characterizingthe mass transfer efficiency when a circumbinary diskis about to form (see Section 6.3), v ∞ /v orbit could beuseful in predicting the angular momentum loss fromthe AGB binary.From Table 5, we can infer that the binary separationin our four binary simulations should not change morethan 5 × − AU. Therefore, our assumption in Sec-tion 3.1 is reasonable. On the other hand, if the binaryseparation change rate is kept constant for 10 years,which corresponds to the average lifetime of an AGBstar, the decrease in binary separation of the two close d γ am δγ am η η em ˙ d/d ∆ d [AU] [%] [ × − / yr] [AU]5 . − . − . . − . − . . . . . .
816 0 . Table 5.
The measured γ am of the four binary models andrelative rate of the change of the binary separation. δγ am =(( γ am − . / . × years if the binaryseparation change rate is kept constant. binary is non-negligible. Considering that we have ig-nored the spin-orbit coupling of the binary and we set asmall sampling shell to get γ am , we are underestimatingthe orbital angular momentum loss. The actual orbitalshrinkage should be higher. As the binary separationdecreases, the WRLOF may become more similar to theRLOF. CONCLUSIONS AND DISCUSSIONSIn this paper, we carry out 3D radiation-hydrodynamic simulations of AGB binaries. We useray-tracing radiation transfer to resolve the opticaldepth in the azimuthal and polar directions. We takeinto account that dust can be destroyed by shocks, ina high temperature close to an LTE environment, andif radiation from the AGB star can sublimate the smalldust grains (see Section 3.2 for detail). We calculate thecooling strength on H , H, and H + ; their number densi-ties are found by solving the Saha equations (Chen et al.2019). We model the AGB star as an inner boundarycondition that has sinusoidally varying radial velocity(Bowen 1988). We use MESA to obtain the temperatureand density that determine the boundary conditions ofour piston model.7The single star model presented in Section 6.1 has anaverage mass-loss rate of 3 . × − M (cid:12) · year − anda terminal wind speed of 15 . · s − . We find thatthe mass-loss rate is reasonable for an M-type AGB starwith an effective temperature of 2874 K and a lumi-nosity of 4384 L (cid:12) (H¨ofner & Olofsson 2018). Observa-tions show that the terminal wind speed should be 7-13km · s − (H¨ofner & Olofsson 2018). We argued that thisdifference would not adversely affect our conclusion. Onthe contrary, lower terminal wind speed may be cap-tured or deflected by the secondary more easily. In thecase of low wind speed, a circumbinary disk may formwhen the binary separation is greater than 6.0 AU (seeSection 6.2 for detail).We find that the accretion rate increases when the bi-nary separation decreases (see Table 3). The mass trans-fer efficiency may increase dramatically to 31% when acircumbinary disk forms. Such a mass transfer efficiencycan be up to 8 times of what the canonical BHL ac-cretion rate predicts. The presence of a circumbinarydisk imposes severe challenges to the simplified masstransfer mechanisms. It enlarges the domain of depen-dence of the secondary to the whole equatorial region.The circumbinary disk may also exert a torque on thebinary and change the eccentricity of the binary sys-tem (Artymowicz & Lubow 1994; Dermine et al. 2013;Moody et al. 2019). The eccentricity problem has beenheavily studied by the proto-planetary disk community,and similar conditions may apply here (Ragusa et al.2018). The formation of a circumbinary disk is a re-sult of the focused fluids in the equatorial plane by thesecondary (Section 6.4). The focused fluids increase theoptical depth in the equatorial plane; thus, the radiationpressure on the fluids decreases. The balance betweenthe radiation pressure and the gravity will be tilted to-ward the gravity when the optical depth becomes largerthan a critical value, and some of the gas may fall back.We point out that the non-local momentum transfer byradiation is crucial in the formation of the circumbinarydisks. It would be difficult to observe the formation ofthe circumbinary disks without the non-local momen-tum transfer even when WRLOF takes place (Mohamed& Podsiadlowski 2012).The accretion rates in the two closest AGB binariesare irregular. We outlined two possible reasons for theirregularity in Section 6.3. It may be worthwhile to in-crease the resolution and allow a subsonic atmosphere(or envelope) to build up around the secondary in futureresearches. The atmosphere can provide a subsonic re-gion that may smooth out the supersonic flow that fallsonto the secondary. The particle can accrete mass fromthe atmosphere when the atmosphere becomes Jeans un- stable (Jeans 1902; Federrath et al. 2010). It is alsonatural to model proper feedback, i.e., radiative cool-ing, from the secondary when an atmosphere presents.In our current model, we have not taken radiative feed-back from the accretor and the radiation transfer in theaccretion disk into consideration. In a realistic circum-stance, the accretion disk will become optically thickwhen the density of the disk becomes high enough. Thetransition from the optically thin to the optically thickstate may incur different accretion modes. A better un-derstanding of the accretion rate may be achieved bymodeling the transition correctly.The two closest binaries have dusty circumbinarydisks, which may resemble the situation AR Puppis (Er-tel et al. 2019). The two widest binaries have spiralstructure outflows. Such structure resembles many post-AGB binaries with wide binary separations or planetarynebulae (Edgar et al. 2008; Maercker et al. 2012; Ram-stedt et al. 2017; Kim et al. 2019).In Section 6.5, we discussed the orbital dynamics ofthe AGB binaries in our simulations. We confirmed thatit is reasonable to ignore the orbital evolution in oursimulations. We also find that the two closest AGB bi-naries may experience an orbital shrinkage, where thebinary with the initial separation of 5 . . γ am , we are underestimating the orbitalangular momentum loss in close binaries. The decreaseof the binary separation in the two close binaries shouldbe faster than what we estimated. The binary with 5 . − − − g · cm − in the middle plane near the bi-nary orbit) in the circumbinary disks can also foster dustgrowth. Large dust grains may process photons with awavelength longer than the K band. If multiple scatter-ing by the dust becomes essential, a more sophisticatednon-gray radiative transfer model is in need. In oursimulations, we find that the optical depth as calculatedalong the line of sight from the AGB star is typically lessthan 10 (except for the accretion disk). Monte Carlo ra-8diative transfer would be an efficient way to model theradiation-hydrodynamics in the circumbinary disk.In summary, the mass transfer efficiency β and angularmomentum loss efficiency γ am are two important quan-tities but with large uncertainty in AGB binaries. Theyhugely affect the binary evolution. In this work, we de-rived β and γ am from 3D radiation-hydrodynamic simu-lations. We discover a huge discrepancy (up to 8 times)in β by comparing our simulations with the canonicalBHL accretion mechanism. The discrepancy is closelyrelated to the presence of the circumbinary disks. It isessential to carry out the non-local radiation transfercalculation to model the formation and dynamics of acircumbinary disk. Software:
AstroBEAR (Carroll-Nellenback et al.2013), matplotlib (Hunter 2007),
MESA v10398 (Pax-ton et al. 2011, 2013, 2015, 2018),
Visit 3.0.1 (Childset al. 2012). ACKNOWLEDGMENTSWe thank the anonymous referee whose suggestionshave greatly improved the quality of this work. Wethank Craig Heinke, Rodrigo Fernandez, Xuening Bai,Falk Herwig, Pavel Denisenkov, Nami Mowlavi, Ue-LiPen, and Dong Lai for inspirational discussions. ZCis grateful to the CITA National Postdoctoral Fellow-ship and the Tsung-dao Lee visiting scholarship. N.I.acknowledges support from CRC program and fund-ing from NSERC Discovery. This research was enabledby the use of computing resources provided by Com-pute/Calcul Canada, and was supported in part by theNational Science Foundation under Grant No. NSFPHY-1748958.REFERENCES
Abate, C., Pols, O. R., Izzard, R. G., & Karakas, A. I.2015, A&A, 581, A22, doi: 10.1051/0004-6361/201525876Abate, C., Pols, O. R., Izzard, R. G., Mohamed, S. S., & deMink, S. E. 2013, A&A, 552, A26,doi: 10.1051/0004-6361/201220007Agertz, O., Moore, B., Stadel, J., et al. 2007, MNRAS, 380,963, doi: 10.1111/j.1365-2966.2007.12183.xArtymowicz, P., & Lubow, S. H. 1994, ApJ, 421, 651,doi: 10.1086/173679Beers, T. C., & Christlieb, N. 2005, ARA&A, 43, 531,doi: 10.1146/annurev.astro.42.053102.134057Bidelman, W. P., & Keenan, P. C. 1951, ApJ, 114, 473,doi: 10.1086/145488Bladh, S., H¨ofner, S., Aringer, B., & Eriksson, K. 2015,A&A, 575, A105, doi: 10.1051/0004-6361/201424917Bladh, S., H¨ofner, S., Nowotny, W., Aringer, B., &Eriksson, K. 2013, A&A, 553, A20,doi: 10.1051/0004-6361/201220590Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273,doi: 10.1093/mnras/104.5.273Boulangier, J., Clementel, N., van Marle, A. J., Decin, L.,& de Koter, A. 2019a, MNRAS, 482, 5052,doi: 10.1093/mnras/sty2560Boulangier, J., Gobrecht, D., Decin, L., de Koter, A., &Yates, J. 2019b, arXiv e-prints, arXiv:1908.09633.https://arxiv.org/abs/1908.09633Bowen, G. H. 1988, ApJ, 329, 299, doi: 10.1086/166378Carroll-Nellenback, J. J., Shroyer, B., Frank, A., & Ding,C. 2013, Journal of Computational Physics, 236, 461,doi: 10.1016/j.jcp.2012.10.004 Cen, R. 1992, ApJS, 78, 341, doi: 10.1086/191630Chen, Z., Blackman, E. G., Nordhaus, J., Frank, A., &Carroll-Nellenback, J. 2018, MNRAS, 473, 747,doi: 10.1093/mnras/stx2335Chen, Z., Coleman, M. S. B., Blackman, E. G., & Frank, A.2019, Journal of Computational Physics, 388, 490,doi: 10.1016/j.jcp.2019.03.016Chen, Z., Frank, A., Blackman, E. G., Nordhaus, J., &Carroll-Nellenback, J. 2017, MNRAS, 468, 4465,doi: 10.1093/mnras/stx680Childs, H., Brugger, E., Whitlock, B., et al. 2012, in HighPerformance Visualization–Enabling Extreme-ScaleScientific Insight, 357–372Colella, P. 1990, Journal of Computational Physics, 87, 171,doi: 10.1016/0021-9991(90)90233-QDahn, C. C., Liebert, J., Kron, R. G., Spinrad, H., &Hintzen, P. M. 1977, ApJ, 216, 757, doi: 10.1086/155518de Val-Borro, M., Karovska, M., & Sasselov, D. 2009, ApJ,700, 1148, doi: 10.1088/0004-637X/700/2/1148de Val-Borro, M., Karovska, M., Sasselov, D. D., & Stone,J. M. 2017, MNRAS, 468, 3408,doi: 10.1093/mnras/stx684Dermine, T., Izzard, R. G., Jorissen, A., & Van Winckel, H.2013, A&A, 551, A50, doi: 10.1051/0004-6361/201219430Dorschner, J., Begemann, B., Henning, T., Jaeger, C., &Mutschke, H. 1995, A&A, 300, 503Edgar, R. 2004, NewAR, 48, 843,doi: 10.1016/j.newar.2004.06.001Edgar, R. G., Nordhaus, J., Blackman, E. G., & Frank, A.2008, ApJL, 675, L101 Ertel, S., Kamath, D., Hillen, M., et al. 2019, AJ, 157, 110,doi: 10.3847/1538-3881/aafe04Escorza, A., Karinkuzhi, D., Jorissen, A., et al. 2019, A&A,626, A128, doi: 10.1051/0004-6361/201935390Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S.2010, ApJ, 713, 269, doi: 10.1088/0004-637X/713/1/269Frank, J., King, A., & Raine, D. J. 2002, Accretion Powerin Astrophysics: Third EditionFreytag, B., & H¨ofner, S. 2008, A&A, 483, 571,doi: 10.1051/0004-6361:20078096Freytag, B., Liljegren, S., & H¨ofner, S. 2017, A&A, 600,A137, doi: 10.1051/0004-6361/201629594Gail, H. P., Keller, R., & Sedlmayr, E. 1984, A&A, 133, 320Gail, H. P., & Sedlmayr, E. 1985, A&A, 148, 183—. 1986, A&A, 166, 225—. 1988, A&A, 206, 153—. 1999, A&A, 347, 594Gezer, I., Van Winckel, H., Manick, R., & Kamath, D.2019, MNRAS, 488, 4033, doi: 10.1093/mnras/stz1967Gingold, R. A., & Monaghan, J. J. 1977, MNRAS, 181,375, doi: 10.1093/mnras/181.3.375Giridhar, S., Rao, N. K., & Lambert, D. L. 1994, ApJ, 437,476, doi: 10.1086/175011Hansen, T. T., Andersen, J., Nordstr¨om, B., et al. 2016,A&A, 588, A3, doi: 10.1051/0004-6361/201527409Helling, C., & Woitke, P. 2006, A&A, 455, 325,doi: 10.1051/0004-6361:20054598Henyey, L. G., & Greenstein, J. L. 1941, ApJ, 93, 70,doi: 10.1086/144246Hillen, M., Kluska, J., Le Bouquin, J. B., et al. 2016, A&A,588, L1, doi: 10.1051/0004-6361/201628125H¨ofner, S., Bladh, S., Aringer, B., & Ahuja, R. 2016, A&A,594, A108, doi: 10.1051/0004-6361/201628424H¨ofner, S., & Freytag, B. 2019, A&A, 623, A158,doi: 10.1051/0004-6361/201834799H¨ofner, S., & Olofsson, H. 2018, A&A Rv, 26, 1,doi: 10.1007/s00159-017-0106-5Homan, W., Richards, A., Decin, L., et al. 2017, A&A, 601,A5, doi: 10.1051/0004-6361/201630340Hoyle, F., & Lyttleton, R. A. 1939, Proceedings of theCambridge Philosophical Society, 35, 405,doi: 10.1017/S0305004100021150Huarte-Espinosa, M., Carroll-Nellenback, J., Nordhaus, J.,Frank, A., & Blackman, E. G. 2013, MNRAS, 433, 295,doi: 10.1093/mnras/stt725Hunter, J. D. 2007, Computing in Science & Engineering, 9,90, doi: 10.1109/MCSE.2007.55Jaeger, C., Mutschke, H., Begemann, B., Dorschner, J., &Henning, T. 1994, A&A, 292, 641 Jeans, J. H. 1902, Philosophical Transactions of the RoyalSociety of London Series A, 199, 1,doi: 10.1098/rsta.1902.0012Jeong, K. S., Winters, J. M., Le Bertre, T., & Sedlmayr, E.2003, A&A, 407, 191, doi: 10.1051/0004-6361:20030693Jorissen, A., Boffin, H. M. J., Karinkuzhi, D., et al. 2019,A&A, 626, A127, doi: 10.1051/0004-6361/201834630Jorissen, A., Van Eck, S., Van Winckel, H., et al. 2016,A&A, 586, A158, doi: 10.1051/0004-6361/201526992Karambelkar, V. R., Adams, S. M., Whitelock, P. A., et al.2019, ApJ, 877, 110, doi: 10.3847/1538-4357/ab1a41Keenan, P. C. 1942, ApJ, 96, 101, doi: 10.1086/144435Kervella, P., Montarg`es, M., Lagadec, E., et al. 2015, A&A,578, A77, doi: 10.1051/0004-6361/201526194Khouri, T., Velilla-Prieto, L., De Beck, E., et al. 2019,A&A, 623, L1, doi: 10.1051/0004-6361/201935049Kim, H., Liu, S.-Y., & Taam, R. E. 2019, ApJS, 243, 35,doi: 10.3847/1538-4365/ab297eKobayashi, H., Kimura, H., Watanabe, S. i., Yamamoto, T.,& M¨uller, S. 2011, Earth, Planets, and Space, 63, 1067,doi: 10.5047/eps.2011.03.012Krumholz, M. R., McKee, C. F., & Klein, R. I. 2004, ApJ,611, 399, doi: 10.1086/421935Lacy, J. H., Knacke, R., Geballe, T. R., & Tokunaga, A. T.1994, ApJL, 428, L69, doi: 10.1086/187395Lamers, H. J. G. L. M., & Cassinelli, J. P. 1999,Introduction to Stellar WindsLi, L.-X., Goodman, J., & Narayan, R. 2003, ApJ, 593, 980,doi: 10.1086/376695Liljegren, S., H¨ofner, S., Nowotny, W., & Eriksson, K. 2016,A&A, 589, A130, doi: 10.1051/0004-6361/201527885Lin, D. N. C. 1977, MNRAS, 179, 265,doi: 10.1093/mnras/179.2.265Liu, Z.-W., Stancliffe, R. J., Abate, C., & Matrozis, E.2017, ApJ, 846, 117, doi: 10.3847/1538-4357/aa8622Lucy, L. B. 1977, AJ, 82, 1013, doi: 10.1086/112164Maercker, M., Mohamed, S., Vlemmings, W. H. T., et al.2012, Nature, 490, 232, doi: 10.1038/nature11511Manick, R., Van Winckel, H., Kamath, D., Hillen, M., &Escorza, A. 2017, A&A, 597, A129,doi: 10.1051/0004-6361/201629125Mastrodemos, N., & Morris, M. 1998, ApJ, 497, 303,doi: 10.1086/305465—. 1999, ApJ, 523, 357, doi: 10.1086/307717Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ,217, 425, doi: 10.1086/155591M¨atzler, C. 2002, IAP Res. Rep, 8, 9McClure, R. D., & Woodsworth, A. W. 1990, ApJ, 352,709, doi: 10.1086/168573 Mohamed, S., & Podsiadlowski, P. 2012, Baltic Astronomy,21, 88, doi: 10.1515/astro-2017-0362Moody, M. S. L., Shi, J.-M., & Stone, J. M. 2019, ApJ,875, 66, doi: 10.3847/1538-4357/ab09eeMowlavi, N., Lecoeur-Ta¨ıbi, I., Lebzelter, T., et al. 2018,A&A, 618, A58, doi: 10.1051/0004-6361/201833366Neufeld, D. A., & Kaufman, M. J. 1993, ApJ, 418, 263,doi: 10.1086/173388Omukai, K., & Nishi, R. 1998, ApJ, 508, 141,doi: 10.1086/306395Oomen, G.-M., Van Winckel, H., Pols, O., & Nelemans, G.2019, A&A, 629, A49, doi: 10.1051/0004-6361/201935853Oomen, G.-M., Van Winckel, H., Pols, O., et al. 2018,A&A, 620, A85, doi: 10.1051/0004-6361/201833816Ortiz, R., & Guerrero, M. A. 2016, MNRAS, 461, 3036,doi: 10.1093/mnras/stw1547Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192,3, doi: 10.1088/0067-0049/192/1/3Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208,4, doi: 10.1088/0067-0049/208/1/4Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS,220, 15, doi: 10.1088/0067-0049/220/1/15Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS,234, 34, doi: 10.3847/1538-4365/aaa5a8Podsiadlowski, P., & Mohamed, S. 2007, Baltic Astronomy,16, 26Press, W. H., & Teukolsky, S. A. 1992, Computers inPhysics, 6, 188Ragusa, E., Rosotti, G., Teyssandier, J., et al. 2018,MNRAS, 474, 4460, doi: 10.1093/mnras/stx3094 Ramstedt, S., Mohamed, S., Vlemmings, W. H. T., et al.2017, A&A, 605, A126,doi: 10.1051/0004-6361/201730934Roulston, B. R., Green, P. J., Ruan, J. J., et al. 2019, ApJ,877, 44, doi: 10.3847/1538-4357/ab1a3eSahai, R., Findeisen, K., Gil de Paz, A., & S´anchezContreras, C. 2008, ApJ, 689, 1274, doi: 10.1086/592559Saladino, M. I., & Pols, O. R. 2019, A&A, 629, A103,doi: 10.1051/0004-6361/201935625Saladino, M. I., Pols, O. R., & Abate, C. 2019, A&A, 626,A68, doi: 10.1051/0004-6361/201834598Saladino, M. I., Pols, O. R., van der Helm, E., Pelupessy,I., & Portegies Zwart, S. 2018, A&A, 618, A50,doi: 10.1051/0004-6361/201832967Springel, V. 2010, MNRAS, 401, 791,doi: 10.1111/j.1365-2966.2009.15715.xStarkenburg, E., Shetrone, M. D., McConnachie, A. W., &Venn, K. A. 2014, MNRAS, 441, 1217,doi: 10.1093/mnras/stu623Tasker, E. J., Brunino, R., Mitchell, N. L., et al. 2008,MNRAS, 390, 1267,doi: 10.1111/j.1365-2966.2008.13836.xTheuns, T., & Jorissen, A. 1993, MNRAS, 265, 946,doi: 10.1093/mnras/265.4.946Van der Swaelmen, M., Boffin, H. M. J., Jorissen, A., &Van Eck, S. 2017, A&A, 597, A68,doi: 10.1051/0004-6361/201628867Van Winckel, H., Waelkens, C., Waters, L. B. F. M., et al.1998, A&A, 336, L17Vlemmings, W., Khouri, T., O’Gorman, E., et al. 2017,Nature Astronomy, 1, 848,doi: 10.1038/s41550-017-0288-9Woitke, P. 2006a, A&A, 460, L9,doi: 10.1051/0004-6361:20066322—. 2006b, A&A, 452, 537, doi: 10.1051/0004-6361:20054202
Figure 8.
An illustrative picture of the spherical mesh of theequatorial region. Angular resolution doubles in both polarand azimuthal coordinates in the outer region. The figure isonly for illustration, the resolution of our simulation is shownin Table 6. A.
OPTICAL DEPTH τ ( x, y, z ) is calculated by interpolating τ ( r, θ, φ ) lin-early after translating ( x, y, z ) to the frame whose originis the center of the primary star. τ ( r, θ, φ ) is the opti-cal depth trace back to an angular dependent surface r surf ( θ j , φ k ) that encloses the photosphere and probablythe chromosphere of the primary star. For simplicity, r surf is defined as the smallest contour of ρ = 10 − g · cm − that encloses the AGB star. τ is defined on adiscrete spherical coordinate τ ( r i , θ j , φ k ) = (cid:40) τ inside r i < r surf , (cid:82) r i r surf κ ( r, θ j , φ k ) ρ ( r, θ j , φ k ) dr r i ≥ r surf , (A1)where θ j ∈ [0 , π ] is the discretized polar angle and φ k ∈ [0 , π ] is the discretized azimuthal angle. τ inside is a large number but is irrelevant to the actual cal-culation because a rad is set to 0 within r surf . To dis-tinguish, we will call a finite volume in Cartesian coor-dinate a ’cell’ and a finite volume in spherical coordi-nate a ’bin’ as shown in Figure 8. The center of thisspherical coordinate is at the center of the primary startherefore the coordinate is translated. κ ( r i , θ j , φ k ) and ρ ( r i , θ j , φ k ) are mapping and regriding of the data inthe translated Cartesian coordinate. To distinguish thetranslated Cartesian coordinate from the untranslatedcoordinate, we use ( x (cid:48) , y (cid:48) , z (cid:48) ) to denote the translatedCartesian coordinate. The density of a bin in the spher-ical coordinate is defined as ρ ( r i , θ j , φ k ) = (cid:80) ρ ( x (cid:48) , y (cid:48) , z (cid:48) ) dv ( x (cid:48) , y (cid:48) , z (cid:48) ) δ ( x (cid:48) , y (cid:48) , z (cid:48) ) (cid:80) dv ( x (cid:48) , y (cid:48) , z (cid:48) ) δ ( x (cid:48) , y (cid:48) , z (cid:48) ) , (A2)where dv is the finite volume in the translated coor-dinate. To achieve a high precision of ρ ( r i , θ j , φ k ), we apply 4 levels of mesh refinement to all the quantitiesthat are defined on the cells, that means each dv has aphysical size of (0 . . δ ( x (cid:48) , y (cid:48) , z (cid:48) ) can be thoughtof as a kernel function, i.e., δ ( x (cid:48) , y (cid:48) , z (cid:48) ) (cid:40) x (cid:48) , y (cid:48) , z (cid:48) ) ∈ dv ( r i , θ j , φ k ) , x (cid:48) , y (cid:48) , z (cid:48) ) (cid:54)∈ dv ( r i , θ j , φ k ) . (A3)The membership relation ∈ can be defined explicitly as, r i-1/2 ≤ r < r i+1/2 , (A4) θ j-1/2 ≤ θ < θ j+1/2 , (A5) φ k-1/2 ≤ φ < φ k+1/2 , (A6)where ( r, θ, φ ) is the spherical coordinate of ( x (cid:48) , y (cid:48) , z (cid:48) ).Similarly, we define the opacity of a bin as massweighted average κ ( r i , θ j , φ k ) = (cid:80) ρ ( x (cid:48) , y (cid:48) , z (cid:48) ) κ ( x (cid:48) , y (cid:48) , z (cid:48) ) dv ( x (cid:48) , y (cid:48) , z (cid:48) ) δ ( x (cid:48) , y (cid:48) , z (cid:48) ) (cid:80) ρ ( x (cid:48) , y (cid:48) , z (cid:48) ) dv ( x (cid:48) , y (cid:48) , z (cid:48) ) δ ( x (cid:48) , y (cid:48) , z (cid:48) ) . (A7)In the mapping and regriding process, the computa-tional domain is divided into two regions, one is thepolar region and the other is the equatorial region. Thepolar region does not resolve the azimuthal angle, there-fore, in the polar region, Equation A3 becomes δ ( x (cid:48) , y (cid:48) , z (cid:48) ) (cid:40) x (cid:48) , y (cid:48) , z (cid:48) ) ∈ dv ( r i , θ j ) , x (cid:48) , y (cid:48) , z (cid:48) ) (cid:54)∈ dv ( r i , θ j ) . (A8)The equatorial region resolves the azimuthal angle. An-gular resolution is increased for both two regions as r gets larger. We illustrate the spherical grid that is usedfor the radiative transfer in Figure 8. The resolutionstructure that we use is listed in Table 6. Inner region Refined region r ∈ [1 . , .
5] AU r ∈ [2 . ,
36] AUdr d θ d φ dr d θ d φ . π/ π/
50 0 . π/ π/ Table 6.
Angular resolution of the equatorial region in ra-diative transfer. The polar region has the same angular res-olution except that azimuthal angle is not resolved.
In Section 2.2, we calculate the dust-free optical depthby setting κ = κ mol = 2 . × − g · cm − everywhere.The τ ( r, θ, φr, θ, φ