Modeling Circumstellar Gas Emission around a White Dwarf Using cloudy

Gaia J0611−6931 has been observed extensively using multiple optical, infrared, and ultraviolet (UV) spectrographs (Dennihy et al. 2020; Melis et al. 2020; Owens et al. 2023; Rogers et al. 2024). In this paper, we focus on data from X-SHOOTER on the Very Large Telescope (VLT; program IDs: 104.C-0107, 105.202N, and 106.2130), the Cosmic Origins Spectrograph (COS) on the Hubble Space Telescope (HST; Rogers et al. 2024; ID: 16204), ¹⁵ and Flamingos-2 at the Gemini Observatory (GS-2021B-Q-244; Owens et al. 2023). The X-SHOOTER data were obtained on three different dates (i.e., 2019 October 15, 2020 December 3, and 2021 August 22 UT), and no significant changes were identified in the shape or the strength of the emission lines. Therefore, all the spectra were combined to have the highest signal-to-noise ratio (S/N). Stare mode was used for the first two epochs, while nod on slit was used for the last epoch, so the spectra from the near-IR (NIR) arm can be extracted. We used the standard ESO data reduction pipeline esoreflex (Freudling et al. 2013) and molecfit (Kausch et al. 2015; Smette et al. 2015). The complete spectrum of Gaia J0611−6931 is shown in Figure 1. The flux calibration is performed using photometry from SkyMapper (Onken et al. 2019). The wavelength is presented in vacuum throughout the paper.

Zoom In Zoom Out Reset image size

In the infrared beyond 1 μm, the X-SHOOTER data have higher spectral resolution than Flamingos-2, which is useful for checking the line identifications presented in Owens et al. (2023). However, the individual exposure time is 300 s for the NIR arm, which may be too long for proper sky subtraction. That might explain the apparent difference in the line strength between the X-SHOOTER and Flamingos-2 data.

As an additional check, we obtained a medium-resolution (R ≈ 25,000) optical spectrum of Gaia J0611−6931 with the Magellan Inamori Kyocera Echelle (MIKE) spectrograph(Bernstein et al. 2003) on the 6.5 m Magellan Clay Telescope on 2021 August 27 (UT). The observation and data reduction were performed using CarPy (Kelson et al. 2000; Kelson 2003). The MIKE data provide a wavelength coverage of 3500–5060 Å for the blue arm and 5000–9400 Å for the red arm. The total exposure time is 3111 s, and the S/N is 12 at 4000 Å and 14 at 6370 Å. We did not notice any difference between the X-SHOOTER and MIKE spectra.

The spectral classification of Gaia J0611−6931 is DAZ, which means hydrogen lines are the strongest features in the optical spectrum, as can be seen in Figure 1. There are also narrow absorption features from elements heavier than helium in the spectra (giving "Z" in the spectral classification), and the abundance analysis is presented in Rogers et al. (2024). In this paper, we focus on the emission features from the circumstellar gas around Gaia J0611−6931.

We performed an independent emission line analysis of Gaia J0611−6931 using all available spectra. The identifications for the optical spectra are straightforward because most emission lines are resolved and not blended with other features. There is one triple-peaked feature that extends between 5160 and 5190 Å, and it has been previously attributed to a blend of Fe ii 5170 Å and the Mg i triplet at 5169, 5174, and 5185 Å (Dennihy et al. 2020; Melis et al. 2020). As discussed in Section 4.2, we think this feature is dominated by the Mg i triplet with little contribution from Fe ii. The centroid of the Mg i 8809 Å line is shifted by about −90 km s⁻¹ compared to all the other emission lines. It is unclear what the cause of this is. There may be a similar wavelength shift for the Mg i 3832 Å line too, though the emission line is a lot weaker, making it difficult to confirm. In the near-UV data, there is a broad emission feature around 2800 Å that comes from four Mg i lines (see Figure 2). The infrared line identification in the Flamingos-2 data is more challenging because many lines are blended, and our result is consistent with those reported in Owens et al. (2023).

Zoom In Zoom Out Reset image size

To characterize the emission features, we measured the full width at zero intensity (FWZI) and the equivalent width (EW). For the EW, the uncertainty is dominated by the choice of continuum. For each line, we did several measurements selecting different wavelength regions for estimating the continuum flux. We report the average and the standard deviation of the measurements as the EW and its uncertainty, respectively. We also assume a 3% minimum uncertainty in the EW of the optical lines given the S/N of the spectra. The Si i 3907 Å, Si i 4104 Å, and Ca i 3935 Å lines sit on the slope of the broad hydrogen lines and therefore have large uncertainties in the EW. The Ca ii 3970 Å line is also detected, but it is difficult to perform any measurements due to its proximity to H 3971 Å. The EWs and FWZIs of the infrared lines are taken from Owens et al. (2023), who used the Flamingos-2 data. All the measurements are reported in Table 2.

We also measure the integrated line flux of all the emission lines, which is the flux above the continuum in absolute flux units. Note that in some cases, the photospheric absorption is visible in the profile, such as Ca ii 3935 Å and O i 8488 Å. However, their overall contribution is small and therefore not subtracted out in the measurement of line flux. In the case of Mg ii 2796 Å, the contribution from the photospheric line is significant, as shown in Figure 2. The integrated line flux of the Mg ii absorption lines from the model spectrum was added back to derive the flux of the emission feature around Mg ii 2976 Å. Similar to the EW measurements, the uncertainty of the line flux is dominated by the choice of continuum. For each line, we did several measurements selecting different wavelength ranges given the uncertainties in the FWZIs. We also assume a 3% minimum uncertainty in the line flux for the optical lines and a 10% minimum uncertainty for the infrared lines based on the S/N of the data.

As discussed in Gentile Fusillo et al. (2021), the EW of an emission line is a relative number, which is normalized by the continuum flux. In comparison, the line flux is an absolute number above the continuum flux, which is a more intrinsic way to characterize the emission line strength. For example, Table 2 shows that the EW of the Ca ii 3934 Å line is approximately 10 times weaker than that of the Ca ii 8544 Å line. However, that is not due to the relative strength of the emission line but rather to the difference in the white dwarf continuum flux. As shown in Figure 1, the white dwarf continuum flux is about a factor of 10 higher around Ca ii 3934 Å than the Ca infrared triplet region. Therefore, the line strengths for the Ca ii 3934 Å and Ca ii 8544 Å lines are actually similar. For the remainder of the manuscript, we will use line flux as the main observable to compare with the model outputs.

Figure 3 shows a compilation of all the unblended and resolved emission features. The overall shape of the lines is a bit different from element to element, but they all have a similar FWZI of 1530 km s⁻¹, half of which corresponds to a Keplerian velocity at a radius of 20R_wd assuming the gas is in a circular orbit and edge-on. If the disk is inclined relative to our line of sight, the projected velocity would correspond to a smaller radius. The half-peak separation is about 350 km s⁻¹, which corresponds to a Keplerian velocity at a radius of 70R_wd. The photospheric component can be seen as the narrow absorption feature around 0 km s⁻¹ in the Ca ii 3935 Å and O i doublet at 8448 Å.

Zoom In Zoom Out Reset image size

Here, we explore modeling the emission lines around Gaia J0611−6931 using the photoionzation code Cloudy (Ferland et al. 2017). Cloudy version 17.03 was used for this work. The setup is very similar to those described in Steele et al. (2021). Our configuration uses the cylinder geometry and assumes that the circumstellar gas has a circular orbit and constant density and is edge-on. Cloudy is a 1D code, and therefore it only computes the radiative transfer in the radial direction. The central white dwarf is the only ionizing source, and its parameters are listed in Table 1. The main free parameters for the Cloudy calculation are the radial extent of the disk, the hydrogen number density, and the abundances of different elements, the latter two affecting both the total mass of the gas and the excitation of the different elements.

For Cloudy outputs, we saved the predicted line intensities for a list of lines, which can be directly compared to the measured line fluxes in Table 2. The line emissivity file, which reports the emergent emissivity as a function of cloud depth, was also saved to calculate the line profile (see Section 3.2). Other output files were also saved following Steele et al. (2021).

We adopted an analytical method to compute the line profile, as illustrated in Figure 4. An additional assumption here is that the gas disk is optically thin; therefore, we can see the radiation all the way to the inner radius of the disk. As shown in Section 4.2, this assumption is not supported by the Cloudy calculations. Using polar coordinates, the disk can be approximated as many small patches along the radial and the azimuthal directions. The Keplerian velocity V_d at a given distance d from the white dwarf is

$$\begin{eqnarray}&&{V}_{{\rm{d}}}=\sqrt{\displaystyle \frac{{{GM}}_{\mathrm{wd}}}{d}}=753\,\mathrm{km}\,{{\rm{s}}}^{-1}\sqrt{\displaystyle \frac{20{R}_{\mathrm{wd}}}{d}},\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant and M_wd and R_wd are the mass and radius of the white dwarf. The velocity along the line of sight is V_d,α = V × sin α, where α is the azimuthal angle of the patch.

Zoom In Zoom Out Reset image size

For each patch, the emissivity is (d), which is a function of the distance from the white dwarf. The emissivity is also a direct output from Cloudy. The line profile of the emission line is calculated as the sum of the emissivity of the patch in both the azimuthal and the radial directions, Σ_d,α (d), which is a function of the line-of-sight velocity (or wavelength).

In this section, our goal is to use Cloudy to better constrain the radial extent of the disk, i.e., the inner disk radius R_in and the outer disk radius R_out. The gas disk is estimated to extend from 20R_wd to 70R_wd based on the profiles of the double-peaked emission lines (see Section 2.2).

The hydrogen number density is another basic input for Cloudy, but this value is difficult to constrain because hydrogen emission lines are not detected. For ease of calculation, we fix the ratio between hydrogen and calcium to be 1 (by number), because it is possible to estimate the density of calcium. Assuming a constant density profile, we can calculate the number density of calcium n(Ca) as

$$\begin{eqnarray}&&n(\mathrm{Ca})=\displaystyle \frac{{M}_{\mathrm{Ca},\mathrm{gas}}}{{u}_{\mathrm{Ca}}}\times \displaystyle \frac{1}{\pi {\left(100{R}_{\mathrm{wd}}\right)}^{2}\times h},\end{eqnarray} \tag{ 2 }$$

where M_Ca,gas is the mass of calcium in the circumstellar gas, u_Ca is the atomic mass of calcium, and $\pi {(100{R}_{\mathrm{wd}})}^{2}\times h$ is the estimated volume for the gas, assuming a constant disk scale height h of 1R_wd. M_Ca,gas can be estimated from the mass of Ca in the white dwarf's photosphere M_Ca,phot (7 × 10¹¹ g from Rogers et al. 2024), the settling time of Ca in the white dwarf's photosphere τ_Ca (10^−2.539 yr from Dufour et al. 2017), and the accretion timescale of the disk τ_diff,

$$\begin{eqnarray}&&{M}_{\mathrm{Ca},\mathrm{gas}}=\displaystyle \frac{{M}_{\mathrm{Ca},\mathrm{phot}}}{{\tau }_{\mathrm{Ca}}}\times {\tau }_{\mathrm{diff}}.\end{eqnarray} \tag{ 3 }$$

Using the α disk model, τ_diff can be estimated using Equation (14) of Goksu et al. (2023). Putting in all the numbers, the calcium number density is

$$\begin{eqnarray}&&n(\mathrm{Ca})=4.2\times {10}^{7}\,{\mathrm{cm}}^{-3}\times \displaystyle \frac{1{R}_{\mathrm{wd}}}{h}\times \displaystyle \frac{{10}^{-2}}{\alpha }.\end{eqnarray} \tag{ 4 }$$

The viscosity parameter α is uncertain by a few orders of magnitude. Therefore, we explore the calcium number density n(Ca) over a range of values.

We ran a few Cloudy models using different values of disk scale height and found it is proportional to the absolute strength of the emission lines. A larger scale height means a larger emitting area and therefore a stronger emission line flux. For this work, the disk scale height is kept at 1R_wd. Note that 1R_wd is just a convenient value to adapt; it is larger than the hydrostatic value, which depends on the disk temperature and radius (Goksu et al. 2023).

We compute a grid of Cloudy calculations with an inner radius R_in of 15R_wd, 20R_wd, and 25R_wd and an outer radius R_out of 40R_wd, 70R_wd, and 100R_wd. The result is shown in Figure 5. The separation of the double peak is very sensitive to the outer radius of the disk, but not as much to the inner disk radius. This is a result of the assumed constant density profile of the gas, where most of the disk mass is in the outer part. Therefore, a different density profile will change the radial extent of the disk. The model that best matches the line profile goes from 20R_wd to 100R_wd.

Zoom In Zoom Out Reset image size

Here, we attempt to constrain the chemical composition of the circumstellar gas. We focus on the isolated lines listed in Table 2. Mg i 11831 Å is excluded from this analysis because it is likely to be blended with other lines due to its large FWZI, even though a specific line has not been identified. For oxygen, there are no isolated lines, so both the O i triplet around 7774 Å and the O i doublet at 8488 Å are used. The disk is fixed at 20R_wd–100R_wd, the optimal parameters from Section 4.1.

For each element, we used Cloudy to compute the curve of growth, which shows the line flux as a function of the number density. In this iteration, we assumed the circumstellar gas has the same abundance as the bulk Earth and varied the total amount of material. The results are shown in Figure 6. The lines that are mostly likely to be optically thin are the O i 7774 triplet, Mg i 17113 Å, Si i 12274 Å, and Ca ii 8500.4 Å. However, given the observed line flux, most of the observed lines are optically thick.

Zoom In Zoom Out Reset image size

As an additional check, we consider the limiting case where the emission lines are completely saturated and the line flux reaches the Planck blackbody curve with a temperature equal to the gas temperature. In that case, the emission line ratios would be equal to the blackbody flux ratios. The result for this limiting case is presented in Figure 7. Within the uncertainties, the observed line ratios for O, Mg, and Ca all agree with the blackbody ratios for a temperature between 4500 and 10,000 K (the computed gas temperature from Cloudy; see Section 4.3). For these elements, only an upper limit to the number density can be derived. The only possible exception is Si. Si i 12274 Å can remain optically thin at a number density of 10^9.8 cm⁻³, the highest of all the detected emission lines. A silicon density may be derived using this line.

Zoom In Zoom Out Reset image size

One more iteration of the Cloudy model was computed to find the abundance that best matches the observed line strengths. However, because the number densities of elements are similar, the abundances of one element can affect the line strength of the other element. More tweaks of the abundances are needed to find the best model. Our best-fit abundance is listed in Table 3, which is able to reproduce the observed line strength within a factor of 2. All the lines have very large optical depths, and therefore only a lower limit can be derived for the density. The hydrogen number density is kept at 10^5.7 cm⁻³, the same as the number density of calcium. The radial extent of the disk is fixed at 20–100R_wd. As an additional confirmation, different Cloudy output files are checked to ensure that no other strong lines (such as hydrogen emission or forbidden lines) appear in the model.

As a test to the line fitting code presented in Section 3.2, we computed a line profile for Si i 12274 Å and compared it with data, as shown in Figure 8. The match is respectable, but the spectral resolution of the data is too low to make a proper comparison.

Zoom In Zoom Out Reset image size

There is a three-peaked feature around 5170 Å, which is a blend of three Mg i lines from 5168.8, 5174.1, and 5185.1 Å and possibly Fe ii 5170 Å (Dennihy et al. 2020). Our current model with only Mg already overpredicts the observed feature, and we conclude that the contribution from Fe ii 5170 Å in this feature should be minimal.

The gas temperature profile from the best-fit Cloudy model is shown in Figure 9. The temperature ranges from 10,000 K around the innermost radius at 20R_wd to 4500 K around the outermost radius at 100R_wd. The radiation from the white dwarf is the main heat source, and cooling is from various emission lines. The gas temperature is similar to the analytical solution proposed by Melis et al. (2010). They described a "Z ii" region model, which is akin to an H ii region, and the gas is heated by UV light from the white dwarf and cools by optically thick emission lines. As mentioned in Steele et al. (2021), Cloudy computes a self-consistent temperature profile, which is an advantage compared to other disk modeling papers assuming an isothermal disk (e.g., Fortin-Archambault et al. 2020; Budaj et al. 2022).

Zoom In Zoom Out Reset image size

Gaia J0611−6931 has a bright infrared excess with a fractional luminosity of 5% (Dennihy et al. 2020). In fact, a flat disk model failed to reproduce such a strong infrared excess, and the disk is either warped or has a significant amount of optically thin dust (Owens et al. 2023). Assuming there is dust around the same region as the circumstellar gas, we compute two dust temperatures using the optically thick model from Jura (2003) and an optically thin model (assuming a blackbody). In both cases, the dust temperature is much cooler than the gas. It is difficult for the dust and gas to occupy the same region given the large temperature difference. Alternatively, the dust and gas may not occupy the same region, as proposed for the eccentric disk around SDSS J1228+1040 (Goksu et al. 2023). Future work that incorporates both the dust and gas is needed to properly model these systems.

Rogers et al. (2024) reported detections of photospheric C, S, P, Al, Fe, and Ni in Gaia J0611−6931 that were not detected in the circumstellar gas. We use the Cloudy model to make some predictions of the strongest emission lines from these elements, as listed in Table 4. The lines are separated into UV, optical, and infrared wavelengths. As discussed in Section 2.2, generally speaking, the UV and optical line fluxes are larger compared to the infrared lines. However, the detectability of the emission lines is better in the infrared, because the white dwarf continuum flux is smaller. Effectively, it means it is possible to detect less luminous lines in the infrared. It will be useful to search for these additional emission lines to further constrain the composition of the circumstellar gas and compare with the photospheric abundances. A more extensive study exploring the detectability of the emission lines as a function of white dwarf and disk parameters will be presented in a future work (A. Steele 2024, in preparation).

We can now compare the abundances of the circumstellar gas around Gaia J0611−6931 with the photospheric abundances. The comparison is shown in Figure 10. The photospheric abundance is taken from Rogers et al. (2024) without any additional correction. The abundance difference between the buildup (the observed abundance in the white dwarf atmosphere equals that of the parent body) and steady state (the observed abundance is the abundance of the parent body modified by differential settling; see details in Dufour et al. 2017) is within 0.2 dex, which does not make a difference in the interpretation here.

Zoom In Zoom Out Reset image size

Only upper limits were obtained for the elements in the circumstellar gas. Our Cloudy model assumes that the hydrogen abundance is the same as calcium, because no hydrogen emission line is detected. In Figure 10, we arbitrarily subtracted 10 dex in the number ratios for circumstellar gas to put the compositions of the photosphere and the gas in a similar scale. For Gaia J0611−6931, we found that the upper limit to Si is the most constraining in the gas.

The mass of the gas disk can be estimated using the number densities in Table 3. The other assumption is that the disk extends from 20R_wd to 100R_wd with a constant scale height of 1R_wd. We can derive a lower limit on the mass of the gas disk mass to be 1.8 × 10¹⁹ g, considering O, Si, Ca, Na, and Mg. On the other hand, the total amount of heavy elements in the atmosphere of Gaia J0611−6931 is 7.1 × 10¹⁶, which is much smaller than the minimum mass of the gas disk. Taking the mass accretion rate of 6.57 × 10⁸ g s⁻¹ (Rogers et al. 2024), the minimum gas disk lifetime is 900 yr.

Our derived minimum gas disk mass of 1.8 × 10¹⁹ g is close to the gas mass estimate of 10²¹ g for another white dwarf, SDSS J1228+1040, using emission line kinematics (Goksu et al. 2023). In other studies to model the circumstellar gas absorption around white dwarfs, the estimated gas mass is about 10¹⁶ g for WD 1145+017 (Fortin-Archambault et al. 2020) and 5 × 10⁸–1.3 × 10¹⁶ g for WD 1124−293 (Steele et al. 2021). Compared to emission line detection, circumstellar absorption lines appear to be much more sensitive to small amounts of circumstellar material around white dwarfs.

In this paper, we developed a simple model for the circumstellar gas around the white dwarf Gaia J0611−6931 using Cloudy. The number densities of O, Si, Ca, Na, and Mg are constrained, and a minimum gas disk mass is also derived. Here, we discuss the limitations of the current work and future improvements.

Modeling: disk geometry, density profile, and scale height. Our Cloudy model takes the simplest assumption that the gas disk is edge-on with a constant density profile and a constant scale height. This is a major uncertainty that affects the projected velocities, radiative transfer, disk temperature, and line flux. In addition, radiative transfer is only computed along the radial direction and not in the vertical direction. Given the large optical depth in the radial direction, the line photons may escape from the disk surface. Future work is needed to explore these possibilities.

Observation: emission lines with resolved profiles. This study demonstrates the advantage of using infrared emission lines, which are more likely to be optically thin and therefore are better probes of the density of the gas. However, due to the lower spectral resolution of the infrared data, the line profiles are not resolved and most of the infrared lines are heavily blended, limiting its usage. Future observations with a higher spectral resolution that resolves the profiles of the infrared emission lines will be particularly useful in constraining the parameters of the disk and the abundances. It will also be useful to search for additional emissions from other elements (such as those listed in Table 4) that are currently not detected.