Spectral Modeling of the Supersoft X-ray Source CAL87 based on Radiative Transfer Codes

CAL87 (Long et al., 1981) is a reasonably bright source located in the Large Magellanic Cloud (LMC) at a distance of 48.1 $\pm$ 2.9 kpc (Macri et al., 2006). It is an eclipsing binary with an orbital period of 10.6 hr (Pakull et al., 1988; Callanan et al., 1989; Cowley et al., 1990; Schmidtke et al., 1993), suggesting a high inclination angle (Schandl et al., 1997). The X-ray emission is considered to come from the extended accretion disk corona to account for the shallow and long X-ray eclipse due to the companion star (Schmidtke et al., 1993; Asai et al., 1998; Ebisawa et al., 2001) and for the lack of changes in the X-ray emission line intensity in and out of the eclipse (Ribeiro et al., 2014).

Ebisawa et al. (2001) proposed a picture, in which the X-ray emission from the WD atmosphere is permanently blocked, but its electron-scattered emission in the corona is observed. The observed O VII and O VIII edges, recognized with the CCD spectral resolution, are footprints of the incident WD atmosphere emission. The corona radius is estimated as $\sim$5 $\times 10^{10}$ cm, which is $\sim$0.4 times the orbital separation, and the inclination angle is $\sim$73 degree to account for the observed X-ray light curve. Emission lines should be observed from the corona, which were not recognized in the CCD spectra, but indeed found later with the improved energy resolution of X-ray grating spectrometers (Orio et al., 2004; Greiner et al., 2004; Ebisawa et al., 2010; Ribeiro et al., 2014). We follow the picture thus established.

CAL87 was observed with the XMM-Newton observatory (Jansen et al., 2001) on 2003-04-18 for 21.8 hours covering two orbital cycles (the observation number 0153250101). We used the Reflection Grating Spectrometer (RGS; Den Herder et al. 2001) to exploit its high-resolution X-ray spectra as well as the Metal Oxide Semiconductor (MOS; Turner et al. 2001) type detector of the European Photon Imaging Camera (EPIC) to constrain the continuum level.

The data were retrieved from the archive and reduced using the Science Analysis System (SAS) software version 13.5.0 with the calibration files available as of December 5, 2014. For the spectral fitting, we used Xspec version 12.12.1. We used the standard processing using the rgsproc task to extract the dispersion of the first order in the RGS data both for the source and background spectra. We combined the RGS1 and RGS2 spectra using rgscombine task to improve the statistics for all the data. The MOS data were processed with the standard processing using the emproc task. We selected the fiducial events with a PATTERN of 0–12 and removed those affected by the background flares. Then, we extracted the source and background events respectively from a circle centered at the source with a radius $40^{\prime\prime}$ and an annulus of the inner and outer radius of $75^{\prime\prime}$ and $100^{\prime\prime}$. Relative normalization between MOS1 and MOS2 was assumed identical, while that between MOS and RGS was determined in the spectral fitting. The spectra are averaged over the exposure time because the spectral shape does not change in and out of the eclipse (Ebisawa et al., 2010).

The X-ray spectra of the two MOS and the combined RGS data are shown in Figure 1. We characterize them by constructing a phenomenological spectral model for later use in physically-motivated models. For all the models in this work, we included the attenuation by the Galactic interstellar medium (ISM) in the direction of CAL87 and the ISM in the LMC using two photoelectric absorption models (tbabs and tbvarabs: Wilms et al., 2000). The column density of the former ($N_{\mathrm{H}}^{\mathrm{Gal}}$) was fixed at 7.58$\times$10²⁰ cm^-2 (Dickey & Lockman, 1990) and that of the latter ($N_{\mathrm{H}}^{\mathrm{LMC}}$) was derived from the fitting. The metal abundance for the LMC absorption was fixed to be half of the Galactic ISM value (Russell & Dopita, 1992; Welty et al., 1999).

The continuum model was constrained using the MOS spectra with the richer statistics (Fig. 1 a). The spectra are characterized by a sharp cut-off at 0.7–0.9 keV despite the monotonically increasing effective area of the telescope. We thus used a single blackbody model with a photoelectric absorption by the O VII and O VIII K shell edges at 0.74 and 0.87 keV, respectively. The same model was used in the previous work using X-ray CCD spectra (Asai et al., 1998; Ebisawa et al., 2001). The bolometric luminosity corrected for the interstellar absorption ($L_{\rm obs}$) is $1.4\times 10^{37}$ erg s^-1 assuming the spherical emission.

Upon the best-fit continuum model, we added line components by using both the MOS and RGS (Fig. 1 b). A handful of emission lines were identified in the RGS spectrum. We used Gaussian models for each line. From the best-fit energies, we identified emission lines of O VIII Ly$\alpha$ and Ly$\beta$, O VII He$\alpha$, N VII Ly$\alpha$ and Ly$\beta$, and four Fe XVII lines for the transition to the ground state (2s)²(2p)⁶ from the excited states shown in Table 2. All of them are electric dipole transitions of a large oscillator strength among the Fe XVII lines (Chen et al., 2003). These identifications are reported in Ebisawa et al. (2010).

Based on the line identifications, we fitted individual line complexes in a narrow (20–30 eV) energy range. For each of the Ly$\alpha$ and Ly$\beta$ line complexes, we included two lines (${}^{2}P_{1/2}$ and ${}^{2}P_{3/2}$). The relative energy and intensity ratio between the two lines were fixed. The line intensity, energy shift, and broadening were fitted collectively. For the He$\alpha$ line complex, we modeled with three lines of resonance ($r$), inter-combination ($i$), and forbidden ($f$). For the Fe XVII line complex at 0.73–0.74 keV, we modeled with two lines (3F and 3G; Table 2). For the Fe XVII line complex at 0.81–0.83 keV, we modeled with two lines (3C and 3D; Table 2). The relative energy of each line was fixed. The line energy shift and broadening were fitted collectively, while the line intensity was fitted individually.

The result of the line complex fitting is shown in Figure 2. For the line shift and broadening, most complexes exhibit a redward shift with a mean of $\sim$900 km s^-1 and broadening of $\sim$400 km s^-1, which have been reported in Orio et al. (2004) and Ribeiro et al. (2014). The line shift of N VII Ly$\beta$ deviates from the trend, which may be biased by the RGS chip gap eliminating a large fraction of the line profile.

We assume a static and spherically-symmetric geometry (Fig. 3), in which a point-like source (WD) at the center is surrounded by the photoionized plasma ($=$corona). The WD has an effective temperature of $T_{\mathrm{eff}}$, a surface gravity of $g$, and a luminosity of $L_{\mathrm{WD}}$. It emits incident photons that photo-ionize the corona. It is assumed that the observer does not directly observe the incident emission from the WD being fully blocked by the disk but does observe the diffuse emission from the corona and the scattered emission originating in the WD. The inner and outer radii of the corona are $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$, respectively. The electron density profile is assumed to be $n_{\mathrm{e}}(r)=n_{\mathrm{e,in}}(r/r_{\mathrm{in}})^{-2}$. The metal abundance in the corona was fixed to 0.5 times the solar value (Wilms et al., 2000). The abundance would be determined from the data after a successful spectral modeling but is left out of scope for this study.

We first estimate the values of the model parameters analytically. The ionization parameter of the corona at $r_{\mathrm{in}}$ is given by

in the unit of erg s^-1 cm. This is almost constant throughout the corona due to the assumed radial dependence of $n_{e}(r)$. The electron column density $N_{\mathrm{e}}$ in the line of sight is

The electron scattering optical depth is given by

in which $\sigma_{\mathrm{T}}=6.65\times 10^{-25}$ cm² is the Thomson cross-section. We make the hydrogen and electron column densities equal as $N_{\mathrm{H}}=N_{\mathrm{e}}$. Finally, the observed luminosity $L_{\mathrm{obs}}$ is approximated as

We now substitute $L_{\mathrm{obs}}=1.4\times 10^{37}$ erg s^-1 (Table 1), $r_{\mathrm{out}}=4.8\times 10^{10}$ cm from the partial eclipse in the X-ray light curve (Ebisawa et al., 2001), and $\xi=10^{2.5}$ erg s^-1 cm from the observed lines, which we will derive later using a numerical calculation in § 3.3.1. The unknown variable $L_{\mathrm{WD}}$ is parameterized by $\alpha=L_{\mathrm{obs}}/L_{\mathrm{WD}}$, in which $0\leq\alpha\leq 1$. As a function of $\alpha$, we solve the other unknown variables $n_{\mathrm{e,in}}$, $N_{\mathrm{H}}$, and $r_{\mathrm{in}}$ normalized by $r_{\mathrm{out}}$ as $f=r_{\mathrm{in}}/r_{\mathrm{out}}$ (Fig. 4). The set of equations solves only for a limited range of $\alpha=0.19-0.29$, and the derived values of the others are limited to a narrow range of $n_{\mathrm{e,in}}=10^{13.4}$–$10^{13.7}$ cm^-3 and $f=$0.38–0.61. Within this range, $\tau_{\mathrm{es}}=$0.21–0.34 and $N_{\mathrm{H}}=10^{23.5}$–$10^{23.7}$ cm^-2.

The WD effective area and gravity are related to the luminosity $L_{\mathrm{WD}}=0.5-0.7\times 10^{38}$ erg s^-1 as

and

respectively. Here, $R_{\mathrm{WD}}$ and $M_{\mathrm{WD}}$ are the radius and mass of the WD, respectively, and $\sigma_{\mathrm{SB}}$ and $G$ are the Stefan-Boltzmann and gravity constants. For $M_{\mathrm{WD}}=1.2~{}M_{\odot}$, $R_{\mathrm{WD}}=3.8\times 10^{8}$ cm from the mass-to-radius relation (Nauenberg, 1972), which corresponds to $\sim$0.2 $r_{\mathrm{in}}$. Then, $T_{\mathrm{eff}}\sim 900$ kK and $\log{g}\sim 9.0$.

Some of these derived values can be verified by other relations not used above. First, the intrinsic luminosity $L_{\mathrm{WD}}$ is consistent with the luminosity to support steady nuclear burning on the WD surface of a mass of $\approx 1M_{\odot}$ (Ebisawa et al., 2001). Second, if the observed energy shift of $\sim$900 km s^-1 represents the radial velocity of the expanding corona:

in which $\dot{M_{\mathrm{w}}}$ is the wind mass loss rate and $m_{\mathrm{p}}$ is the proton mass. By substituting $\dot{M_{\mathrm{w}}}\sim 10^{-7}M_{\odot}$ yr^-1 (Ablimit & Li, 2015), we obtain $n_{\mathrm{e}}r^{2}\sim 4\times 10^{33}$ cm^-1, which agrees with the value derived from equation (1) within an order.

We next verify that the analytically estimated parameters of the model (§ 3.2) are consistent with the numerical calculation of the photo-ionization plasma. We model the corona using the 1D radiative transfer code xstar v2.58e (Kallman et al., 2004). The code solves the radiative transfer in the scheme using the two-stream (radially inward and outward directions) solver and the escape probability approximation. It calculates the ionization/recombination, excitation/de-excitation, and heating/cooling balance at each zone of the spherically symmetric plasma and generates the charge state distributions and the level populations at each zone as well as synthesized X-ray spectra inward and outward of the plasma.

The heating and cooling sources are only the radiations. The heating source, or the incident emission, is given by the spectral energy distribution (SED) and the intensity parameterized with $\xi$. The plasma is described by its density profile $n_{e}(r)$ and total thickness ($N_{\mathrm{H}}$). The turbulent velocity is fixed to 400 km s^-1 (§ 2.3). We left H and He at the solar metacility and C, N, O, Ne, Mg, Si, S, Ar, Ca, and Fe at a half solar value. We ignored other elements. The electrons have a Maxwellian energy distribution of a temperature ($T_{\mathrm{e}}$) for the NLTE condition, which is calculated to reach a thermal balance between the heating and cooling.

For the purpose of a qualitative assessment, we use nearly a slab geometry of a uniform density and a blackbody spectrum for the incident emission using the best-fit values in the phenomenological model (§ 2.3). We use the code in an open geometry, in which the emission from the photoionized plasma can escape both the inward and outward directions, and the effects of radiative excitation are taken into account. The radiative excitation is known to be important for X-ray spectra from photoionized plasmas without the incident emission in the beam such as in Seyfert 2 galaxies (e.g., Kinkhabwala et al., 2002), which should also apply to the accretion disk coronae of X-ray binaries.

A major discrepancy in the xstar fitting (Fig. 7 a) is too weak continuum emission in the model. This makes sense considering that most of the continuum emission comes from the scattered emission by electrons in the corona, not the diffuse emission. In the 1D calculation of xstar, the electron scattering is implemented as a pure continuum absorption, and no scattered emission is produced.

In MONACO, the electron scattering is considered as scattering, including multiple scattering in the 3D space. We decompose the MONACO model spectrum into the scattered and diffuse components (Fig. 7 b). A single photon cannot be attributable to either one of them; a photon can experience both electron scattering and resonance scattering before escaping from the system. Here, the resonance scattering is not scattering but is absorption immediately followed by emission, thus should be included in the diffuse component. We divided all the photons based on their last interactions in the system. As expected, the scattered emission accounts for a large part of the continuum emission.

Another discrepancy in the xstar fitting (Fig. 7a) is the too strong line emission. The lack of scattered emission in the model accounts for this partially, but not fully. To illustrate this, we compare the xstar and MONACO spectra for the same parameters in Figure 8. For the transmitted spectra in panel (a), the continuum levels agree well between the two. For the diffuse spectra in panel (b), the radiative recombination continua at C VI, N VII, O VII, and O VIII agree in general. However, the line emission is systematically stronger for xstar than MONACO.

We argue two possibilities for the difference between the two solvers. One is that MONACO assumes that all excitation occurs only collisionally and does not include radiative excitation, which is known to enhance the emission lines (Chakraborty et al., 2021). The other is the assumption on the escape probability employed in xstar. The escape probability is an assumption necessary to solve the radiative transfer equations in reasonable computational resources in two-stream solvers, in which line photons escape from the system at a probability of $e^{-\tau(\nu)}$. Here, $\tau(\nu)$ is the opacity in the direction of the photon with a frequency of $\nu$. It has been claimed that solvers based on the escape probability assumption yield line emission stronger than the accelerated lambda iteration solver, which is considered more reliable for the lines (Hubeny, 2001). In fact, Dumont et al. (2003) argued that the escape probability method yields an overestimation of the line intensity, in particular for resonance lines like O VIII Ly$\alpha$ in optically thick cases, by an order.

The opacity at the line center is very large for the conspicuous lines. Taking the strongest line O VIII Ly$\alpha$ as an example, it is estimated by

in which $\Delta\nu_{\mathrm{D}}=\frac{E_{0}}{hc}\sqrt{\frac{2k_{\mathrm{B}}T}{m_{\mathrm{O}}}}$, $E_{0}$ is the energy of the line, $gf$ is the weighted oscillator strength, $A_{\mathrm{O}}$ is the oxygen abundance relative to hydrogen, $A_{8+}$ is the charge fraction of O⁸⁺, $A_{\mathrm{1s}}$ is the level population fraction of the ground state, $m_{\mathrm{O}}$ is the Oxygen atomic mass, $T$ is the plasma temperature, and other symbols follow their conventions. For $E_{0}=653$ eV, $f=0.14$, $N_{\mathrm{H}}=10^{23.5}$ cm^-2, $A_{\mathrm{O}}=0.5\times 4.9\times 10^{-4}$, $A_{8+}=0.5$ (Fig. 5), $A_{\mathrm{1s}}=1$, and $T=10^{6}$ K, we obtain $\tau_{\mathrm{Ly}\alpha}\sim 6.7\times 10^{4}\;\tau_{\mathrm{es}}\sim 1.4\times 10^{4}$.

Line photons scatter numerous times locally due to their large opacity at the line center until their wavelength diffuses into the dumping wing of the line profile. Then, the photons escape from the system in a single long flight with a significantly reduced opacity. This is a premise for the escape probability approximation implemented in xstar. In such a situation, the line profile is known to be double-peaked around the center. The analytical solution is obtained for a plane-parallel (Harrington, 1973; Neufeld, 1991) and spherical (Dijkstra et al., 2006) geometry of a uniform density. Despite the premise, xstar synthesizes emission lines with the Voigt profile as it is computationally challenging to calculate the diffusion over wavelengths.

Figure 9 shows a close-up view of the O VIII Ly$\alpha$ profile. The Voigt and Dijkstra profiles of no turbulence are shown with red and blue broken curves. They are smoothed with a Gaussian of the RGS energy resolution to compare to the data as shown with the solid curves of the same color. The smoothed Dijkstra profile matches better to the data than the smoothed Voigt profile. The turbulence velocity required in the phenomenological fitting (Fig. 2) can be explained by the line profile distortion by the radiative transfer effects. MONACO calculates the energy shift in every resonance scattering, thus the emergent line profile exhibits a distortion as shown with a cyan curve in Figure 2. This is too broad compared to the data, which we come back to later in § 5.5.

We discuss discrepancies found in several line ratios. First is the O VII He$\alpha$ triplet. The line ratios are different between xstar and MONACO most notably in their synthesized diffuse spectra (Fig. 8 b). Among the resonance ($r$), inter-combination ($i$), and forbidden ($f$) lines, $f$ is the weakest in both models, which is expected (§ 3.3.2) and agrees with the observation (Fig. 1). However, the $i$/$r$ ratio in the MONACO model is too strong compared to xstar and the observation. This stems from the limitation of MONACO, in which the level populations are calculated assuming that they are purely dominated by collisional processes. In the present application, both collisional and radiative processes matter, which are taken into account in the xstar calculation.

For the line ratios, we should focus on the xstar results. We examine the O VIII Lyman decrement (Ly$\beta$/Ly$\alpha$). The observed decrement is $\sim$1/3 (Fig. 2), while the calculated value is $\sim$0.05 (Fig. 6). The observed Ly$\beta$ is too strong with respect to Ly$\alpha$. The calculated decrement hardly changes if we change the model parameters. Another line ratio discrepancy is found in the four Fe XVII lines (Fig. 7 a), in which the observed 3C/3D line pair is too strong with respect to the 3F/3G pair. Both are primarily due to the attenuation of incident photons beyond the O VII edge.

A possible solution would thus be to introduce another spectral component such as the collisionally ionized emission independent of the O VII edge attenuation. A plausible origin of the collisionally ionized plasma is the shock of the expanding shell. If we assume that it is produced by the shock of the expanding velocity of $v\sim 900$ km s^-1 (Fig. 2), the temperature is given by $T=\frac{3}{16}\frac{\mu m_{\mathrm{p}}}{k_{\mathrm{B}}}v^{2}=1.0~{}\mathrm{keV},$ in which $\mu$ is the mean molecular weight, $m_{\mathrm{p}}$ is the proton mass, and $k_{\mathrm{B}}$ is the Boltzmann constant. Such a spectral component has been found in other SSS (Ness et al., 2022).

Based on the discussion above, we modify the spectral model for the final fitting. We start with the MONACO model (Fig. 7b) and add an apec model to account for the additional collisionally ionized plasma. We ignore the energy range of 0.565–0.569 and 0.715-0.740 keV to avoid the known caveats of the MONACO modeling for the O VII He$\alpha$ ($i$) and Fe XVII 3F/3G lines. The plasma temperature ($T$) and the normalization of the apec model were fitted, while the abundance was fixed to half the solar value. Both the MONACO and apec models were red-shifted collectively. The result is given in Figure 10.

The fitting is finally good enough (reduced $\chi^{2}=$1.5) to derive the best-fit parameter ranges as $\log{\xi}$ (erg s^-1 cm) $=$2.52–2.55, $\log{N_{\mathrm{H}}}$ (cm^-2) $=$23.49–23.50, $T_{\mathrm{eff}}$ (kK) $=$1000–1005, $v_{\mathrm{redshift}}$ (km s^-1) $=$ 895–941, and $T$ (keV) $=$ 0.24–0.25. The corona model parameters agree well with the analytical solution. The O VIII Lyman decrement is now better reproduced since the Ly$\beta$ is contributed by the apec component. The Fe XVII 3C line is also reproduced by the apec model, which was not accounted for in the corona model only with MONACO. The O VIII Ly$\alpha$ profile was better reproduced by a combination of the broad corona profile and the narrow collisionally ionized plasma profile. We still see some discrepancies such as Fe XVII 3D line at 0.81 keV and C VI radiative recombination continuum at 0.49 keV. These, along with the identified caveats and justification of the additional apec model, are left for future improvements of the codes.