Formation of the Asymmetric Accretion Disk from Stellar Wind Accretion in an S-type Symbiotic Star

In order to simulate the wind accretion processes in S-type symbiotic stars, we adopt the adaptive mesh refinement code FLASH (Fryxell et al. 2000) in version 4.5, with modifications similar to those in the works carried out by Kim & Taam (2012b) and Kim et al. (2013, 2017, 2019). In this work, the fluid is assumed to be inviscid and non-self-gravitating. The Eulerian hydrodynamic equations are used in the center-of-mass frame, with the Cartesian coordinates.

Based on the continuity equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

the initial density distribution is set by

$$\begin{eqnarray}&&{\rho }_{0}({\boldsymbol{r}})=\displaystyle \frac{{\dot{M}}_{{\rm{G}}}}{4\pi | {\boldsymbol{r}}-{{\boldsymbol{r}}}_{{\rm{G}}}{| }^{2}\,{V}_{\mathrm{ej}}},\end{eqnarray} \tag{ 2 }$$

assuming a steady-state stellar wind originating from the intrinsically spherical symmetric mass loss of a nonpulsating and nonrotating giant star situated at r = r _G. The initial density distribution of the background is relaxed within an orbit of the binary in the simulation. Here, the parameter V_ej indicates the stellar wind velocity at ∣ r − r _G∣ = r_*, which we define as the surface of the giant star. The radius r_*, at which the wind density and velocity conditions are reset for every simulation step, is set to 0.25 au, similar to the photospheric radius of an evolved giant star of an S-type symbiotic system (Skopal 2005; see also Dumm & Schild 1998). Because the radius of the giant star is much less than the distance to the inner Lagrangian point (∼1 au) in the binary system employed in our simulations, mass transfer through the Roche lobe overflow does not occur. The adopted mass-loss rate of the giant star is ${\dot{M}}_{{\rm{G}}}={10}^{-7}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}$. Because the mass loss of the giant star during the total simulation time is negligible, in terms of the orbital evolution of the binary system, a simplification is made, by fixing the stellar masses and their orbits throughout the simulation.

The masses of the two stars are fixed to M_G = 1.5 M_⊙ and M_WD = 0.6 M_⊙. The two stars orbit the center of mass of the binary system, following the assumed circular trajectories with a fixed separation of 2 au (=∣ r _G − r _WD∣). The corresponding orbital period is ∼2 yr. The orbital velocity of the WD is V_WD = 21.8 km s⁻¹, with respect to the center of mass of the binary, and the relative velocity between the stars is V_rel = 30.5 km s⁻¹. In Table 1, we summarize the fixed parameters adopted for the models exhibited in this paper.

Table 1. Fixed Model Parameters Adopted in the Hydrodynamic Simulations in This Work

Parameter	Value	Description
MG	1.5 M⊙	Mass of the giant star
${\dot{M}}_{{\rm{G}}}$	10−7 M⊙ yr−1	Mass-loss rate of the giant star
MWD	0.6 M⊙	Mass of the WD
a	2 au	Binary separation
r* & rsoft,G	0.25 au	Photospheric and gravitational softening radii of the giant star
Teff	4,000 K	Effective temperature of the giant star
f	1	Acceleration factor

Download table as:  ASCIITypeset image

The acceleration force on the wind material is determined by the equation of motion:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t}+({\boldsymbol{v}}\cdot {\rm{\nabla }}){\boldsymbol{v}}=-\displaystyle \frac{1}{\rho }{\rm{\nabla }}P-{\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{G}}}^{\mathrm{eff}}-{\rm{\nabla }}{{\rm{\Phi }}}_{\mathrm{WD}}.\end{eqnarray} \tag{ 3 }$$

An adiabatic equation of state with a specific heat ratio, γ = 5/3, for a monatomic ideal gas is used to update the pressure P. We also set the effective temperature of the giant star, T_eff = 4000 K, which is proper for the giant star of S-type symbiotic stars, as presented by Skopal (2005).

The specific gravitational forces $-{\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{G}}}^{\mathrm{eff}}$ and − ∇Φ_WD, due to the presence of the giant star and the WD, are defined by the Plummer models with the softening radii r_soft,G and r_soft,WD, respectively (Binney & Tremaine 2008):

$$\begin{eqnarray*}&&\begin{array}{l}-{\rm{\nabla }}{{\rm{\Phi }}}_{{\rm{G}}}^{\mathrm{eff}}=-\displaystyle \frac{{{GM}}_{{\rm{G}}}(1-f)}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{{\rm{G}}}{| }^{2}+{r}_{\mathrm{soft},{\rm{G}}}^{2}},\\ -{\rm{\nabla }}{{\rm{\Phi }}}_{\mathrm{WD}}=-\displaystyle \frac{{{GM}}_{\mathrm{WD}}}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}_{\mathrm{WD}}{| }^{2}+{r}_{\mathrm{soft},\mathrm{WD}}^{2}}.\end{array}\end{eqnarray*}$$

Here, G, M_G,WD, and r _G,WD are the gravitational constant, the individual masses of the giant star and the WD, and the position vectors of the stars at the center-of-mass coordinates, respectively. The softening radius r_soft,G for the giant star is set to 0.25 au, being comparable to the expected photospheric radius of the giant star (see above). The size of the WD (on the order of 10⁻⁵ au) cannot be resolved in our simulations, due to the limitations on computation time; therefore, we set the softening radius r_soft,WD as the smallest length scale at a given resolution of the simulation model, within which we require more than a minimum of four grid cells (see Table 2). The acceleration factor f of unity is assumed in all of our simulations (see Section 2.1 of Kim & Taam 2012a for details of the wind model).

The orbital velocity of the WD relative to the wind velocity in situ is an important factor regulating the detailed properties of the accretion disk. The simulations are carried out by varying the intrinsic wind velocity, V_w, defined as the velocity at the distance to the WD (i.e., 2 au) in the corresponding single-star simulation, in which the WD is intentionally absent. Alternatively speaking, V_w is determined merely by the balance between the gas pressure and the effective gravitational potential of the giant star with an inclusion of the acceleration factor, mimicking the radiation pressure onto the hypothetical dust particles. The velocity V_w for each model is listed in Table 2.

In the energy equation, we include the radiative cooling as two separate terms, following Saladino et al. (2018):

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot [(E+P){\boldsymbol{v}}]=-\rho \left[\displaystyle \frac{3k}{2\mu {m}_{{\rm{p}}}}\displaystyle \frac{(T-{T}_{\mathrm{eq}})\rho }{C}+{\dot{Q}}_{\mathrm{sh}}\right].\end{eqnarray} \tag{ 4 }$$

Here, E, T, P, v , and ρ are the basic quantities that describe the fluid, representing the total energy, temperature, pressure, velocity vector, and mass density, respectively. The constants k and m_p are the Boltzmann constant and proton mass, respectively, and μ = 1.29 is the mean molecular weight appropriate to the solar metallicity (e.g., Anders & Grevesse 1989).

The first cooling term on the right-hand side of Equation (4) regulates the gas temperature, being balanced with the equilibrium temperature (Chandrasekhar 1934),

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}={T}_{\mathrm{eff}}{\left(W+3/4{\tau }^{{\prime} }\right)}^{1/4},\end{eqnarray} \tag{ 5 }$$

within the equilibrium timescale, C/ρ, where the constant C = 10⁻⁵ g s cm⁻³ is derived from an approximate density-dependent expression of the radiative cooling rate in the atmosphere of giant stars (Bowen 1988). The effective temperature T_eff represents the temperature at the photosphere of the giant star. The parameter $W={(1/2)[1-(1-{r}_{* }^{2}/{s}^{2})}^{1/2}]$ indicates the geometrical dilution factor at a distance s = ∣ r − r _G∣ from the giant star, while the optical depth ${\tau }^{{\prime} }={\int }_{s}^{\infty }({\kappa }_{{\rm{g}}}+{\kappa }_{{\rm{d}}})\rho \,{r}_{* }^{2}/{s}^{2}\ {ds}$ and the factor 3/4 in Equation (5) are determined by assuming a plane-parallel atmosphere. The opacities for gas and dust, κ_g and κ_d, are 2 × 10⁻⁴ and 5 cm² g⁻¹, respectively. Notice that when the gas temperature drops below T_eq, this first cooling term behaves like a heating source with regard to the radiation of the giant star, preventing unrealistic reductions in the gas temperature.

The second term on the right-hand side of Equation (4) describes the radiative cooling, taking into account neutral hydrogen and heavier chemical elements in a gaseous medium with solar abundance. The cooling rate is defined as

$$\begin{eqnarray}&&{\dot{Q}}_{\mathrm{sh}}=\displaystyle \frac{{{\rm{\Lambda }}}_{\mathrm{hd}}{n}_{{\rm{H}}}^{2}}{\rho },\end{eqnarray} \tag{ 6 }$$

with the cooling function Λ_hd provided by Schure et al. (2009) in their Table 2. The hydrogen number density n_H = X ρ/m_p is calculated with a fixed mass fraction of hydrogen in the gas (X = 0.707). When the gas temperature rapidly increases at the shock fronts that are established at the spiral's arms or in the vicinity of the WD, under the strong gravitational potential well, the second cooling term plays an important role in efficiently reducing the internal energy of the gas. It therefore provides a favorable condition for disk formation surrounding the WD. It should be noted that the low-temperature regime (${\mathrm{log}}_{10}T\lt 3.8$) is not considered in this radiative cooling term, and thus the gas temperature in such a regime stays near the equilibrium temperature, as determined by the first cooling term in Equation (4).

The WD is treated as an accretion sink, with the sink radius r_sink, in order to prevent the unrealistic penetration of flows through the WD and to obtain legitimate values for the accretion rate. In our simulation, the radius r_sink, within which the accretion sink occurs, is set to be the same as the gravitational softening radius r_soft,WD. When the flow approaches the WD within the radius r_sink and its energy condition satisfies the gravitationally bound state, the flow is regarded as having entered the sink, and no further trace is made. Instead, the density within the sink radius is replaced by an arbitrary small value, the sink density ρ_sink, which is predefined to be much smaller than the background density, practically as the in situ density in the absence of the WD. The velocity within the sink radius is replaced by the orbital velocity of the WD. The accretion rate, ${\dot{M}}_{\mathrm{acc}}={\sum }_{i}{\rho }_{i}\,{{dU}}_{i}/{dt}$, is recorded at the end of each simulation time step, of the time duration dt, where the density and volume elements of the ith grid cell within the accretion sink, satisfying the sink condition, are denoted by ρ_i and dU_i . Several tests have been performed to verify that our choice for ρ_sink does not affect the overall morphology of the accretion flows and their physical properties. We have also performed tests to confirm the sufficiency of the numerical resolution in assigning at least four cells within the accretion sink radius (see also Truelove et al. 1997). In tests for the sink radius reduced by half (and accompanied by proper grid resolution), negligible changes are observed in the disk radius and height measurements, and in the density and temperature ranges within the accretion disk.

The computational domain is set to be a rectangular cube, with the dimensions 12.8 × 12.8 × 6.4 au³, and with the center of mass of the binary stars (x = y = z = 0) at the center of the first and second coordinates, but on the boundary of the third coordinate, assuming a symmetry about the equatorial plane (z = 0). The boundary condition for the plane coinciding with the equatorial plane reflects all quantities in the z > 0 space to the z < 0 space through the mirror at the z = 0 plane. The boundary conditions for the other five boundary planes are set to diode, defined in the FLASH code, which allows the outflow but prohibits the inflow of fluid. At each simulation time step, the meshes for the high-density regions are refined through the adaptive mesh refinement technique, with a predefined maximum refinement level of 6 (or a level of 7 in one model, for the resolution test) and 64 × 64 × 32 cells being included in each level of refinement. To interpolate the physical quantities of each grid cell, the Piecewise Parabolic Method is adopted in this work.

In the analyses for the accretion disk of the WD, we introduce the coordinates (x', y', z'), defined in the rest frame of the WD. We also note that the notations r and R indicate the radii in spherical and cylindrical coordinates, respectively. The angle Θ is defined about the WD in the counterclockwise direction, starting from the opposite side to the location of the giant star.

We have carried out eight simulations, Sim1–Sim8, with different wind conditions, and their input parameters and the resulting properties of the accretion disks are tabulated in Table 2. Our experiments are designed to investigate the physical properties of the outflows and the accretion disks surrounding the WD, according to the increase in V_w from Sim1 to Sim6. The intrinsic wind velocity of Sim1, measured from an unperturbed single-star simulation, starts from V_ej =6 km s⁻¹ at the surface of the giant star (i.e., at radius r_*) and reaches V_w = 13.7 km s⁻¹ at the distance corresponding to the binary separation. This wind model mimics a transonic wind, as the sound speed at the surface of the giant star is ∼7.4 km s⁻¹. The winds of the giant star in the other seven models are intrinsically supersonic. A slow wind with V_w smaller than the orbital velocity V_WD = 21.8 km s⁻¹ (Sim1–Sim5) facilitates the capture of stellar wind material, resulting in favorable conditions for the formation of an accretion disk surrounding the WD.

A fast wind with the velocity V_w exceeding the orbital velocity V_WD tends to pass over the WD, instead of effectively building a disk around it. Under such unfavorable conditions for disk formation, the Sim6 model failed to create a disk. With reduced sink radii, the Sim7 and Sim8 models were computed in order to assess the effects of the sink radius and the grid resolution on the formation of the accretion disk. We assured our minimum grid criterion, four cells within radius r_sink, for these models by increasing the maximum refinement level.

In the comparison of all the simulated models, the Sim5 model exhibits a transitional morphology (see Figure 1). The bow shock surrounding the accretion disk is well developed in the forward direction of the orbital motion of the WD in the slower Sim1–Sim4 models, while the material above the disk radius and beneath the bow shock is stripped out in the faster Sim5–Sim8 models. As a result, the accretion disks in the Sim5–Sim8 models are exposed to winds, and the accretion tails (or accretion columns) are revealed in the manner that is demonstrated in the BHL theory (Bondi & Hoyle 1944). Notice that morphological transition occurs with the in situ intrinsic wind velocity being similar to the orbital velocity of the WD (V_w ∼ V_WD).

Zoom In Zoom Out Reset image size

Figure 2 presents, as an example, the detailed structures of the Sim2 model at large and small scales at the time when the stars have completed their four orbits. The first and second columns of the figure exhibit the density and temperature distributions, respectively, viewed face-on (upper panels) and edge-on (lower panels). The third and fourth columns show zoomed-in images of the same distributions near the WD. The binary stars are rotating in the counterclockwise direction around their center of mass, at (x, y, z) = (0, 0, 0). The giant star and the WD are located at (x_G, y_G) = (−0.57 au, −0.01 au) and (x_WD, y_WD) = (1.42 au, 0.02 au), respectively, at the current epoch.

Zoom In Zoom Out Reset image size

It is clearly identified that a spiral structure attached to the WD coils around the central stars (see also Theuns & Jorissen 1993; Kim & Taam 2012b; Huarte-Espinosa et al. 2013; Saladino et al. 2018, 2019). The inner and outer edges of the spiral structure are well defined by the narrow regions showing high temperatures (T > 8000 K; see Figure 2(b)). The outer edge is extremely turbulent, and it creates a thick wall filled with high-temperature substructures. The inner edge is relatively smoothly distributed, in both density and temperature. This inner edge is tightly wound around the giant star, terminating at the meeting point with the outer edge of the spiral, (x, y) ∼ (1.8 au, 0.2 au), as seen in Figures 2(a) and (c). Toward this position, the temperature of the inner edge is no longer high (see Figure 2(b)). The region around the giant star confined by the inner spiral edge is relatively low in density.

Most of the wind material entering into the spiral structure eventually escapes from the binary potential, while some fraction of the material remains gravitationally captured by the WD, forming an accretion disk. The blue lines in the third and fourth columns of Figure 2 demonstrate the 3D morphology of the accretion disk (as defined in Section 3.2). Interestingly, the accretion disk is vertically not razor-thin, but has a flare shape (Figures 2(g)–(h)).

The formation criterion for an accretion disk adopted in this work is the presence of circularized accretion flows around the WD that last longer than an orbital period. If formed, the shape and mass of the accretion disk are computed as follows.

• 1.  
  Disk height H_disk. We first draw 10 density profiles along vertical lines passing the points of intersection with the equatorial plane at an arbitrary radius R from the WD, starting from r_sink and increasing with intervals of ΔR. The gas density decreases, following an exponential function up to a certain height, and beyond this height the density drops abruptly. We define the scale height H(R) at radius R as the average of the 10 values for these heights that characterize the exponential density declines. The disk height H(R_disk), or simply H_disk, is the scale height at the radius of the disk, R = R_disk, defined as below.
• 2.  
  Disk mass M_disk. We estimate the mass M(R) within the volume of the flare-shaped disk, bounded by a cylinder with an arbitrary radius R and scale height H(R), with its inner boundary at the sphere of the radius r_sink. We repeat this mass estimation by gradually raising the radius of the outer boundary by ΔR, until the mass increment is negligible (less than 1%), i.e., M(R + ΔR) − M(R) <0.01M(R + ΔR), and the final value is defined as the disk mass M_disk, or M (R_disk).
• 3.  
  Disk radius R_disk. This is correspondingly defined as the radius at which the above mass condition is satisfied.

We estimate these three quantities at each simulation time step. The values show some fluctuations in the early orbits, and become stabilized in the later orbits. In Figure 3, the R_disk, H_disk, and M_disk of, for example, the Sim2 model are displayed as a function of time. The disk radius R_disk increases slightly beyond the noise level of variation, until it eventually reaches a near-constant value of ∼0.15 au at t > 4 orbit. The disk height H_disk shows a ∼40% increase in the early evolution (t < 4 orbit) and attains a near-constant value of ∼0.03 au. The disk mass reaches up to M_disk ∼ 4.6 × 10⁻⁸ M_⊙ at t ∼ 5.4 orbit, and stops increasing afterward. The final disk mass corresponds to ∼3% of the mass expelled from the giant star over 14 yr (seven orbits). We note that the growth of the disk mass within the nearly stationary volume of the disk at t = 4–5.4 orbit enhances the density of the disk (see the change in the minima of the density profiles with time in Section 4.1). The mean values of these quantities over the last orbit of the simulation are listed in Table 2.

Zoom In Zoom Out Reset image size

Disks are always formed, regardless of the initial parameter sets, if the simulation mesh resolution and the accretion sink radius are properly assigned. Under the same spatial conditions as the Sim1–Sim5 models, the Sim6 model fails to create an accretion disk (see Figure 1). It turns out to be a numerical artifact, as the Sim7 model with a sink radius r_sink (=r_soft,WD) reduced by half successfully forms an accretion disk with the radius of 0.1 au. From these experiments, we note that a sink radius greater than ∼25% of the (potential) disk radius causes the artificial destruction of the disk. Here, we did not change the resolution of the simulation grids, as they satisfied the minimum requirements for resolving r_sink (i.e., four cells within r_sink). The Sim8 model is for the disk formation in an environment where the adjacent wind velocity is faster than the orbital velocity of the WD. In this model, by reducing r_sink and r_soft,WD to 0.0125 au, we successfully produced the accretion disk with its radius of 0.058 au.

As noted in Table 2, the representative disk radius (i.e., the average disk radius over the temporal oscillation during the last simulation orbit) does not change with wind velocity in the Sim1–Sim4 models, but starts to decrease with faster winds in the Sim5–Sim8 models. The disk radius <R_disk>, tabulated in Table 2 and indicated by the filled circles in Figure 4, is smaller than the accretion-line impact parameter b (the black solid line in Figure 4) derived by Huarte-Espinosa et al. (2013) for the BHL flow toward the retarded position of the WD:

$$\begin{eqnarray}&&b=2\,{q}^{2}{\left(1+q\right)}^{3}{\left(\displaystyle \frac{1}{a}\right)}^{5/2}{\left(\displaystyle \frac{1}{a}+\displaystyle \frac{1}{{a}_{w}}\right)}^{-7/2},\end{eqnarray} \tag{ 7 }$$

where q = M_WD/M_G is the stellar-mass ratio, a is the binary separation, and ${a}_{w}={{GM}}_{{\rm{G}}}/((1+q){V}_{{\rm{w}}}^{2})$ is a specific length scale indicating the binary separation when the wind velocity V_w equals the orbital velocity of the WD about the center of mass of binary stars, i.e., V_w = V_WD(a = a_w ). The disk radii of the Sim1–Sim4 models are almost the same, and significantly smaller than the b parameter, while the disk radii of the Sim5–Sim8 models decrease, following the decreasing trend of the b line with wind velocity. On the other hand, the bars above the filled circles in Figure 4 present the distances to the outermost rotating component of the material around the WD, measured along the line toward the giant star; in most cases, it represents the stand-off distance of the bow shock. Figure 4 indicates that the b parameter actually limits this distance, which, by definition, is larger than the disk radius. For reference, the Hill radius and the sink radius of each model are also drawn in Figure 4, in order to show the definitely forbidden area for the disk radius.

Zoom In Zoom Out Reset image size

The mass-accretion rate ${\dot{M}}_{\mathrm{acc}}$ is measured by integrating the mass entering into the sink sphere with its energy condition that forbids its escape from the sink. For comparison, we also measure the mass transfer rate ${\dot{M}}_{\mathrm{inflow}}$, by integrating the mass fluxes inflowing through the surface of the flared disk (defined in Section 3.2). The disk mass M_disk is independently measured in Section 3.2, by integrating the density distribution within the flared disk, excluding the sink sphere. In our fiducial Sim2 model, the averaged rates over the entire simulation time are ${\dot{M}}_{\mathrm{inflow}}\sim 4.8\times {10}^{-8}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}$ and ${\dot{M}}_{\mathrm{acc}}\sim 2.3\times {10}^{-8}\ {M}_{\odot }\,{\mathrm{yr}}^{-1}$, implying that the inflowing matter to the disk is twice the matter swallowed by the WD. Because ${\dot{M}}_{{\rm{i}}{\rm{n}}{\rm{f}}{\rm{l}}{\rm{o}}{\rm{w}}}$ exceeds ${\dot{M}}_{\mathrm{acc}}$, the inflowing matter is gradually accumulated within the disk, and therefore the disk mass increases with the evolution, which is consistent with Figure 3(c).

Figure 5 presents the time-averaged value of the mass-accretion efficiency β, which is defined as the mass-accretion rate into the sink of the WD divided by the mass-loss rate from the giant star. To avoid spikes in the measurements of the accretion rates, the averages and errors are estimated in the median base. The accretion efficiency is higher than 10% in the Sim1–Sim3 models, and it rapidly decreases with the wind velocity, up to ∼18 km s⁻¹. The accretion efficiency remains at the level of a few percent in the Sim4–Sim8 models. These measurements in our simulation models are compared with the accretion efficiency under the BHL assumption (Boffin & Jorissen 1988):

$$\begin{eqnarray}&&{\beta }_{\mathrm{BHL}}=\displaystyle \frac{{\alpha }_{\mathrm{BHL}}}{2}{\left(\displaystyle \frac{q}{1+q}\right)}^{2}{\left(\displaystyle \frac{{V}_{\mathrm{rel}}}{{V}_{{\rm{w}}}}\right)}^{4}{\left[1+{\left(\displaystyle \frac{{V}_{\mathrm{rel}}}{{V}_{{\rm{w}}}}\right)}^{2}\right]}^{-3/2},\end{eqnarray} \tag{ 8 }$$

where V_rel = (1 + q)V_WD is the orbital velocity of the WD in the rest frame of the giant star, and the efficiency coefficient α_BHL = 1 is adopted. It is found that the accretion efficiencies in the slow-wind models (Sim1–Sim3) are considerably larger than β_BHL, while those in the faster-wind models (Sim4–Sim8) are similar to the BHL prediction.

Zoom In Zoom Out Reset image size

Espey & Crowley (2008) pointed out the poor understanding of the mass-loss processes in cool giants, quoting the wind terminal velocities <100 km s⁻¹ (see also, e.g., Dupree & Reimers 1987). In the case of S-type symbiotic stars with wind velocity exceeding 20 km s⁻¹, our simulations imply that the accretion efficiency will be a few percent, comparable with that of the BHL accretion. However, a higher accretion efficiency of ∼10% or more is expected for S-type symbiotic systems with a lower wind velocity than the orbital velocity of the WD. The orbital parameters for most D-type symbiotic stars are only poorly known, so that a reliable estimate of the accretion efficiency based on the current work is not obvious (e.g., Schmid & Schild 2002; Matthews & Karovska 2006; Hinkle et al. 2013).

Of particular interest is the presence of two inflowing spirals within the accretion disk (see, e.g., Figure 2(c)). In cataclysmic variables, the disk spirals are excited by tidal interaction, as has been confirmed by many observations through indirect imaging techniques, such as Doppler tomography, and by inviscid hydrodynamic calculations at several nonadiabatic equations of state (Matsuda et al. 2000, and references therein). In symbiotic stars, however, an observational confirmation of the disk spirals has not been made, and a theoretical approach to the disk spirals has only been performed by Bisikalo et al. (1997). In this work, we stress that disk spirals can be formed even in symbiotic stars through the 3D wind accretion calculation for the first time.

We present, for example, maps from the Sim2 model in Figure 6, with panel (a) showing the density in linear color scale, panel (b) showing the radial velocity V_R , and panel (c) showing the rotational velocity V_ϕ of the matter. The magnitudes of the directional velocities in the maps in panels (b) and (c), and the fluid velocity vectors denoted by the gray arrows, are all recalculated in the rest frame of the WD at the current time t = 4 orbit. One spiral entering into the disk radius at Θ ∼ 30°, i.e., from the forward direction of the orbital motion of the WD (hereafter, the front spiral), is winding the WD for more than one lap, before reaching the sink radius (see the black solid curve in Figures 6(a)–(c)). The other spiral presents a shorter trajectory, starting from R = R_disk at Θ ∼ −40° (hereafter, the rear spiral; see the solid gray curve).

Zoom In Zoom Out Reset image size

In order to measure the physical quantities along the spirals, we first define the shapes of the disk spirals in a functional form of the logarithmic spiral: $R/{R}_{\mathrm{disk}}=\exp [-({\rm{\Theta }}-{{\rm{\Theta }}}_{a})/{{\rm{\Theta }}}_{b}]$, where the Θ_a parameter satisfies Θ_a = Θ(R = R_disk) and the polar slope of the spiral is a constant $-{{\rm{\Theta }}}_{b}^{-1}$. After unrolling the density map of the accretion disk in the polar coordinates, we find the positions of the local density maxima along the radius axis, within a 5 pixel margin from the test function of the logarithmic spiral. The χ² goodness-of-fit test provides the coefficients (Θ_a , Θ_b ) = (33° ± 4°, 500° ± 10°) for the front spiral and (Θ_a , Θ_b ) = (−41° ± 3°, 332° ± 8°) for the rear spiral, where the radial distance R from the WD is defined from the disk radius R_disk = 0.15 au to the sink radius r_sink = 0.05 au.

Figure 6(d) shows that the densities along these spirals are ∼5 × 10⁻¹² g cm⁻³ at R ∼ R_disk (marked by the open circles), which are similar to each other. The density along the rear spiral increases up to ∼2 × 10⁻¹¹ g cm⁻³ at Θ ∼ 180°, and then steeply decreases until the spiral reaches R = r_sink. On the other hand, the density along the front spiral increases up to ∼2.5 × 10⁻¹¹ g cm⁻³ at Θ ∼ 220°–250°, and then decreases until the spiral reaches the sink radius.

It is noticeable that the density along the front spiral ceases the decreasing trend at Θ ∼ 300°–360° (see Figure 6(d)). Interestingly, in a similar range of Θ, the radial velocity of the matter along the front spiral becomes positive at Θ ∼ 220°–370° (see Figure 6(e)), and the rotational velocity ceases its increasing trend by staying at V_ϕ ∼ 60 km s⁻¹ at Θ = 200°–360° (see Figure 6(f)). Our best interpretation for these abnormal behaviors of the front spiral is that the flows entering the disk from the front side, with respect to the orbital motion of the WD, are not well circularized, but are slightly overshot toward the rear side (Θ ∼ 270°), perhaps due to the complex process during the convergence of the streamlines (see Section 3.5). The overshot matter raises the neighboring density beyond its ordinary value, and then the matter is diverted onto the outer disk region, before falling back onto its ordinary path beyond Θ ∼ 360°. In contrast, the rear spiral does not show such a drifting event, and its clear inflowing feature is implied by the gray-colored profiles in Figures 6(d)–(f).

A close inspection of the streamlines reveals that there are two major streams acting as the origin of the disk material. For a streamline analysis, we trace the streamlines, starting from 400 random points far from the disk. When the flows arrive at the disk radius with negative (i.e., incoming) radial velocities, we record those entering spots at the disk radius R_disk in Figure 7(a). The arrows in this figure show the magnitudes and directions of the fluid velocities in the rest frame of the WD.

Zoom In Zoom Out Reset image size

We find that the flows do not enter the disk through the entire outer rim of the disk, but that the entering spots are instead limited to the two arc segments marked by the blue and red colors in Figure 7(a). The angular extent of the blue arc segment ranges from 36° to 216°, spanning a half circle, whereas the red arc segment has a small angular extent of ΔΘ ∼ 3°, limited to the range between 336° and 339°.

Figure 7(b) shows six representative accretion streamlines that arrive at the outer rim of the disk. It is apparent that the three red streamlines converge at the red arc segment shown in Figure 7(a). The three blue streamlines represent wind components originating directly from the giant star, which subsequently get deflected as they pass through the bow shock. The entering spots of this blue group of flows are spread over ΔΘ ∼ 180°, as marked by the blue arc segment in Figure 7(a). The inflows from the rear side amount to ∼50% of the total incoming mass, which means that a similar amount of mass enters the disk from the front side through the bow shock.

The two groups of inflows appear to be assimilated into the two spiral features after entering the disk. The rear spiral is likely composed of the inflowing matter entering through the small arc segment, marked by the red spots in Figure 7(a), whose angular position is similar to Θ_a of the rear spiral. The inflows from the front side through the large arc segment (the blue spots) would compose the front spiral, and the incessant injection of matter throughout a large fraction of the spiral would be responsible for the abnormal behavior of the front spiral described in the previous subsection.

One of the notable features in the disk shape is the asymmetric density distribution. Figure 8 presents the time evolutions of the density maps of the Sim2 disk model in the equatorial plane. We find that the overall density distribution in the disk is conspicuously asymmetric and that the highest density is always achieved near the direction toward the giant star. In addition, the entering spots of the accretion streamlines do not show significant changes in the azimuthal angular coordinate, Θ, measured with respect to the direction toward the giant star being Θ = 180°. The two spirals repeat the acts of merging (at ∼4.65 orbit and ∼5.40 orbit) and splitting (at ∼4.00 orbit and ∼5.05 orbit) along their evolution, which shows a tendency for the highest degree of asymmetry when they merge.

Zoom In Zoom Out Reset image size

Figure 9 shows the density profiles averaged over the radius in the range r_sink < R < R_disk for the Sim2 model. Overall, the density profiles are single-peaked at Θ ∼ 218°, 175°, 229°, and 189° at t = 4, 5, 6, and 7 orbits, respectively. As we analyzed in Section 3.4, such a density enhancement at the direction toward the giant star, or on the slightly rare side with respect to the orbital motion, probably originates from the overshooting of the accretion flows composing the front spiral. The maximum densities are ∼50% higher than the minima on the opposite sides of the disk. This implies that, if the emissivity depends only on the squares of density, the emissivity ratio between the maximum and the minimum in a disk exceeds ∼2, which is consistent with the best-fitting result for the accretion disk of AG Dra hypothesized by Lee et al. (2019).

Zoom In Zoom Out Reset image size

We investigate the model dependence of the degree of asymmetry, defined as the ratio between the maximum and the minimum of the radially averaged density profile along the azimuthal angle of the disk, ρ_Max/ρ_Min. In the upper panel of Figure 10, the filled circle presents the mean value of the time variation of this quantity during the last orbit of the simulation, and the error bar shows its standard deviation. The resulting density ratio is nearly independent of the wind models, being ∼1.5 on average.

Zoom In Zoom Out Reset image size

The large temporal variations in the degree of density asymmetry, as indicated by the error bars, have an implication for an observation at a marginal spatial resolution: the time variation of the emissivity from the dominant part of the disk may be a natural consequence and contribute to the variability, in addition to the pulsation of the giant star, the oscillation during the binary orbital period, and the intrinsic variation in the stellar wind (e.g., Gális et al. 2015; Richie et al. 2020).

The bottom panel of Figure 10 shows that the angular direction of ρ_Max tends to be toward the giant star (Θ ∼ 180°) in the relatively slow-wind models (Sim1–Sim4), in which relatively large disks are formed. The deviation from Θ ∼ 180° with time is, however, quite large, probably due to the turbulent motions of the inflowing matter, as expected from the fluctuations in the outer outflowing spiral enveloping the two stars. For the faster-wind models (Sim5–Sim8), a more refined numerical approach, substantiated with increased computing power, should be adopted in order to clarify the physical origin of the large uncertainties.

In addition to the matter distribution around the WD component, the kinematic property of the accretion flow is essential to properly model the spectroscopic observations. Lee et al. (2019) investigated the line profile of Raman-scattered O vi features at 6825 and 7082 Å, based on the assumption that the O vi emission region of AG Dra is characterized by Keplerian motion around the WD component (e.g., Heo et al. 2016, 2021). However, one would naturally expect the real accretion flow to be much more complicated than a simple Keplerian motion, due to the presence of the giant component and the contributions from the thermal and turbulent motions.

An analysis of Sim2 reveals that the azimuthal velocity component V_ϕ of the accretion stream is about 90% of the Keplerian velocity V_K. More specifically, at R = 0.135 au, the azimuthal speed of the accretion flow V_ϕ is measured to be 53 km s⁻¹, where the expected Keplerian speed is V_K =61 km s⁻¹ (see Figure 6(f)). The sub-Keplerian rotation is probably attained as a result of the combined effect of the turbulent component and the pressure gradients associated with the adiabatic equation of state, which act against the gravity of the WD. The sub-Keplerian accretion flow highlights an important caveat in interpreting the spectroscopic observations of symbiotic stars, including AG Dra. That is, one may overestimate the size of the Raman O vi emission region when it is deduced under the simple assumption that the flow is Keplerian.

For example, in the spectrum of AG Dra investigated by Lee et al. (2019), the separation of the peaks in the double-peak profile of Raman-scattered O vi at 6825 Å was observed to be ≃8.2 Å, which was translated to a Keplerian velocity V_K ≃ 27 km s⁻¹ of the main O vi emission region around the WD. With the adopted value of M_WD = 0.6 M_⊙, the distance to the O vi emission region from the WD is 0.7 au, assuming that the flow is Keplerian. If the accretion flow is sub-Keplerian, as the result of the Sim2 model indicates, then the O vi emission region would be located at a distance of r_OVI = 0.5 au from the WD, which is smaller by a factor of ∼70% than the value deduced under the Keplerian flow assumption.

However, the disk radius for the Sim2 model is <R_disk> = 0.16 au, which is only one-third of the above estimate of r_OVI. This substantial discrepancy may indicate other possibilities for radiative transfer modeling, with O vi emission regions and neutral regions more sophisticated than those of Lee et al. (2019). For example, based on the spectra obtained with the Far Ultraviolet Spectroscopic Explorer, Young et al. (2005) proposed that the O vi emission-line region is located in the illuminated part of the giant atmosphere. It has also been proposed that the O vi emission region may be found in a slowly expanding region from the hot WD (Schmid et al. 1999). It should be noted that the radiative transfer of Raman-scattered O vi by Lee et al. (2019) was based on an additional assumption, that all the neutral matter was present near the giant star. Furthermore, Schmid (1996) proposed that the kinematics of the neutral wind from the giant component in symbiotic stars is very important in the line formation of Raman-scattered O vi features, pointing out that the receding part of the neutral wind contributes to the red-enhanced line profile exhibited in Raman O vi features.