Spectral Modeling of the Supersoft X-Ray Source CAL87 Based on Radiative Transfer Codes

Supersoft X-ray sources (SSSs) are white dwarf (WD) binaries that radiate almost entirely below ∼1 keV. Their X-ray spectra are often complex when viewed with the X-ray grating spectrometers, where numerous emission and absorption features are intermingled and hard to separate. The absorption features are mostly from the WD atmosphere, for which radiative transfer models have been constructed. The emission features are from the corona surrounding the WD atmosphere, in which incident emission from the WD surface is reprocessed. Modeling the corona requires different solvers and assumptions for the radiative transfer, which has yet to be achieved. We chose CAL87, an SSS in the Large Magellanic Cloud, which exhibits emission-dominated spectra from the corona, as the WD atmosphere emission is assumed to be completely blocked by the accretion disk. We constructed a radiative transfer model for the corona using two radiative transfer codes: xstar for a one-dimensional two-stream solver and MONACO for a three-dimensional Monte Carlo solver. We identified their differences and limitations in comparison to the spectra taken with the Reflection Grating Spectrometer on board the XMM-Newton satellite. We finally obtained a sufficiently good spectral model of CAL87 based on the radiative transfer of the corona plus an additional collisionally ionized plasma. In the coming X-ray microcalorimeter era, it will be required to interpret spectra based on radiative transfer in a wider range of sources than what is presented here.


INTRODUCTION
Super Soft X-ray Sources (SSS) represent a particular class of the white dwarf (WD) binary systems characterized by their extremely soft X-ray emission below ∼1 keV.The X-ray luminosity reaches up to ∼10 38 erg s −1 (Greiner 1996) close to the Eddington limit for a solar-mass WD.They are observed either as a persistent SSS with steady nuclear burning on the WD Corresponding author: Masahiro Tsujimoto tsujimot@astro.isas.jaxa.jp2013).In the corona, incident photons from the WD atmosphere are reprocessed.Ness et al. (2013) conducted a systematic study of high-resolution X-ray spectra of SSS and found that they are categorized into two types; those mainly dominated by absorption features (SSa) and others by emission features (SSe).Ness et al. (2013) proposed an interpretation that the X-ray spectra of SSa and SSe respectively have a lower and higher fraction of the corona emission with respect to the WD surface emission.The difference is caused by different degrees of the obscuration of the WD surface emission by the companion star or the accretion disk.Ness et al. (2013) further pointed out that SSS with high inclination angles (≳ 70 • ) tends to have the SSe spectra.This interpretation is reinforced in the partially eclipsing system V5116 Sgr, in which the SSa spectrum is seen out of the eclipse and the SSe spectrum is seen during the eclipse (Sala et al. 2008(Sala et al. , 2010)).
The next obvious step is to interpret the observed spectra with physical models and constrain their physical quantities.This has been explored intensively for SSa and was applied to explain some aspects of the observed data.Lanz et al. (2005) extended an industrystandard radiative transfer code for the stellar atmosphere TLUSTY into the WD regime and explained the CAL 83 spectra.Rauch et al. (2010) developed another model called TMAP (Tübingen Model Atmosphere Package) and fitted the SSS phase spectra of the nova V4743 Sgr.Both models solve the radiative transfer for the non-local thermal equilibrium (NLTE) condition under the plane-parallel geometry and hydrostatic density distribution using the accelerated lambda iteration solver.van Rossum & Ness (2010) further explained the V2491 Cyg spectra using the PHOENIX code calculating the WD atmosphere with an expanding envelope.
The SSe spectra, on the contrary, are less explored because of the increased complexity.In addition to the WD atmosphere model, we need to construct the corona model.The corona is photo-ionized by the WD atmosphere emission (hereafter called "incident emission").Reprocessed photons are emitted through radiative deexcitation and recombination ("diffuse emission") and electron scattering of the incident photons ("scattered emission").NLTE radiative transfer calculation is undoubtedly needed, but different solvers and assumptions may be required from the ones used for the WD atmosphere models.
The purpose of this paper is to construct a physical model for SSe.We choose CAL87, a representative SSe source, which is simple enough so that we can constrain the system parameters based only on analytical calcula-tions.We construct corona models using two different codes using different solvers and compare their results with observations and discuss their differences and limitations.
The structure of this paper is as follows.We describe the data in § 2 and the model in § 3.In the data section ( § 2), we give properties of the target ( § 2.1), the observation and data reduction ( § 2.2), and the characterization of the spectra using phenomenological models ( § 2.3).In the model section ( § 3), we set up the assumption for the geometry and the parameterization ( § 3.1), based on which we give an analytical estimation of the parameter ranges ( § 3.2) and numerical verification of the parameter values by comparing the data and the synthesized spectra ( § 3.3).In § 4, we construct physically motivated models using two different codes: one is xstar using a one-dimensional (1D) two-steam solver ( § 4.1) and the other is MONACO using a three-dimensional (3D) Monte Carlo solver ( § 4.2) of the radiative transfer.We identify discrepancies between the two models against the observation and discuss their possible origins in § 5.The quoted errors are 1σ statistical uncertainty throughout the paper.
2. DATA 2.1.Target CAL87 (Long et al. 1981) is a reasonably bright source located in the Large Magellanic Cloud (LMC) at a distance of 48.1 ± 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 ∼5 ×10 10 cm, which is ∼0.4 times the orbital separation, and the inclination angle is ∼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.

Observations and Data reduction
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 ′′ and an annulus of the inner and outer radius of 75 ′′ and 100 ′′ .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).

Spectral characterization
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 Gal H ) was fixed at 7.58×10 20 cm −2 (Dickey & Lockman 1990) and that of the latter (N LMC H ) 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.9keV 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 obs ) is 1.4 × 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α and Lyβ, O VII Heα, N VII Lyα and Lyβ, and four Fe XVII lines for the transition to the ground state (2s) 2 (2p) 6 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α and Lyβ 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α 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.74keV, we modeled with two lines (3F and 3G; Table 2).For the Fe XVII line complex at 0.81-0.83keV, we modeled with two lines (3C and 3D; Table 2).The The RGS spectrum in the 14-31 Å band (corresponding to ∼0.4-0.9 keV), where the vertical axis is in the linear scale.The identified emission lines are marked with vertical lines and labeled in black at their rest wavelengths.Those absent lines which would have been expected at slightly different ionization states (see Fig. 5) are also indicated in grey.The cyan-shaded region and the red solid line indicate the data and the best-fit model, respectively.
b Weighted oscillator strength.
c Emissivity for the collisionally ionized plasma at a temperature of 0.6 keV (Smith et al. 2001).
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 ∼900 km s −1 and broadening of ∼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β deviates from the trend, which may be biased by the RGS chip gap eliminating a large fraction of the line profile.

Assumptions
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 eff , a surface gravity of g, and a luminosity of L 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 in and r out , respectively.The electron density profile is assumed to be n e (r) = n e,in (r/r 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.

Analytical estimation
We first estimate the values of the model parameters analytically.The ionization parameter of the corona at r 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 e in the line of sight is The electron scattering optical depth is given by in which σ T = 6.65 × 10 −25 cm 2 is the Thomson crosssection.We make the hydrogen and electron column densities equal as N H = N e .Finally, the observed luminosity L obs is approximated as We now substitute L obs = 1.4 × 10 37 erg s −1 (Table 1), r out = 4.8 × 10 10 cm from the partial eclipse in the X-ray light curve (Ebisawa et al. 2001), and ξ = 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 WD is parameterized by α = L obs /L WD , in which 0 ≤ α ≤ 1.As a function of α, we solve the other unknown variables n e,in , N H , and r in normalized by r out as f = r in /r out (Fig. 4).The set of equations solves only for a limited range of α = 0.19 − 0.29, and the derived values of the others are limited to a narrow range of n e,in = 10 13.4 -10 13.7 cm −3 and f =0.38-0.61.Within this range, τ es =0.21-0.34 and N H = 10 23.5 -10 23.7 cm −2 .
The WD effective area and gravity are related to the luminosity L WD = 0.5 − 0.7 × 10 38 erg s −1 as and respectively.Here, R WD and M WD are the radius and mass of the WD, respectively, and σ SB and G are the Stefan-Boltzmann and gravity constants.For M WD = 1.2 M ⊙ , R WD = 3.8 × 10 8 cm from the mass-to-radius relation (Nauenberg 1972), which corresponds to ∼0.2 r in .Then, T eff ∼ 900 kK and log g ∼ 9.0.Some of these derived values can be verified by other relations not used above.First, the intrinsic luminosity L WD is consistent with the luminosity to support steady nuclear burning on the WD surface of a mass of ≈ 1M ⊙ (Ebisawa et al. 2001).Second, if the observed energy shift of ∼900 km s −1 represents the radial velocity of the expanding corona: in which Ṁw is the wind mass loss rate and m p is the proton mass.By substituting Ṁw ∼ 10 −7 M ⊙ yr −1 (Ablimit & Li 2015), we obtain n e r 2 ∼ 4 × 10 33 cm −1 , which agrees with the value derived from equation ( 1) within an order.

Numerical calculation
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 twostream (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 ξ.The plasma is described by its density profile n e (r) and total thickness (N 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 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 bestfit 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.

Presence and absence of lines
We first confirm the presence or absence of several lines expected at particular ionization states (Fig. 1b), which gives a constraint on the charge state distribution and the ionization parameter ξ.Using parameters in the conceivable range as N H = 10 23.5 cm −2 and n e = 10 14 cm −3 , we calculated the fraction of ionization states of C, N, O, and Fe averaged over the corona volume (Fig. 5).In the observed spectrum (Fig. 1b), we identified N VII lines (Lyα and Lyβ) but not N VI lines (Heα) produced by the recombination cascade respectively of N 7+ and N 6+ ions; these conditions give the lower ionization limit log ξ ≳ 2. produced in the same process of O 8+ and O 7+ ions; these conditions in contrast gives the upper-lmit log ξ ≲ 3.5.Over this ξ range, most Fe stays in the Ne-like (Fe 16+ ), which is consistent with the presence of Fe XVII features produced by radiative decay from excited levels to the ground state of Fe 16+ ions.The absence of Fe XVI and Fe XVIII features is due to the small fraction of Fe 15+ and Fe 17+ ions, respectively.These conditions are consistent with the estimated ξ range of 10 2.2 − 10 3.5 .

Line ratios
We next examine the intensity ratio of the identified lines.A well-established method is to use the Heα triplet line intensity ratios of O VII, namely the resonance (r), inter-combination (i), and forbidden lines (f ).G ≡ (f + i)/r and R ≡ f /i ratios are often used as temperature and density indicators when collisional processes dominate (Gabriel & Jordan 1969) or the ion-ization and column density indicators when radiative processes dominate (Porter & Ferland 2007).We refrain from using these diagnostics to constrain model parameters, though, given the poor constraints by the data (Fig. 1a inset).We only argue that it is reasonable to have the r line to be the strongest and the f line to be the weakest.The reasons for the former are that (1) T e = 0.2 − 2.0 MK from the thermal balance is close to the maximum formation temperature of the line in the collisional process and (2) the line is enhanced by the continuum pumping in the radiative process (Chakraborty et al. 2021).The reason for the latter is that n e ∼ 10 14 cm −3 ( § 3.2) is much higher than the critical density of 3.4×10 10 cm −3 for O (Blumenthal et al. 1972).

PHYSICAL MODELING
Based on the comparison of the observed spectra ( § 2.3) against the analytical ( § 3.2) and numerical ( § 3.3) calculations, we have verified that our assumed model and the parameterization are appropriate.We have also constrained the range of parameters.In order to contain the parameters more accurately, we now construct physical models using two complementary methods: one is 1D two-steam solver (xstar; § 4.1) and the other is 3D Monte Carlo solver (MONACO; § 4.2) of the radiative transfer in the corona.

Setup
For the corona model, we use xstar in the same setup with § 3.3 except for two improvements; inclusion of the density profile n e (r) = n e,in (r/r in ) −2 and use of a WD atmosphere model spectrum for the incident emission.For the atmosphere model, we used TMAP (Rauch et al. 2010), which provides spectra1 at 10-8000 Å for T eff =405-1050 kK and log g (cm s −2 )=9 with varying abundances for major elements.(Fig. 5) and, resultantly, the emission line ratio of different charge states; e.g., O VIII Lyα and O VII Heα.
For log N H , the emission line ratios are not influenced, but the spectral shape above the edge energies change due to different columns in the line of sight.For log n H , the O VII triplet changes most significantly ( § 3.3.2).

Result
We ran the xstar calculations for a grid of three parameters: T eff , ξ, and N H .We skip n H as the constraining power of the O VII triplet in the data is poor.The Galactic and Magellanic interstellar absorption columns were fixed to the value derived in the phenomenological fitting in Table 1.The turbulent velocity is fixed to 400 km s −1 ( § 2.3).The redshift (v shift ) was fitted collectively.All the other parameters, including n H , are derived from the free parameter values based on the equations in § 3.2.
Figure 7 (a) shows the best-fit model and the residuals to the fit.The fitting is very poor.The best-fit parameters are away from those in the analytical estimates ( § 3).We identify several major discrepancies: emission lines are too strong, continuum emission is too weak, and the line profiles and ratios are poorly reproduced.Many of these are expected due to the known limitations of this solver, which we will discuss in § 5.

Setup
We next use the other radiative transfer code MONACO based on the Monte Carlo solver.MONACO is a photon tracking simulator based on the Monte Carlo framework (Watanabe et al. 2006a;Odaka et al. 2011a;Odaka 2012) and is applied to many astrophysical applications (e.g., Hagino et al. 2015;Tomaru et al. 2018;Tanimoto et al. 2019;Mizumoto et al. 2021).We used MONACO version 1.7 release candidate and the database version 1.7.
A number of photons are emitted from the source and individual photons are tracked for their interactions with the matter placed in a 3D space.Each photon is tracked until it either escapes from the system or is destructed inside the system by absorption.Interactions between the photons and matter are calculated, including radiative excitation/de-excitation, ionization/recombination, photoelectric absorption, and scattering by free and bound electrons.However, because the framework is based on photon tracking, emission processes initiated by electrons such as free-free emission, are currently not included.We collected photons that experienced at least one interaction in the corona for the secondary emission, which includes both the scattered and diffuse emission.
The design of the framework allows one to calculate electron scattering in a 3D space in a consistent manner.
The code does not calculate the thermal, ionization, and excitation balances unlike xstar.The electron density and temperature (n e and T e ) as well as the charge state distributions of ions at each location need to be given as inputs.We used the result of the xstar simulation with the best-fit parameters ( § 4.1).
In the current version, MONACO calculates the level populations determined purely by the collisional processes based on the given n e and T e .This is not correct for the photo-ionized plasmas, in which the radiative processes also contribute.Unlike SKIRT version 9.0, a widely used Monte Carlo solver (Vander Meulen et al. 2023), MONACO calculates the photon-ion interactions for the H-like and He-like ions of major metals.For this work, we made two modifications to the code and the database: (1) inclusion of lower ionized Fe ions than Fe 24+ down to Fe 16+ and (2) recording both the positive and the negative deposit energy for the resonance scattering.

Result
We used the same approach with the 1D two-stream solver ( § 4.1) for the fitting.The geometry (Fig. 3) is set up in a 3D space.The simulation was run at the same grid of the same free parameters (T eff , ξ, and N H ) with the turbulent velocity fixed, the bulk velocity collectively fitted, other corona parameters derived from the grid parameters based on the analytical relations ( § 3.2), and the Magellanic and Galactic absorption fixed to the value from the phenomenological fitting.
For each grid parameter, we launched half a million photons from a point source at the center.The photons have a direction selected randomly over a uniform distribution and an energy selected randomly over a uniform distribution in 0.35-1.24keV.The output was convolved with the input SED.For the grid of log ξ = 2.5, log N H = 23.5, and T eff = 10 3 kK as an example, (a) 4% photons are destructed in the corona by photo-electric absorption, (b) 72% photons escape from the system with no interactions, and (c) the rest escape with some interactions in the corona.We use (c) for the scattered plus diffuse emission.
Figure 7 (b) shows the best-fit model and the residuals to the fit.The fitting has improved with the best-fit parameters of log ξ = 2.6, log N H = 23.2, and T eff = 1015 kK.They are in reasonable agreement with the analytical solution, yet the fitting is still too poor to derive their uncertainties.Some discrepancies remain, in particular, in the line profile of O VIII Lyα, the line ratio within O VII Heα, and Fe XVII lines.We will discuss these issues in § 5. We constructed the spectral models of the corona based on two different solvers of the radiative transfer calculation and compared them to the observed spectra ( § 4).The xstar model yielded a poor fitting with the best-fit parameters far from those by the analytical estimates ( § 3.2).The MONACO model yielded a better fitting with the best-fit parameters consistent with the analytical estimates ( § 3.2).Both of them exhibited some discrepancies against the observations.We discuss possible reasons for four major discrepancies ( § 5.1- § 5.4) and present the final fitting by considering the model limitations ( § 5.5).

Continuum emission
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.

Line intensity
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 −τ (ν) .Here, τ (ν) is the opacity in the direction of the photon with a frequency of ν.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α 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α as an example, it is estimated by 3) is shown with the black curve.The phenomenological model has two unresolved Gaussian lines representing Lyα 1 and Lyα 2 with a fixed energy and intensity ratios, which are collectively shifted redward by 1×10 3 km s −1 and smoothed by 360 km s −1 (Fig. 2) and convolved with the detector line spread function (LSF).Upon this, compared are the Voigt (blue) and Harrington (red) profiles.The profiles are broadened only for the natural and the thermal broadening of 10 6 K and not by turbulence.The compared profiles are convolved with the Gaussian core of the LSF but no underlying continuum emission, thus the tail parts should be ignored.The cyan curve shows the profile calculated by MONACO with finite but negligible turbulence (50 km s −1 ) to stabilize the calculation.
N H = 10 23.5 cm −2 , A O = 0.5 × 4.9 × 10 −4 , A 8+ = 0.5 (Fig. 5), A 1s = 1, and T = 10 6 K, we obtain τ Lyα ∼ 6.7 × 10 4 τ es ∼ 1.4 × 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α 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.

Line ratios
We discuss discrepancies found in several line ratios.First is the O VII Heα 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β/Lyα).The observed decrement is ∼1/3 (Fig. 2), while the calculated value is ∼0.05 (Fig. 6).The observed Lyβ is too strong with respect to Lyα.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 ∼ 900 km s −1 (Fig. 2), the temperature is given by T = 3 16 µmp kB v 2 = 1.0 keV, in which µ is the mean molecular weight, m p is the proton mass, and k B is the Boltzmann constant.Such a spectral component has been found in other SSS (Ness et al. 2022).

Final fitting
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.569and 0.715-0.740keV to avoid the known caveats of the MONACO modeling for the O VII Heα (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 χ 2 =1.5) to derive the best-fit parameter ranges as log ξ (erg s −1 cm) =2.52-2.55,log N H (cm −2 ) =23.49-23.50,T eff (kK) =1000-1005, v 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β 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α 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.

SUMMARY AND CONCLUSION
We presented result of the spectral modeling of CAL87, a super-soft source dominated by emission features originating from the extended accretion disk corona around the WD.First, we presented a phenomenological fitting of the observed spectra with the EPIC-MOS and RGS instruments onboard XMM-Newton.Based on the fitting results and the system geometry, an analytical solution was obtained.
We then performed radiative transfer calculations using two codes employing different solvers and assumptions: xstar for a 1D, two-streamer solver and MONACO for a 3D, Monte Carlo solver.We fitted the spectral model to the data and identified discrepancies between the codes against the observation.
We discussed possible causes of these discrepancies in terms of the continuum emission, line intensity, line profile, and line ratio.We found that xstar excels in the level population calculations, while MONACO excels in the scattering (both electron and resonant line) calculations.We also argued for the presence of the additional collisionally ionized plasma.Finally, we presented the detailed spectral model that yields consistent parameters with the analytical model and agrees reasonably well with the observation.
In the coming microcalorimetry era in X-ray astronomy, the interpretation based on radiative transfer will be more seriously required than before.This is particularly the case for sources with hard X-ray spectrum dominated by reprocessed emission in the surrounding medium such as X-ray binaries during eclipse, accretion disk corona sources, and Seyfert 2 galaxies.Comparisons between different solvers and implementations, and comparisons against observations should be made in a wider range of applications than the one presented here.
We are grateful for our reviewer Jan-Uwe Ness for his insightful comments that not just improved the manuscript but gave us inspirations for the present work.We appreciate expert advice from Timothy Kallman at NASA GSFC for the implementation in xstar and Atsushi Tanimoto at Kagoshima University in MONACO.This work made use of the JAXA's supercomputing system JSS3.This work was supported by JSPS KAKENHI Grant Numbers JP14J11810, JP19H01906, JP19H05185, JP18H05861, JP21K13958, and JP24740190.

Figure 1 .
Figure 1.Fitting results with a phenomenological model.(a) MOS1 (black), MOS2 (green), and combined RGS (red) spectra in the 0.3-2.0keV band.Two lines indicate the photoelectric absorption K edges by O VII and O VIII.The upper panel shows the binned data with crosses and the best-fit models with solid lines.The inset shows a closeup view of the local fitting for O VII Heα lines.The middle panel shows the residuals to the fit, and the bottom panel shows the effective areas of the individual detectors.One of the two RGS units (RGS2) does not cover the 20.0-24.1 Å range, where N VI Lyβ and O VII Heα lines are located.(b)The RGS spectrum in the 14-31 Å band (corresponding to ∼0.4-0.9 keV), where the vertical axis is in the linear scale.The identified emission lines are marked with vertical lines and labeled in black at their rest wavelengths.Those absent lines which would have been expected at slightly different ionization states (see Fig.5) are also indicated in grey.The cyan-shaded region and the red solid line indicate the data and the best-fit model, respectively.

Figure 2 .
Figure 2. Results of the fitting of individual line complex for (a) energy shift, (b) line broadening (Gaussian standard deviation), and (c) intensity.The weighted mean values are shown with dashed lines in (a) and (b).The 1 σ statistical error bar is given.

Figure 3 .
Figure 3. Setup of the geometry.The notations are given in the text.

Figure 4 .
Figure 4. Analytical solution of (a) f = rin/rout, (b) column density through the corona NH, and (c) electron density at rin (solid) and volume mean (dashed) as a function of α = L obs /LWD.

Figure 5 .
Figure 5. Charge state distribution as a function of log ξ (erg s −1 cm) calculated by xstar for C, N, O, and Fe averaged over the corona volume.

Figure 6 Figure 6 .
Figure 6.Synthesized spectra of (a) incident emission from the WD atmosphere and (b-d) diffuse emission from the corona with different physical parameters.Note that the incident emission for the corona model is a blackbody emission shown with the dotted grey curve in (a-d).Spectra in (a) are normalized at 0.4 keV.For each of the model parameters (T eff , log ξ, log NH, and log nH), three representative values are chosen; one is a value close to the estimates ( § 3) in black (thus the black line spectra in the panels (b-d) are identical) and a lower (magenta) or a higher (cyan) value than that.Absorption edges are labeled in (a) and shown with dashed vertical lines.Lines are labeled in (d) and shown with dotted vertical lines.

Figure 7 .
Figure 7. Fitting result with the physical model based on (a) xstar and (b) MONACO.Symbols follow Figure 1.In (b), the model spectrum is decomposed into the scattered emission (yellow-green) and the diffuse emission (orange) based on the last interaction of the photons escaping from the system.

Figure 8 .
Figure 8.Comparison between the xstar and MONACO model spectra for the same setup (80 eV blackbody emission for the input with ξ = 10 2.5 erg s −1 cm and NH = 10 22.5 cm −2 ).(a) Incident (dashed) and transmitted (solid) emission, (b) diffuse emission, and (c) diffuse and scattered emission combined, which is available only for MONACO.

Figure 9 .
Figure 9. Profile of the O VIII Lyα line.The data are shown with black error bars, while the best-fit phenomenological model ( § 2.3) is shown with the black curve.The phenomenological model has two unresolved Gaussian lines representing Lyα 1 and Lyα 2 with a fixed energy and intensity ratios, which are collectively shifted redward by 1×10 3 km s −1 and smoothed by 360 km s −1 (Fig.2) and convolved with the detector line spread function (LSF).Upon this, compared are the Voigt (blue) andHarrington (red)  profiles.The profiles are broadened only for the natural and the thermal broadening of 10 6 K and not by turbulence.The compared profiles are convolved with the Gaussian core of the LSF but no underlying continuum emission, thus the tail parts should be ignored.The cyan curve shows the profile calculated by MONACO with finite but negligible turbulence (50 km s −1 ) to stabilize the calculation.

Figure 10 .
Figure 10.Final fitting using the corona model of MONACO (cyan) plus a collisionally ionized plasma model of apec (lime).The energy range of 0.565-0.569and 0.715-0.740keV are ignored to avoid the known caveats of the MONACO modeling for the O VII Heα (i) and Fe XVII 3F/3G lines.

Table 1 .
Best-fit parameters of the phenomenological model (continuum).

Table 2 .
Fe XVII lines