Stellar Wind Accretion and Raman Scattered O VI Features in the Symbiotic Star AG Draconis

TL;DR

This study uses high-resolution spectroscopy and Monte Carlo simulations to analyze Raman scattered O VI features in the symbiotic star AG Draconis, revealing asymmetric accretion flows and bipolar outflows related to the white dwarf's accretion process.

Contribution

It provides the first detailed modeling of Raman O VI features in AG Draconis, linking line profiles to accretion flow asymmetry and bipolar outflows, and estimates the mass loss rate of the giant component.

Findings

Raman features exhibit double profiles with enhanced red parts.

Accretion flow is azimuthally asymmetric with more matter on the entering side.

Mass loss rate estimated at approximately 2 x 10^{-7} solar masses per year.

Abstract

We present high resolution spectroscopy of the yellow symbiotic star AG Draconis with ESPaDOnS at the {\it Canada-France-Hawaii Telescope}. Our analysis is focused on the profiles of Raman scattered \ion{O}{VI} features centered at 6825 \AA\ and 7082 \AA, which are formed through Raman scattering of \ion{O}{VI}$\lambda\lambda$1032 and 1038 with atomic hydrogen. These features are found to exhibit double component profiles with conspicuously enhanced red parts. Assuming that the \ion{O}{vi} emission region constitutes a part of the accretion flow around the white dwarf, Monte Carlo simulations for \ion{O}{VI} line radiative transfer are performed to find that the overall profiles are well fit with the accretion flow azimuthally asymmetric with more matter on the entering side than on the opposite side. As the mass loss rate of the giant component is increased, we find that the flux ratio $F(6825)/F(7082)$ of Raman 6825 and 7082 features decreases and that our observational data are consistent with a mass loss rate $\dot M\sim 2 \times 10^{-7} {\rm\ M_{\odot}\ yr^{-1}}$. We also find that additional bipolar components moving away with a speed $\sim 70{\rm\ km\ s^{-1}}$ provide considerably improved fit to the red wing parts of Raman features. The possibility that the two Raman profiles differ is briefly discussed in relation to the local variation of the \ion{O}{VI} doublet flux ratio.

Tables

Table 1. Table 1: Model parameters of AG Dra adopted in this work.

Parameter	Value	Description	Reference
$M_G$	1.5 M⊙	Mass of giant	Skopal (2005)
$R_G$	35 R⊙	Radius of giant	Hric et al. (2014)
$M_{\mathrm{WD}}$	0.5 M⊙	Mass of white dwarf	Mikołajewska et al. (1995)
a	400 R⊙	Binary separation	Derived from Fekel et al. (2000)
$L_1$	140 R⊙	Inner Lagrange point
v∞	10, 15, 20 $\mathrm{km}{s}^{-1}$	Wind terminal velocity
$\dot{M}$	10-5 - ${10}^{-7}$ $M_{\odot }{\mathrm{yr}}^{-1}$	Mass loss rate

Table 2. Table 2: O  vi local emissivity ϵ ​ ( 𝐫 ) italic-ϵ 𝐫 \epsilon({\bf r}) as a function of the azimuthal angle in the accretion disk.

${\varphi }_j$	$0$	$\pi /4$	$\pi /2$	$3\pi /4$	$\pi$	$5\pi /4$	$3\pi /2$	$7\pi /4$
$A_j$	1.69	2.19	2.5	2.19	1.69	1.13	1	1.13

Table 3. Table 3: Flux ratio F ​ R ​ ( ϕ ) = F ​ ( 1032 ) / F ​ ( 1038 ) 𝐹 𝑅 italic-ϕ 𝐹 1032 𝐹 1038 FR(\phi)=F(1032)/F(1038) as a piecewise step function of azimuthal angle ϕ italic-ϕ \phi in the O  vi annular region.

${\varphi }_j$	$0$	$\pi /4$	$\pi /2$	$3\pi /4$	$\pi$	$5\pi /4$	$3\pi /2$	$7\pi /4$
$B_j$	1.3	1.14	1	1.14	1.3	1.6	2	1.6

Equations

$${\lambda }_i^{-1}={\lambda }_o^{-1}+{\lambda }_{\alpha }^{-1},$$

$$\frac{\Delta {\lambda }_o}{{\lambda }_o}=\left(\frac{{\lambda }_o}{{\lambda }_i}\right)\frac{\Delta {\lambda }_i}{{\lambda }_i},$$

$$r_i=7\times 10^{11}\left(\frac{M_{\mathrm{WD}}}{0.5{\mathrm{M}}_{\odot }}\right)\mathrm{cm},$$

$$X=\frac{4\pi {\mu }^2{m}_{\mathrm{p}}^2}{\alpha }a{L}_{\ast }(v_{\infty }/\dot{M}{)}^2.$$

$$\mathbf{v}(r)=v_{\infty }{\left(1-\frac{R_{\ast }}{r}\right)}^{\beta },$$

$$n(r)=n_{\ast }{\left(\frac{R_{\ast }}{r}\right)}^{-2}{\left(1-\frac{R_{\ast }}{r}\right)}^{-1}.$$

$$n_{\ast }=\frac{\dot{M}}{4\pi \mu m_{\mathrm{p}}{R}_{\ast }^2{v}_{\infty }}=1.5\times 10^{10}\left(\frac{{\dot{M}}_{-7}}{v_{20}{R}_{35}^2}\right)\mathrm{c}{\mathrm{m}}^{-3},$$

$$r_1=[r_i^2+(r_o^2-r_i^2)u_r{]}^{1/2}.$$

$$DF\equiv \frac{{\hat{\mathbf{k}}}_{\mathbf{i}}\cdot {\mathbf{V}}_{\mathrm{Kep}}}{c},$$

$$D{F}^{\prime }=\frac{{\hat{\mathbf{k}}}_{\mathbf{f}}\cdot \mathbf{v}({\mathbf{r}}_1)}{c}.$$

$$\frac{{\lambda }_{\mathrm{f}}-{\lambda }_o}{{\lambda }_o}=\left(\frac{{\lambda }_{\mathrm{o}}}{{\lambda }_{\mathrm{i}}}\right)DF+D{F}^{\prime }.$$

$${\tau }_{\ast }(\lambda )=n_{\ast }{R}_{\ast }{\sigma }_{\mathrm{tot}}(\lambda ).$$

$${\tau }_{12}={\int }_{{\mathbf{r}}_1}^{{\mathbf{r}}_2}n(\mathbf{r}){\sigma }_{\mathrm{tot}}(\lambda )ds,$$

$$\tau (\rho )={\tau }_{\ast }{\int }_{\rho }^{\infty }\frac{d\rho }{(\rho -1)\sqrt{{\rho }^2-{\tilde{b}}^2}},$$

$$\rho =\tilde{b}\cosh \theta (\tau )-1.$$

$$\theta (\tau )$$

$$ϵ(\mathbf{r})=ϵ(\varphi )=A_j,\varphi \in ({\varphi }_j-\pi /8,{\varphi }_j+\pi /8).$$

$$FR(\varphi )=B_j,\varphi \in ({\varphi }_j-\pi /8,{\varphi }_j+\pi /8),$$

Full text

††pubyear: 2019††pagerange: Stellar Wind Accretion and Raman Scattered O VI Features in the Symbiotic Star AG Draconis–Stellar Wind Accretion and Raman Scattered O VI Features in the Symbiotic Star AG Draconis

1 Introduction

A symbiotic star is a wide binary system consisting of a white dwarf and a late type giant undergoing heavy mass loss (e.g., Kenyon, 1986). The cool giant component loses its material in the form of a slow stellar wind, some fraction of which is gravitationally captured by the white dwarf component giving rise to a variety of activities including X-ray emission, outbursts and prominent emission lines. The mass loss rate of the cool component in a typical symbiotic star has been proposed to be in the range $10^{-8}-10^{-6}\ \rm{M_{\odot}\ yr^{-1}}$ (e.g., Dupree, 1986).

AG Draconis is known to be a yellow symbiotic star having an early K type giant as mass donor (e.g. Leedjärv et al., 2016). The orbital period is known to be 550 days (e.g. Fekel et al., 2000). The light curve of AG Dra shows that it underwent major outbursts in intervals of 12-15 years with more frequent minor outbursts with a time scale of $\sim 1$ year (Hric et al., 2014). Sion et al. (2012) proposed the effective temperature of the hot component $T_{\rm eff}=80,000{\rm\ K}$ based on their analysis of spectra obtained with the Far Ultraviolet Spectroscopic Explorer (FUSE). Smith et al. (1996) investigated the chemical abundances of AG Dra to show that it is a metal poor star belonging to a halo population with overabundant heavy $s$-process elements.

The presence of an accretion disk is well established in cataclysmic variables, where the primary white dwarf is accreting material from the red dwarf companion filling its Roche lobe. However, in the case of symbiotic stars, there is lack of direct evidence indicating the presence of an accretion disk. Leibowitz et al. (1985) suggested the presence of an accretion disk in AG Dra based on the flat continuum. Robinson et al. (1994) proposed that an accretion disk is responsible for the double-peak profiles of H$\alpha$ observed in a number of symbiotic stars including AG Dra. However, they also noted the mismatch of H$\alpha$ profiles with the accretion disk model suggesting an alternative model, where H$\alpha$ is formed in the Strömgren sphere and suffers self-absorption in the outlying neutral region of the red giant wind.

A number of symbiotic stars show evidence of bipolar outflows that may be closely related to the formation of an accretion disk around the hot white dwarf component. The radio interferometric observations using the Multi-Element Linked Interferometer Network, (MERLIN) by Mikołajewska (2002) revealed a bipolar structure that is well aligned with the binary plane of AG Dra. A similar result was presented by Torbett & Campbell (1987). Despite the indication of a bipolar structure in AG Dra, observational verification of the presence of an accretion disk in this system still remains a difficult issue.

The formation of an accretion disk can be addressed through hydrodynamic calculations. Mastrodemos & Morris (1998) carried out Smoothed-Particle Hydrodynamics (SPH) calculations to propose that a stable accretion disk may form around a white dwarf that accretes material from a giant companion. The two dimensional hydrodynamical simulations performed by de Val-Borro et al. (2009) revealed that stable accretion disks can form with mass loss ranging $10^{-9}-10^{-6}{\rm\ M_{\odot}\ yr^{-1}}$. In particular, they noted that the accretion flow is characterized by eccentric streamlines and Keplerian velocity profiles.

Raman scattered O vi features at 6825 Å and 7082 Å that have been found only in bona fide symbiotic stars including AG Dra are unique tools to probe the accretion flow in symbiotic stars. These mysterious spectral features were identified by Schmid (1989), who proposed that they are formed through Raman scattering of O vi$\lambda\lambda$1032 and 1038 by atomic hydrogen. When a far UV O vi line photon is incident on a hydrogen atom in the ground state, the hydrogen atom is excited to one of infinitely many $p$ states. There are two channels of de-excitation available for the excited hydrogen atom. One is de-excitation into the ground $1s$ state, which corresponds to Rayleigh scattering. The other channel is Raman scattering, where the final de-excitation is made into the excited $2s$ state with an emission of an optical photon redward of H$\alpha$.

Many Raman scattered O vi features show double or triple peak profiles. Lee & Kang (2007) undertook profile analyses of the two symbiotic stars V1016 Cygni and HM Sagittae to show that the double peak profiles in these objects are well fit by adopting a Keplerian accretion disk around the hot white dwarf component (see also Heo & Lee, 2015; Heo et al., 2016). The fitted profiles are consistent with the O vi emission region in Keplerian motion with $\sim 50{\rm\ km\ s^{-1}}$ implying the physical dimension of the accretion disk $\leq 1{\rm\ au}$. In addition, they attributed the red enhanced asymmetric profiles to higher concentration of emitting material on the entering side of the accretion flow than on the opposite side.

In this paper, we present our high resolution spectrum of AG Dra showing Raman O vi features. In view of their great usefulness as a probe of the O vi emission region, we perform Monte Carlo profile analyses to propose that the Raman O vi profiles are consistent with the presence of an accretion disk in this binary system.

2 Spectroscopic Observation of AG Dra

2.1 CFHT Spectrum of AG Dra

We carried out high resolution spectroscopy of AG Dra using the Echelle SpectroPolarimetric Device for the Observation of Stars (ESPaDOnS) installed on the 3.6 m Canada-France-Hawaii Telescope (CFHT) in queued service observing mode on 2014 September 6. The total integration time was 2,000 seconds and the spectral resolution is 81,000. Data reduction was performed using the CFHT reduction pipeline Upena, which is based on the package Libre-ESpRIT developed by Donati et al. (1997). A standard Th-Ar lamp was used for wavelength calibration, according to which the pipeline wavelength should be shifted by $-0.18{\rm\ km\ s^{-1}}$ for heliocentric correction.

In Fig. 1, we show the part of the spectrum in the wavelength range between 6460 Å and 7140 Å, in which we find Raman scattered O vi features as well as H$\alpha$, two He i emission lines at 6678 Å and 7062 Å and many absorption features due to the giant atmosphere. The vertical axis shows the relative flux, where the normalization is made in such a way that the continuum level is of unit value.

The upper left and right panels of Fig. 2 show the parts of the spectrum around Raman O vi features centered at 6825 Å and 7082 Å, respectively. In these panels, the continuum is subtracted and the flux is rescaled so that the red peak of the Raman 6825 is of unit strength.

We note that the spectrum of AG Dra is plagued with many absorption features attributed to the atmosphere of the giant component. In order to assess the effect of the stellar atmospheric absorption lines on the overall profiles of the Raman O vi features, we invoke the BT-Settl model atmosphere (Allard, 2014). The results are presented in Fig. 2. In this figure, the top spectra are the blow-up version around the Raman features of the CFHT data shown in Fig. 1. The second and third spectra show the BT-Settl model atmospheres of a yellow giant with $T_{\rm eff}=4,300$ and $3800{\rm\ K}$, respectively. As is found in the figure, the effect of molecular absorption features is more serious in the case of $T_{\rm eff}=3,800{\rm\ K}$ than $T_{\rm eff}=4,300{\rm\ K}$.

In the bottom we show the resultant spectra obtained by subtracting the model atmosphere with $T_{\rm eff}=4,300{\rm\ K}$ from the observational data. The Raman 6825 feature is quite strong and affected by relatively weak absorption features, exhibiting an overall smooth profile. In contrast, in the case of the Raman 7082 feature, the subtraction is reasonably good except a number of absorption/emission features in the Doppler factor region $-25{\rm\ km\ s^{-1}}\leq\Delta V<0{\rm\ km\ s^{-1}}$, which can be attributed to imperfect removal of continuum. The imperfection of subtraction makes it very difficult to infer the exact amount of the Raman 7082 flux that is absorbed.

2.2 Profile Comparison

Raman and Rayleigh scattering of a far UV photon with an atomic electron is described by time dependent second order perturbation theory (e.g. Schmid, 1989; Lee & Lee, 1997). The cross sections for O vi$\lambda\lambda$1032 and 1038 are $\sigma_{\rm Ram}(1032)=7.5\sigma_{\rm T}$ and $\sigma_{\rm Ram}(1038)=2.5\sigma_{\rm T}$, respectively, where $\sigma_{\rm T}=0.665\times 10^{-24}{\rm\ cm^{2}}$ is the Thomson scattering cross section (e.g., Schmid, 1989; Nussabumer et al., 1989; Lee et al., 2016). The Rayleigh scattering cross sections are found to be $\sigma_{\rm Ray}(1032)=34\sigma_{\rm T}$ and $\sigma_{\rm Ray}(1038)=6.6\sigma_{\rm T}$.

The energy conservation principle sets a relation between the wavelength $\lambda_{i}$ of an incident photon and that $\lambda_{o}$ of Raman scattered one given by

[TABLE]

where $\lambda_{\alpha}$ is the wavelength of Ly$\alpha$. According to Schmid (1989), the fiducial atomic line centers in air for Raman scattered O vi are 6825.44 Å and 7082.40 Å. Differentiation of Eq. (1) yields

[TABLE]

from which it may be noted that the Raman O vi features will be broadened by a factor $(\lambda_{o}/\lambda_{i})\sim 6$ with respect to the incident far UV O vi. This relation dictates that Raman O vi profiles reflect relative motion between O vi and H i and that they are almost independent of the observer’s line of sight. In view of this property, no consideration of transforming to the binary center of mass frame is made in this work.

In the spectrum shown in Fig.2, the Raman 6825 and 7082 features are found to lie in the ranges 6810 Å $<\lambda<$ 6840 Å and 7068 Å $<\lambda<7095$ Å, respectively, which translates into the full width at zero intensity (FWZI) of $200{\rm\ km\ s^{-1}}$ with respect to the neutral scattering region around the giant. If we interpret this velocity as the twice of the Keplerian speed at the innermost O vi emission region, then it implies that

[TABLE]

where $r_{i}$ is the distance of the innermost O vi emission nebula from the white dwarf and $M_{\rm WD}$ is the mass of the white dwarf.

In Fig. 3, we overplot the two Raman spectral features in the Doppler factor space in order to make a comparison of the two profiles. The two upper panels show our double Gaussian fit to the profiles, where the sum of two Gaussian functions are represented by gray lines. Taking the emission line He i$\lambda$7065 with the atomic line center at 7065.190 Å as velocity reference, we determine the centers of the two Raman O vi features, which are shown by a vertical arrow in the top two panels of Fig. 3. The upper horizontal axis shows the Doppler factors computed in the parent far UV O vi lines. Both Raman 6825 and 7082 features show a main peak shifted redward from the profile center by a Doppler factor of $\sim 25{\rm\ km\ s^{-1}}$.

The double Gaussian decomposition is made first for the Raman 6825 feature. We find that the center separation between the two Gaussian components is measured to be $\sim 57{\rm\ km\ s^{-1}}$. The full widths at half maximum of the blue and red Gaussian functions are $\Delta\lambda=8.15$ Å which corresponds to the O vi velocity width of $54{\rm\ km\ s^{-1}}$. The red Gaussian is 1.19 times stronger than the blue counterpart. We divide the two Gaussian functions by 2.5 and apply to the Raman 7082 feature.

If the double component profiles are attributed to the accretion disk around the white dwarf, the outer radius of the accretion disk is roughly $r_{o}\sim 0.3{\rm\ au}$ with the adopted white dwarf mass of $M_{\rm WD}=0.5{\rm\ M_{\odot}}$ (Mikołajewska et al., 1995). It is also noted that the location of Doppler factor zero appears to be shifted blueward with respect to the center of zero intensities for both the Raman features, implying an additional contribution to the red part.

For a clear comparison, we multiply the flux of Raman 7082 feature by a factor 2.5 in the bottom right panel. It is found that the two profiles coincide within the uncertainty. Whereas the red parts of the two profiles match excellently, the Raman 7082 appears slightly weaker than the 6825 counterpart in the velocity range $-25{\rm\ km\ s^{-1}}<\Delta V<0{\rm\ km\ s^{-1}}$. the absorption features near $\Delta V\sim-20{\rm\ km\ s^{-1}}$ and insufficient data quality prevent a definite conclusion .

3 Formation of Raman Scattered O VI at 6825 Å and 7082 Å

3.1 Ionization Front

The slow stellar wind region in the vicinity of the giant is illuminated and photoionized by the far UV radiation from the hot white dwarf. The ionization front is formed where the photoionization rate balances the recombination rate. The former is closely associated with the mass accretion rate onto the white dwarf component whereas the latter is mainly determined by the mass loss rate $\dot{M}$ of the giant. Seaquist et al. (1984) introduced an ionization parameter $X$ to characterize the ionization front, which is

[TABLE]

Here, $L_{*},\mu,\alpha,m_{\rm p}$ and $a$ are the hydrogen ionizing luminosity, the mean molecular weight, the recombination coefficient, the proton mass and the binary separation, respectively.

In this work, we adopt a simple spherically symmetric stellar wind model so called a $\beta$-law given by

[TABLE]

where $v_{\infty}$ is the wind terminal speed, $R_{*}$ is the distance of the wind launching place. In this work, we set $R_{*}=35R_{\odot}$ which Hric et al. (2014) proposed as the radius of the giant component of AG Dra.

We also choose $\beta=1$ for simplicity. The density $n(r)$ distribution corresponding to the $\beta$-law is given by

[TABLE]

Here, the characteristic number density is

[TABLE]

where we set $\mu=1.4$ and the scaled parameters $\dot{M}_{-7}=\dot{M}/(10^{-7}{\rm\ M_{\odot}\ yr^{-1}})$, $v_{20}=v_{\infty}/(20{\rm\ km\ s^{-1}})$ and $R_{35}=R_{*}/(35R_{\odot})$.

In this stellar wind model, neutral hydrogen column density $N_{HI}$ diverges along lines of sight from the white dwarf toward the giant. On the other hand, if the mass loss rate is moderate, complete photoionization results along lines of sight with large impact factors from the giant. Therefore, one finds that the ionization front is well approximated by a hyperboloid when the accretion luminosity is sufficient to fully ionize except the neighboring region of the giant component (e.g., Lee & Kang, 2007; Sekerás & Skopal, 2015). In this work, we consider mainly the case where the half-opening angle $\theta_{o}=90^{\circ}$ of the ionization front with respect to the white dwarf. Fig. 4 shows a schematic illustration adopted in this work.

Taking the mass $M_{\rm G}=1.5{\rm\ M_{\odot}}$ and radius $R_{\rm G}=35{\rm\ R_{\odot}}$, we show the position and size of the giant component to the scale in Fig. 4 (Skopal, 2005; Hric et al., 2014). In Table 1, we summarize the model parameters of AG Dra adopted in this work. The inner Lagrange point is located at 140 $R_{\odot}$ from the white dwarf, allowing an accretion disk with a Keplerian speed $\sim 30{\rm\ km\ s^{-1}}$ to be accommodated inside the Roche lobe of the white dwarf. The $x-$axis is chosen to coincide with the binary axis. The giant and the white dwarf are located at $x=400\ R_{\odot}$ and at the center, respectively. We put an additional H i region along $x-$axis to mimic the bipolar structure of AG Dra.

3.2 Monte Carlo Approach

In this subsection, we briefly describe our Monte Carlo code to simulate Raman O vi profiles. The Monte Carlo simulation starts with generation of an O vi line photon in the emission region around the white dwarf. In view of the FWZI and the velocity separation between the two Gaussian components of the Raman features, we identify the O vi emission region as the part of a Keplerian disk around the white dwarf with the inner and outer radii denoted by $r_{i}$ and $r_{o}$, where these two parameters are described in Subsection 2.2. Using a uniform random number $u_{r}$ between 0 and 1, an O VI line photon is generated initially at a position with radial coordinate $r_{1}$ by the prescription

[TABLE]

As a check of the code, we consider a point-like and a ring-like O vi emission source by setting $r_{i}=r_{o}=0$ and $r_{i}=r_{o}$ to a positive constant, respectively. For these cases, the O vi emissivity is assumed to be isotropic. However, the observed red enhanced profiles require anisotropic O vi emissivity, which is described in subsection 5,1.

The initial unit wavevector ${\bf\hat{k}_{i}}$ is chosen isotropically. To the initial photon at ${\bf r}_{1}$ with ${\bf\hat{k}_{i}}$, we assign a Doppler factor $DF$ given by

[TABLE]

where ${\bf V}_{\rm Kep}$ is the Keplerian velocity vector associated with the O vi line emitter. In the case of a point-like emission region, ${\rm V}_{\rm Kep}$ is set to be zero.

We trace photons entering the H i region, where Rayleigh or Raman scattering may take place. According to the total scattering cross section $\sigma_{\rm tot}=\sigma_{\rm Ray}+\sigma_{\rm Ram}$, a next scattering site is determined in a probabilistic way. The scattering type is determined according to the branching ratios $b_{\rm Ram}=\sigma_{\rm Ram}/\sigma_{\rm tot}$ and $b_{\rm Ray}=1-b_{\rm Ram}$. If the scattering is Raman, then the Raman photon is presumed to escape from the region to reach the observer. Otherwise, we generate a new unit wavevector ${\bf\hat{k}_{f}}$ along which we find a new scattering site corresponding to a scattering optical depth $\tau$. A new Doppler factor is assigned to the new ${\bf\hat{k}_{f}}$ and the velocity of the scatterer, which is given by

[TABLE]

It may occur that no such scattering site is found when $\tau$ exceeds the optical depth corresponding to the infinity, in which case the O vi line photon escapes through Rayleigh scattering. A scattering loop ends either with Raman escape or with Rayleigh escape. For a Raman-escaping photon, the final wavelength $\lambda_{\rm f}$ is given by

[TABLE]

The probabilistic determination of a next scattering site ${\bf r}_{2}$ from a starting position ${\bf r}_{1}$ and a unit wavevector ${\bf\hat{k}}$ involves the inversion of the physical distance $s$ in terms of a total scattering optical depth $\tau$, which we describe now.

The characteristic total scattering optical depth $\tau_{*}$ is defined as

[TABLE]

The numerical values for O vi are $\tau_{*}(1032)=1.02$ and $\tau_{*}(1038)=0.22$. Given a unit wavevector ${\bf\hat{k}}$, the scattering optical depth between two positions ${\bf r}_{1}$ and ${\bf r}_{2}$ is given by the line integral

[TABLE]

where the parameter $s$ measures the physical distance from ${\bf r}_{1}$ to ${\bf r}_{2}$, that is, ${\bf r_{2}}={\bf r}_{1}+s{\bf\hat{k}}$ . We introduce the impact parameter $b$ that measures the perpendicular distance from the giant to the photon path. In terms of $\tau_{*}$, the integral of interest can be written as

[TABLE]

where the dimensionless parameters are $\rho=r/R_{*}$ and $\tilde{b}=b/R_{*}$ (Lee & Lee, 1997).

The inverse relation for the physical distance corresponding to a scattering optical depth $\tau$ is

[TABLE]

Here, the parameter $\theta(\tau)$ is

[TABLE]

when $\tilde{b}<1$. For $\tilde{b}>1$, the corresponding formula is

[TABLE]

The red enhanced profiles shown in Fig. 2 imply that the accretion disk is asymmetric in the azimuthal direction. Hydrodynamic simulations show that the accretion flow tends to converge on the entrance side as the wind material from the giant is captured by the white dwarf. However, on the opposite side the flow tends to diverge so that the O vi density diminishes leading to weaker O vi emission than the entrance side (e.g. de Val-Borro et al., 2009).

4 Monte Carlo Radiative Transfer and Basic Results

4.1 Raman Conversion Efficiency

In this subsection, we compute Raman conversion efficiency (RCE) for O vi$\lambda\lambda$ 1032 and 1038 as a function of the mass loss rate of the giant and the half-opening angle of the ionization front. Here, we define RCE as the number ratio of far UV incident photons generated in the O vi emission region and Raman scattered photons formed in the H i region. With the availability of far UV spectroscopy, Raman conversion efficiency is a good measure of the mass loss rate of the giant component.

Birriel et al. (2000) proposed RCE of 15 percent for O vi$\lambda$1032 from their near simultaneous observations of far UV using the Hopkins Ultraviolet Telescope and optical spectroscopy. Schmid et al. (1999) used the data obtained with the Orbiting & Retrievable Far and Extreme Ultraviolet Spectrometer to deduce much higher efficiency of 50 percent with quoted error of $\sim$ 20 percent. Birriel et al. (2000) discussed possible origins for the discrepancy including the simultaneity of far UV and optical spectroscopy and interstellar extinction.

Schmid (1996) carried out Monte Carlo simulations to study basic properties of Raman scattering of O vi with atomic hydrogen and estimate the RCE of O vi doublets in symbiotic stars. Lee et al. (2016) investigated the Raman conversion efficiencies of O vi$\lambda\lambda$ 1032 and 1038 for a few simple scattering geometries to find that the flux ratio $F(6825)/F(7082)$ may vary in the range from 1 to 6.

In Fig. 5, we show the RCE for O vi$\lambda\lambda$ 1032, 1038 and the flux ratio $F(6825)/F(7082)$. The binary separation is set to $a=400{\rm\ R_{\odot}}$ and the mass loss rate of the giant is taken to be in the range $10^{-8}-10^{-6}{\rm\ M_{\odot}\ yr^{-1}}$. We present our results for the two values $\theta_{o}=60^{\circ}$ and $90^{\circ}$ of the half-opening angle of the ionization front.

The overall RCE increases monotonically as $\dot{M}$ increases. The maximum Raman conversion efficiency of O vi$\lambda$1032 is found to be 0.23 and 0.11 for $\theta_{o}=90^{\circ}$ and $\theta_{o}=60^{\circ}$, respectively. For O vi$\lambda$1038, smaller values of 0.13 and 0.07 are obtained. When the mass loss rate is high, RCE becomes very sensitive to the half-opening angle. However, when the mass loss rate is low $\leq 10^{-7}{\rm\ M_{\odot}\ yr^{-1}}$, RCE is mainly determined by $\dot{M}$ and almost independent of $\theta_{o}$. In the case of low mass loss rates, most of Raman photons are formed near the giant, where the optical depth diverges irrespective of $\dot{M}$ in our model. It is interesting to note that the RCE of 15 percent of O vi$\lambda$1032 proposed by Birriel et al. (2000) is consistent with a mass loss rate $\dot{M}\simeq 4\times 10^{-7}{\rm\ M_{\odot}\ yr^{-1}}$ with $\theta_{o}=90^{\circ}$.

4.2 Pure $\beta$-Law for the Giant Stellar Wind

In Fig. 6, we show the profiles of Raman scattered O vi features formed in a stellar wind with $\beta=1$ for various values of the mass loss rate $\dot{M}$ and the giant wind terminal speed $v_{\infty}$. In this figure, the O vi emission region is taken to be point-like and monochromatic. The vertical arrow indicates the atomic Raman line center, that is the wavelength a line center photon would acquire if it is Raman scattered photons in a static H i region.

We consider three values of the terminal wind velocity $v_{\infty}=10,15$ and $20{\rm\ km\ s^{-1}}$. As the wind terminal speed is increased, the profile width increases. At a specified neutral region around the giant, the density is proportional to $\dot{M}$ and inversely proportional to $v_{\infty}$ so that the Raman flux decreases as $v_{\infty}$ increases given mass loss rate.

We note that the profile peak is found significantly redward of the Raman line center. This is also found in the work of Lee & Lee (1997), who also presented the contour of the Raman scattering optical depth in their Fig. 2. The contour of the Raman scattering optical depth of unity takes approximately a hyperboloidal shape, which shows that O vi line photons are scattered more efficiently in the stellar wind region receding from the white dwarf. The red peaks are found at $\Delta\lambda=+1.3,1.9$ and 2.6 Å from the atomic line center for Raman scattered 6825 feature for $v_{\infty}=10,15$ and $20{\rm\ km\ s^{-1}}$, respectively. These red peaks turn out to correspond to $v=0.75v_{\infty}$.

In addition, there appears a blue shoulder in the cases where $\dot{M}\leq 1\times 10^{-6}{\rm\ M_{\odot}\ yr^{-1}}$. This is due to contribution from the part of the giant stellar wind facing the white dwarf, where Raman scattering optical depth increases without limit for lines of sight hitting the giant surface. For $\dot{M}>10^{-6}{\rm\ M_{\odot}}$, Raman scattering occurs quite far from the giant surface, where the Raman scattering optical depth is moderate. Therefore, the neutral region near the giant surface facing the white dwarf contributes little, and no shoulder features are found in the bottom panels of Fig. 6.

Fig. 7 shows the profiles and relative fluxes of Raman 6825 and 7082 features where the O vi emission region is assumed to form a uniform ring around the white dwarf. We consider four values of the Keplerian rotation velocity ${\bf V}_{\rm Kep}$ of the O vi emission region, which are 40, 50 ,60 and $70{\rm\ km\ s^{-1}}$. In this figure, we fix the mass loss rate of the giant $\dot{M}=4\times 10^{-7}\ {\rm M_{\odot}\ yr^{-1}}$, the half-opening angle of the ionization front $\theta_{o}=90^{\circ}$ and the binary separation of $400{\rm\ R_{\odot}}$. The terminal velocity of the giant wind is set to $v_{\infty}=10{\rm\ km\ s^{-1}}$.

The two main peaks correspond to the rotation speed of the O vi emission ring. As in the case of Fig. 6, the red peaks are slightly stronger than the blue counterparts. A close look reveals that the full width at zero intensity is slightly larger than the peak separation due to convolving effect of the neutral wind around the giant. We notice that there appears slight excess in the relative flux near the blue and red ends, which is also attributed to the effect of convolution. It should be noted that the overall profile is mainly determined by the kinematics of the emission region if the giant wind speed is small compared to the profile width. The conspicuous red enhancement apparent in the observational data strongly implies that the O vi emission region is significantly azimuthally asymmetric.

5 Raman O VI Profiles and Comparisons with Observation

5.1 Raman Flux Ratio and the Mass Loss Rate

In this subsection, we consider an annular O vi emission region depicted schematically in Fig.4. For simplicity, we assume that the O vi local emissivity $\epsilon({\bf r})$ is a function of only azimuthal angle $\phi$ measured from the line toward the giant component. We further assume that the flux ratio $F(1032)/F(1038)$ is constant throughout the annular O vi region.

Dividing the entire azimuthal angle domain into 8 equal intervals, each being centered at $\phi_{j}=j\pi/4$, $(j=0,1,2,\dots,7),$ we try piecewise step functions taking a constant value in each interval. That is,

[TABLE]

Once we find a set of parameters $A_{j}$ fitting the observed Raman 6825 data, the same set of parameters are used to fit the Raman 7082 feature. In the center panel of Fig. 8, we show our best fitting profiles for Raman 6825 for three different values of $\dot{M}=1\times 10^{-7},\ 2\times 10^{-7}$ and $4\times 10^{-7}{\rm\ M_{\odot}\ yr^{-1}}$. The criterion for best fit adopted in this subsection is the least chi-square in the wavelength range between 6815 Å and 6830 Å that encompasses the blue shoulder and the red main peak.

The best fitting values of $A_{j}$ are listed in Table 2 and we also illustrate $\epsilon(\phi)$ in the left panel of Fig. 8. In the figure, it is noted that the O vi emissivity peaks at around $\phi=90^{\circ}$, where it is 2.5 times larger than at $\phi=270^{\circ}$.

In the right panel of Fig. 8, we present three profiles simulated for Raman O vi$\lambda$7082. Our CFHT spectra are overplotted to illustrate the quality of the fit in the center and right panels. The main effect of the mass loss rate $\dot{M}$ is found in the flux ratio of $F(6825)/F(7082)$, where the flux ratio tends to decrease as $\dot{M}$ increases.

The simulated profiles of Raman 6825 features fail to fit the red part, which is also clearly linked to the blueward shift of the atomic line center. It appears that the observed Raman 7082 feature is best fit with $\dot{M}=2\times 10^{-7}{\rm\ M_{\odot}\ yr^{-1}}$. However, the difference between simulated profiles for various values of $\dot{M}$ is not large compared to the uncertainty of the observational data around 7082 Å, from which it is very difficult to single out a definite value of $\dot{M}$. Further studies including polarimetric observations and spectroscopic monitoring will provide additional constraints on the mass loss rate and mass transfer processes.

5.2 Bipolar Structure

The red excess apparent in the Raman O vi profiles indicates the operation of Raman scattering between O vi and H i components that are receding from each other. Torbett & Campbell (1987) presented their radio 6 cm observation of AG Dra to reveal a bipolar structure present in this system. It appears that bipolar structures are commonly found in symbiotic stars, where collimated outflows are formed in association with the mass accretion on to the white dwarf component (e.g., Angeloni et al., 2011; Skopal et al., 2013). In this subsection, we further consider the presence of an H i component moving away from the binary system, from which an additional contribution is made to the red wing flux of Raman scattered O vi.

The asymmetric wind accretion model provides a reasonable fit to the observational data with the exception that the red wing parts are observed stronger than model fluxes. In order to augment the deficit in the red wing parts for both simulated Raman features, it is necessary to introduce an additional H i region that is moving away from the two stars in the direction perpendicular to the orbital plane.

A significant fraction of planetary nebulae exhibit bipolar morphology, of which the physical origin is still controversial (e.g., Sahai et al., 2007). Possible mechanisms include the binarity of the central star system and magnetic fields. Akras et al. (2017) reported the detection of near IR emission from ${\rm H}_{2}$ in the bipolar planetary nebula K4-47, which delineates the low-ionization structures. They invoked the interaction of collimated flows or bullets that are moving away from the central star along the bipolar directions. Interestingly enough, the receding velocity of the low-ionization structures is typically of order $100{\rm\ km\ s^{-1}}$.

In this subsection, we introduce bipolar H i regions moving away in the directions perpendicular to the orbital planes. The left panel of Fig.9 shows a schematic illustration of the scattering geometry with an addition of bipolar H i regions. The half-opening angle $\theta_{o}$ of the polar H i region is chosen to be $45^{\circ}$ with a receding speed of $70{\rm\ km\ s^{-1}}$. The H i column density is set to $N_{HI}=10^{21}{\rm\ cm^{-2}}$.

In the center and right panels of Fig. 9, we show our simulated profiles with an additional contribution from a receding H i region. The additional fluxes are shown in gray solid lines, constituting about 20 percent of the total flux in both the Raman 6825 and 7082 features. One notable feature is that the additional Raman fluxes exhibit a single peak profile. This is due to the fact that the velocity component parallel to the orbital plane does not contribute to the Doppler factor along the perpendicular direction. With this additional component, substantial improvement is achieved in the quality of fit in the red parts.

5.3 Local Variation of O VI Flux Ratio in the Accretion Flow

In Subsection 2.2, we discussed the possibility that the blue part of the Raman 7082 feature may be relatively more suppressed than in the Raman 6825 counterpart. The disparity of the two Raman O vi profiles in symbiotic stars including V1016 Cyg has been pointed out by a number of researchers (e.g., Schmid et al., 1999). The disparate profiles are attributed to locally varying flux ratios of $F(1032)/F(1038)$ in addition to the asymmetric accretion flow around the hot white dwarf component (Heo & Lee, 2015).

Resonance doublet lines including C iv$\lambda\lambda$1548, 1551 and O vi$\lambda\lambda$1032, 1038 arise from transitions between $S_{1/2}-P_{3/2,1/2}$. The short wavelength component $S_{1/2}-P_{3/2}$ has twice larger statistical weight than the long wavelength component, leading to the flux ratio of 2 in an optically thin nebula. Far UV spectroscopy of symbiotic stars and planetary nebulae reveals that resonance doublet lines rarely show the flux ratio of 2. Instead their flux ratios are observed to vary from one to two (Feibelman, 1983). A similar result was reported by Vassiliadis et al. (1996), who investigated a dozen planetary nebulae in the Magellanic Clouds using the Hubble Space Telescope.

In the accretion flow model shown in Fig. 4, we pointed out that the entering side is denser than the opposite side. Nonuniform O vi density distribution may lead to flux ratio $F(1032)/F(1038)$ that vary locally in the O vi region. As an exemplary case, we assume that flux ratio $FR=F(1032)/F(1038)$ is given as a function of only azimuthal angle $\phi$.

As in the previous section for $\epsilon(\phi)$, we try piecewise step functions for $FR(\phi)$ given as

[TABLE]

where $B_{j}\ (j=0,1,2,\dots,7)$ are constants between 1 and 2. In our model, $FR(\phi)$ is close to 2 in the blue region with $180^{\circ}<\phi<360^{\circ}$ whereas it is near one in the red region with $0<\phi<180^{\circ}$.

In the left panel of Fig. 10, we illustrate $FR(\phi)$ and the values of $B_{j}$ are listed in Table 3. In the center and right panels, we show our Monte Carlo profiles, which are overplotted on the observational data. Here, the dashed line represents the profiles shown in Fig. 9 for comparison, in which the flux ratio is uniform throughout the O vi region. As is expected, the resultant Raman 7082 profile shows significant deficit in the region $-30{\rm\ km\ s^{-1}}\leq\Delta V\leq 0$ compared to the best fit profile obtained from the uniform flux ratio model, achieving enhancement in quality of profile fitting. Future high resolution spectroscopy with higher signal to noise ratio is expected to resolve the issue of the profile disparity and the local variation of the flux ratio.

6 Summary and Discussion

In this work, we present our high resolution spectrum around Raman scattered O vi bands of the yellow symbiotic star AG Dra obtained with the CFHT. Both the Raman 6825 and 7082 features exhibit asymmetric double component profiles. We find that in the reference frame of He i$\lambda$7065, the atomic line centers of both Raman 6825 and 7082 features are shifted blueward, implying an additional contribution from components that recede from the two stars.

We have carried out profile analyses by performing Monte Carlo simulations. The best fit model includes an accretion disk with a size $<0.3$ au around the white dwarf, a slow stellar wind from the giant with a mass loss rate $\sim 2\times 10^{-7}{\rm\ M_{\odot}}\ {\rm yr^{-1}}$ and neutral regions moving away from the binary orbital plane in the perpendicular directions with a speed $\sim 70{\rm\ km\ s^{-1}}$.

The observed line profiles of Raman O vi 6825 and 7082 features exhibit overall redward shift with respect to their expected atomic line centers. With this in mind we introduced bipolar neutral regions that recede from the binary system along the direction perpendicular to the orbital plane. The bipolar neutral component provides an additional contribution to the red wing part of Raman O vi features, offering significant improvement in overall quality of profile fitting. This additional component is naturally inferred from radio observations indicative of the bipolar morphology of AG Dra (Torbett & Campbell, 1987).

Lee et al. (2016) pointed out that the Raman conversion efficiency is also a sensitive function of wavelength so that the efficiency of a blue 1032 photon is slightly larger than that of a red 1032 photon. A similar effect is also expected of O vi$\lambda$1038. Our current Monte Carlo code adopts a set of fixed values for Rayleigh and Raman cross sections for O vi$\lambda$1032 and another set for O vi$\lambda$1038, neglecting the slight difference of the cross section depending on the wavelength. A careful incorporation of this effect into the simulation may strengthen slightly the blue parts of Raman O vi features, leading to further slight distortion of the resultant profiles.

It should be noted that a number of factors may affect Raman O vi including the variation of $\dot{M}$ due to pulsational instability of the giant. Mass transfer may also be subject to instability associated with the accretion disk. It is also quite likely that the part of the H i region illuminated by O vi line radiation and visible to the observer will be dependent on the binary orbital phase. Our current Monte Carlo code does not include these subtle points, which will be considered in the future works.

Raman-scattered features consist of purely scattered photons without any dilution from unpolarized incident radiation. Due to this inelasity and the dipole nature of Raman scattering, the 6825 and 7082 features may develop strong linear polarization in the direction perpendicular to the scattering plane. Harries & Howarth (1996) carried out spectropolarimetry of many symbiotic stars to show that Raman O vi features are generally polarized with large degree of polarization. They also showed that the position angle in the red wing parts differs by an amount of $90^{\circ}$ from the main parts for both Raman O vi features, which is consistent with the bipolar structure inherent to symbiotic stars (e.g., Lee & Park, 1999).

According to Schmid & Schild (1997), AG Dra exhibits rotation of polarization angle by an amount of 45 degrees in the Raman O vi features. That is, the incident direction for red wing photons may differ by the same angle to the line connecting the white dwarf and the giant. Our model strongly implies that the position angle differs by $90^{\circ}$. A more sophisticated scattering geometry including additional H i regions receding in oblique directions, which is deferred to future works. Further high resolution spectropolarimetry will shed much light on the scattering geometry and mass loss process of AG Dra.

Acknowledgements

We are grateful to the anonymous referee for constructive comments that helped improve the paper. This research was supported by the Korea Astronomy and Space Science Institute under the R&D program(Project No. 2015-1-320-18) supervised by the Ministry of Science, ICT and Future Planning. This work was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MSIT) (No. NRF-2018R1D1A1B07043944). This work was also supported by K-GMT Science Program (PID: 14BK002) funded through Korea GMT Project operated by Korea Astronomy and Space Science Institute. RA acknowledges financial support from the DIDULS Regular PR17142 by Universidad de La Serena.