© The Authors 2024

1 Introduction

Despite the strong gravitational acceleration at their surface, 25%–50% of galactic white dwarfs (WDs) show evidence of metal pollution (Zuckerman et al. 2003; Koester et al. 2014). The effective temperature of these metal-polluted, so-called DZ type, WDs ranges from 5000 to 23 000 K, implying a cooling age of 40 million to 2 billion years (Hollands et al. 2017; Koester et al. 2014). This is at least an order of magnitude higher than the settling times of the pollutants, which ranges from days to a few million years (Paquette et al. 1986; Metzger et al. 2012). Metal pollution has been observed to persist for at least twelve years (Gänsicke et al. 2008; Hartmann et al. 2014, 2016; Manser et al. 2016b), and the fading of the phenomenon has not been reported.

Most studies assume that interplanetary matter, such as planets, moons, asteroids, comets, dust, and gas, is the source of metal pollution (Zuckerman et al. 2003; Dufour et al. 2007; Zuckerman et al. 2010; Kleinman et al. 2013; Koester et al. 2014; Gentile Fusillo et al. 2015; Kepler et al. 2015, 2016). The amount of accreted material can range from the mass of Phobos to that of Pluto (refer to Fig. 6 of Veras 2016).

The idea of interplanetary matter polluting WDs is supported by the fact that rocky minor planets have been discovered around DZ WDs. The first of such bodies was found in the WD 1856+534 system along with five additional small bodies (Croll et al. 2017; Gänsicke et al. 2016; Rappaport et al. 2016; Vanderburg et al. 2015; Xu et al. 2016). All of these minor planets orbit the WD on circular and coplanar orbits with periods of 4.5–4.9 h (corresponding to semi-major axes of 5 × 10⁻³ to 5.6 × 10⁻³ au) and have radii of 1–100 km.

Minor planets on such tight orbits have certainly crossed the Roche radius of the WD and have thus begun to shed mass. Other newly discovered rocky planets around DZ WDs also have in common an orbit that lies at or below the Roche limit (Manser et al. 2019; Vanderbosch et al. 2020; Guidry et al. 2021). Furthermore, the observation of asymmetric transit light curves strongly suggests a cometary tail of volatiles around these objects. Both observations imply that terrestrial planets are being seen during the process disruption around DZ WDs. We note that gas giants of with masses 2.5–14 times that of Jupiter have also been discovered around WDs; however, not all of their stars are metal polluted (Luhman et al. 2011; Gaensicke et al. 2019; Vanderburg et al. 2020; Mullally et al. 2024).

Over 50 DZ WDs have been observed to host dust discs (Bergfors et al. 2014), 21 of which also host a gaseous counterpart (Gänsicke et al. 2006, 2007, 2008; Gänsicke 2011; Farihi et al. 2012; Melis et al. 2012; Wilson et al. 2014; Dennihy et al. 2020; Gentile Fusillo et al. 2021; Melis et al. 2020). Currently, all DZ WDs that are known to host a gas disc also share a dusty component. However, many WDs have been observed that are only accompanied by a dust disc. We emphasise that both dust and gas discs are present in the systems that host rocky minor planets. The stability of discs during the giant phase of the host star is not possible due to strong stellar winds; thus discs must have formed after the star had become a WD (Veras 2016).

It has been observed that both the dust and the gas discs extend below the Roche radius of the WD. Dust discs are typically observed to lie between 0.6 R_⊙ and 1.2 R_⊙ (Farihi 2016), with a mass ranging from 10¹⁹ g to 10²⁴ g (Manser et al. 2016b). Gas discs, on the other hand, have been found found at 0.5 R_⊙−1.2 R_⊙ (Gänsicke et al. 2008; Farihi 2016) and have a mass ranging from 10¹⁴ g to 10¹⁹ g (Hartmann et al. 2016). The inner boundary of the gas discs should coincide with a temperature of 5000 K since no disc emission lines have been observed in the far-ultraviolet region of the spectra (Hartmann et al. 2016).

Gas discs are typically discovered through the infrared emission of the Ca II triplet. If the discs were Keplerian, the emission lines would show symmetrical double peaks (Horne & Marsh 1986). However, this is not the case for the gas discs of DZ WDs, which show an asymmetry in the intensity of the redshifted and blueshifted peaks. This is likely due to a non-Keplerian density or velocity distribution in the disc. The asymmetry has been observed to persist for as long as twelve years (Wilson et al. 2014; Manser et al. 2016b,a; Cauley et al. 2018; Dennihy et al. 2018, 2020), but its disappearance has also been observed (Gänsicke et al. 2008; Hartmann et al. 2014).

Additionally, a periodic fluctuation in the asymmetry of the Ca II lines has been detected. This is hypothesised to be caused by the precession of an eccentric disc (e.g. Cauley et al. 2018; Gänsicke et al. 2008; Hartmann et al. 2014; Manser et al. 2016b). The period of the change in asymmetry ranges from 123 minutes to 37 years (Hartmann et al. 2016; Manser et al. 2019). However, to accurately map these changes, it is necessary to measure at a cadence that is a fraction of the precession period (great examples are Dennihy et al. 2018; Manser et al. 2019).

The fact that the WDs at the heart of these systems are metal polluted means that the material of the discs is lost to the WD through accretion. Taking into account the timescale on which anomalous metallicity has been observed and the settling times of the pollutants, one must assume a continuous replenishment of the disc material (Koester 2009). This rules out the direct infall of a planet, asteroid, moon, or comet onto the star (Debes et al. 2012a; Veras et al. 2014, 2016; Payne et al. 2016; Bhatt 1985) as a pollution mechanism. Furthermore, the above scenario is unable to explain the origin of the discs around DZ WDs.

On the other hand, continuous replenishment requires the presence of a constantly disturbed asteroid belt (Bonsor et al. 2011; Debes et al. 2012a; Frewen & Hansen 2014; Veras et al. 2016; Nixon et al. 2020) or a disintegrating planet, asteroid, or comet (Jura 2003, 2008; Debes et al. 2012b; Veras et al. 2013). The latter scenario is considered to be the most likely origin of the dust discs (Debes & Sigurdsson 2002; Jura 2003; Malamud & Perets 2020a,b; Trevascus et al. 2021; Veras et al. 2014).

Numerical simulations suggest that bodies with semi-major axes of a few au can migrate to the Roche limit due to tidal dissipation or perturbations from other planets in the system and that terrestrial planets are more likely to migrate inwards than gas giants (Veras et al. 2016; O’Connor & Lai 2020). This, together with the high density and small size of the planets observed around DZ WDs and the spectra of the contaminated stars, strongly supports the hypothesis that Earth-like bodies are being disrupted around DZ WDs. The observed discs are devoid of hydrogen (Hartmann et al. 2011), which is also consistent with the above scenario.

With regard to the origin of the gas discs, some studies assume that the disintegrating body contains intrinsic gas so the disc is produced immediately as the body passes the Roche limit (Trevascus et al. 2021; Fortin-Archambault et al. 2020). Others investigate the effect of a primordial dust or gas disc, which is later perturbed by a disrupting body (Farihi et al. 2018; Kenyon & Bromley 2017; Malamud et al. 2021; Swan et al. 2020). Rafikov (2011) and Metzger et al. (2012) assumed that the gas is produced at a specific sublimation radius: dust particles are pulled into tighter orbits by the radiation of the host star (Poynting-Robertson drag), where they sublimate due to the higher temperatures. The authors modelled the formation of dust and gas discs simultaneously. However, they only considered disc formation in one dimension.

Multiple studies have found that in order for the disc morphology to be consistent with the observed structures, the disrupting body must have an eccentric orbit (Rafikov 2011; Metzger et al. 2012; Kenyon & Bromley 2017; Fortin-Archambault et al. 2020; Trevascus et al. 2021). As such, several studies have modelled a series of eccentric, precessing gas rings and have calculated synthetic Ca II emission lines using radiative transfer models (Gänsicke et al. 2006, 2008). Some have attempted to fit the observations with Monte Carlo methods (Goksu et al. 2024) or have assumed gas discs that are not in LTE (Hartmann et al. 2014, 2016).

Our aim is to develop a model that can simultaneously explain the anomalous metal pollution of WDs and the asymmetric Ca II lines over at least a 1.2-year time span while being consistent with the observations regarding line morphology and line asymmetry fluctuation periods. This is made possible through combining the results of two-dimensional hydrodynamic simulations with a simplified line radiative transfer model. The applied grid-based hydrodynamic model allows us to follow the evolution of gas discs over a timescale of 1.2 years (21 000 orbits), which has not been previously matched in the literature.

In our model, the gas discs originate from a disrupting asteroid in an eccentric orbit at the sublimation radius of the WD. Synthetic Ca II lines are then calculated using an analytic simplified line radiative transfer model. We give a detailed analysis of the morphology of the lines and investigate the persistence and precession timescales of the discs. We show that both the DZ nature of the WDs and the asymmetry of the Ca II lines can persist for over 1.2 years and that our model can replicate the observed precession periods and line morphologies.

This paper is structured as follows. Section 2 describes the numerical models used in this study and the examined parameter space. Our results are presented in Sect. 3, with each subsection addressing the effect of a single simulation parameter. Section 4 compares our results to the observed systems. Finally, our conclusions are given in Sect. 5.

2 Numerical methods

2.1 Hydrodynamic simulations

In the hydrodynamic simulations, an asteroid in an eccentric orbit loses material as its orbital distance falls below the sublimation radius of the WD. We used the grid-based, two-dimensional hydrodynamic code GFARGO2 to describe the gas dynamics. GFARGO2 is a GPU-based derivative of FARGO (Masset 2000), further developed in Regály (2020) and Regály et al. (2021).

This study focuses solely on gas disc evolution. The hydrodynamic differential equations used are the continuity equation, describing mass conservation, and the Navier-Stokes equation, describing momentum conservation, $$\frac{\mathrm{\partial }\mathrm{\Sigma }}{\mathrm{\partial }t}+\mathrm{\nabla }(\mathrm{\Sigma }\mathbf{v})={\dot{\mathrm{\Sigma }}}_{\mathrm{a}},$$(1) $$\frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }t}+(\mathbf{v}\mathrm{\nabla })\mathbf{v}=-\frac{1}{\mathrm{\Sigma }}\mathrm{\nabla }P+\mathrm{\nabla }\mathbf{T}-\mathrm{\nabla }\mathrm{\Phi }+\left({\mathbf{v}}_{\mathbf{a}}{\dot{\mathrm{\Sigma }}}_{\mathrm{a}}\right),$$(2)

where Σ and v are the surface density and the velocity vector of the gas, while Σ_a is the density added during the disruption of the asteroid and v_a is the velocity of the asteroid.

We calculated Σ_a at each time step of the hydrodynamic simulations in a way where material is added to the gas disc based on the mass loss rate of the asteroid derived from a disruption timescale, τ_dis, see later. The new material was distributed in a two-dimensional Gaussian distribution. The half width at half maximum of the distribution was calculated as (r_ah)/2, where r_a is the radial distance of the asteroid from the star and h is the aspect ratio of the disc. The density of the new matter was then added to the disc density (see Eq. (1)), and its momentum was calculated from the asteroid’s momentum (see Eq. (2)) at each time step. Over time, due to Keplerian shear and gas viscosity, the matter spread out to form an extended disc. The process of disc formation is discussed in detail in Sect. 3.1.

Assuming local thermodynamic equilibrium (LTE), the pressure in the disc can be expressed as $P=c_{\mathrm{s}}^2\mathrm{\Sigma }$, where c_s is the local sound speed. The gravitational potential (Φ) used in the Navier-Stokes equation is: $$\mathrm{\Phi }=-\frac{G{M}_{\mathrm{a}}}{\sqrt{r^2+r_{\mathrm{a}}^2-2r{r}_{\mathrm{a}}\cos (\varphi -{\varphi }_{\mathrm{a}})+(ϵH{)}^2}}-\frac{G{M}_{\ast }}{r}-{\mathrm{\Phi }}_{\text{ind\hspace{0.17em}}},$$(3)

where G is the gravitational constant, r, ϕ, and r_a, ϕ_a are the radial and azimuthal coordinates of a given grid point and that of the asteroid, respectively. M_* and M_a denote the mass of the WD and that of the asteroid, respectively.

The star was positioned at the centre of the grid, which does not coincide with the centre of the inertial frame. Therefore, the acceleration of the star must be taken into account as an indirect potential term, represented by Φ_ind in Eq. (3). Additionally, the gravitational potential of the asteroid must be smoothed to avoid a cell-asteroid distance singularity. This was achieved through the parameter ϵH in Eq. (3), adopting a value of 0.6 based on Kley et al. (2012). The viscous stress tensor is given by $$\mathbf{T}=\nu (\mathrm{\nabla }\mathbf{v}+\mathrm{\nabla }{\mathbf{v}}^T-\frac{2}{3}\mathrm{\nabla }\mathbf{v}\mathbf{I}),$$(4)

where I is the 2D identity matrix. Using the Shakura & Sunyaev (1973) α viscosity parameter, the kinematic viscosity can be expressed by ν = αc_sH. In the above expression H is the local pressure scale height, H = hr^1+γ, where h is the aspect ratio, and γ is the flaring index of the disc. For of geometrically thin discs, γ = 0, which means that H/r = h, where h is a constant.

Numerical algorithms in grid-based hydrodynamic codes are not suitable to describe areas with zero density. It is therefore necessary to model a background disc in the simulations. However, since the mass of the background disc was chosen to be one millionth of the mass of the asteroid, its properties did not affect the simulation results. During the simulations, material from the discs falls onto the star. This causes their mass to decrease over time. If the mass of a disc is less than that of the background, it is considered depleted.

At the inner and outer boundaries, an open boundary condition was applied in our simulations. This let the material flow out from the disc towards the star, but prevented it from flowing into the disc at the inner edge.

Given that the mass of the disc is less than that of the embedded asteroids, there isn’t enough torque acting on the asteroids to cause them to migrate (Veras et al. 2023). Therefore, we have disregarded the back-reaction of the gas on the asteroid, so it orbits with a constant semi-major axis.

During the simulations, we monitored the surface density (Σ), radial velocity (v_r), and azimuthal velocity (v_t) of the gas. Analysis of the models also included the examination of the temporal evolution of the asymmetry of the Ca II line, the disc mass, the total emitted flux and the accretion rate of the WD.

We calculated the Ca II emission line profiles by implementing a simplified line radiative transfer model, which is detailed in Appendix A. Our assumption that discs are in LTE and are optically thin means that the radial temperature profile in the disc cannot be correctly modelled using canonical blackbody radiation. In order to see the effect of an altered source function of blackbody radiation, we conducted a study where the Planck function is perturbed in either temperature or frequency space. The results of our analysis are presented in Appendix B. We conclude that perturbing the Planck function in either the temperature or the frequency domain has little to no effect on the Ca II lines. As such, we proceeded by assuming an unperturbed source function for blackbody radiation.

Line morphology is described by four parameters: degree of asymmetry (ξ, the difference in normalised flux between the two peaks of the line); central dip; shift or offset of line centre in velocity space (v_shift, the average distance of the two peaks from 0 km s⁻¹); and the half width of the line at base (HWB). These parameters are visualised in Fig. 1.

The temperature profile of the gas was given analytically. The blackbody temperature profile is defined as $$T_{\mathrm{B}\mathrm{B}}(r)={\left(\frac{R_{\ast }}{r}\right)}^{1/2}{T}_{\ast },$$(5)

where r is the radial distance from the star, while R_* and T_* are the radius and the effective temperature of the WD¹.

Our work assumed that the gas disc is formed by debris falling under the sublimation radius, which can be defined using the above-mentioned temperature profile. According to Kobayashi et al. (2011), dust is sublimated at 1500 K.

Computational limits based on the Roche radius of the WD and the temperature of the disc were used in the hydrodynamic modelling. The external computational limit was set at 5/3 times the Roche radius. The internal computational limits of the simulations in our model corresponded to the points at T = 5000 K in accordance with Hartmann et al. (2014). Figure 2 shows the process of disc formation.

The effects of general relativity are significant at distances less than 1000 times the Schwarzschild radius of the star (Sobolenko et al. 2017). However, the inner radius of the gas discs modelled in this paper is over 10 000 times the Schwarzschild radius of the WD. Therefore, in line with the studies of Manser et al. (2016b,a, 2019), the effects of general relativity in the simulations are neglected.

Fig. 1 Summary of the parameters describing line morphology: asymmetry ξ, central dip, shift of the line centre in velocity space (vshift) and half width at base (HWB).

Fig. 2 Process of disc formation. When the asteroid passes under the Roche limit, it loses material, which initially exists as dust until the asteroid passes the sublimation radius (T = 1500 K), within which gas is produced. The gas then spreads out viscously and forms a disc. The inner limit of Ca II emission is T = 5000 K. We note that the figure is not to scale.

2.2 Examined parameter space

The mass of the WD was set to M_* = 0.7 M_⊙, and its radius to R_* = 4.6 × 10⁻⁵ au in all models in accordance with the observations of Gänsicke et al. (2008) and Koester et al. (2014) The effective temperature of the star was set to 17000 K, which is consistent with the average temperatures derived by Gänsicke et al. (2007) and Melis et al. (2010). The mass of the disrupting asteroid was assumed to be M_a = 10⁻¹² Earth masses, or 6 × 10¹⁵ g, which corresponds to about the mass of comet 67P/Churyumov-Gerasimenko.

In our model the total mass of the asteroid was converted into gas. We note that the Ca II content is a free parameter in the simplified, optically thin line radiative transfer model used in this study. This means that the total emitted Ca II flux is directly proportional to the Ca II content. It is important, however, that the line morphology is independent of the injected mass as long as the gas is optically thin. The Ca II content can thus be modified in accordance to the presumed composition of the asteroid, which might be more complex, containing an iron core or even retaining a significant amount of water (Malamud & Perets 2016, 2017, a,b).

Assuming a stellar density of ρ_* = 3 × 10⁶ g cm⁻³ (Chandrasekhar 1994), and an asteroid bulk density of ρ_a = 3 g cm⁻³, the Roche limit is $$R_{\text{Roche\hspace{0.17em}}}=R_{\ast }{\left(3\frac{{\rho }_{\ast }}{{\rho }_{\mathrm{a}}}\right)}^{1/3}\approx 6.6\times {10}^{-3}\text{au}.$$(6)

The asymmetry and width of the Ca II emission lines increases with the inclination of the disc. When the disc is viewed edge-on (i = 90°), the line-of-sight velocities are maximal. Our study assumed an inclination of i = 45° for all models.

Using Eq. (5), the sublimation radius, below which asteroids lose mass, is located at R_sub = 6 × 10⁻³ au. In this case, given the star’s effective temperature of 17 000 K, the Roche and sublimation radii almost coincide.

The semi-major axis of the asteroids was chosen to always coincide with the sublimation radius, which resulted in an orbital period of 4.84 hours. The asteroids were modelled to be moderately eccentric, with e_a = 0.4 and 0.6. The semi-major axis of asteroid was kept constant. This meant that the time the asteroid spent below the sublimation radius varied depending on its eccentricity. Figure 3 illustrates the orbits of such asteroids.

In contrast to the simulations conducted by Trevascus et al. (2021), the asteroid’s disintegration was not assumed to be instantaneous. The chosen mass loss rate in one set of models ensured that if the asteroid were orbiting on a circular orbit, it would lose all of its mass during one revolution (τ_dis corresponds to 1 orbit). Other models investigated the effect of disruption that happens on a timescale equal to the length of the simulations (i.e. τ_dis corresponds to 1.2 yr ≈ 21 000 orbits). It is important to note that since mass loss only occured below R_sub, disruption on an eccentric orbit took longer than τ_dis would dictate.

Viscosity in the SPH simulations of Trevascus et al. (2021) depends on the number of particles modelled: higher particle numbers result in lower viscosity. However, modelling a large number of particles is computationally expensive, making SPH simulations unsuitable for modelling low viscosities (α ≤ 0.001). In contrast, the GFARGO2 code used in this work is grid-based and the viscosity is defined in the initial conditions, allowing a wider range to be investigated (0.001 ≤ α ≤ 0.1).

The aspect ratio of the disc was set to either h = 0.05 or h = 0.02 in our models, the latter of which coincides with the value adopted by Trevascus et al. (2021). The discs modelled were assumed to be geometrically thin (γ = 0), in accordance with the assumptions of Jura (2003) and Veras (2016). The parameters of all models are summarised in Table 1.

Fig. 3 Arrangement of asteroid orbits (marked by red lines) at different eccentricities relative to the sublimation radius (white dashed lines). We note that eccentric asteroids spend only part of their orbits below the sublimation radius.

3 Results

3.1 Disc formation

Let us start by making some general statements about the process of disc formation by analysing the model B1 (see parameters in Table 1). Three phases can be identified during disc formation. These are shown in Fig. 4, where each Panel displays the surface density of the disc at different times. During the first phase, the gas is pulled apart in a spiral tail along the orbit of the disintegrating asteroid by the Keplerian shear. In the second phase, the gas forms a ring that begins to precess in the retrograde direction. At this point, the Ca II triplet exhibits asymmetry due to the asymmetric density distribution (the effect of the density and velocity distributions is discussed further in Sect. 3.2). Emphasise that the first and second phases are transient phenomena, meaning that the chance to observe a disc during these phases is low. However, the second phase can last longer if the viscosity of the disc is lower (see a further discussion in in Sect. 3.3). In the third and longest phase, the gas ring viscously spreads to form a disc and starts to accrete onto the WD. If the asteroid had an eccentric orbit, the resulting disc will also be eccentric. Eventually the disc depletes and reaches a mass that is equal to the background disc mass. We note that there is a continuous transition between each phase.

Changes of the disc density and velocity distributions in phase three of disc formation are shown in Fig. 5 along with the calculated background-normalised Ca II line profiles. The asymmetry of line peaks change on the timescale of the period of the retrograde precession of the disc. The greatest asymmetry can be seen when viewing the disc from the major axis and the asymmetry is absent when viewing from the minor axis. The 1.2-year evolution of the Ca II line profiles in each examined model are shown on spectrograms in Appendix C.

Steep wings are common to all of the calculated Ca II lines. This is due to the fact that at the internal computational boundary (i.e. at T = 5000 K) Ca II is converted into Ca III, as noted by Gänsicke et al. (2006). Therefore, Ca II is not observable within this limit.

Fig. 4 Phases of disc formation: snapshots of the density distribution in model B1. Grey dashed lines mark the sublimation radius. A black arrow shows the orbital direction of the asteroid.

3.2 Origins of the asymmetric Ca II lines

In this section we analyse the line profile of model B1 in two cases: a) when only the density distribution of the hydrodynamic simulations is used when calculating the Ca II emission, while the velocity distribution is assumed to be Keplerian²; and (b) when both the density and velocity distributions are extracted from the hydrodynamic simulations.

The distribution of the line-of-sight velocity of the gas in the two cases is shown in the top two Panels of Fig. 6. The left Panel shows the Keplerian velocity distribution, while the right Panel shows the one extracted from the hydrodynamic simulation. In both cases, the density distribution is the same as the one in Fig. 5 at t = 132.3 days. It is clear that the line-of-sight velocity distribution is symmetric in the Keplerian case. In contrast, the velocity distribution of the hydrodynamic simulation is perturbed.

The line profiles calculated for each case are shown in the lower Panel of Fig. 6. Despite the symmetry of the Keplerian case, the Ca II profile still shows some degree of asymmetry at ξ = 0.06. Here, the asymmetric density distribution resulting from the asymmetric material distribution of the eccentric disc is responsible for the asymmetry. On the other hand, if the velocity distribution is also non-Keplerian, the asymmetry of the Ca II line becomes significant, with ξ = 0.22.

For a Keplerian velocity distribution, the centre of the line is located at 0 km s⁻¹ in all cases. However, if the velocity distribution is also non-Keplerian, the line centre is shifted (refer to the vertical dashed and solid grey lines in Fig. 6). If an offset asymmetric line profile is detected, it can be concluded that the gas disc is eccentric and the velocity distribution is asymmetric (Regály et al. 2011, 2014). In the model discussed here, the line centre is located at 26.8 km s⁻¹, which corresponds to an offset of 0.74 nm in wavelength space.

The width of the Ca II lines and location of the line centres are also dependent on the inclination of the gas discs. Figure 7 displays a line profile of model B1 at various angles of inclination (t = 132.3 days). The line becomes narrower and less offset as the inclination decreases. At an inclination of i = 90° (i.e. looking at the disc edge-on), the asymmetry is at its maximum. Additionally, in this case, the line is as wide as possible since the line-of-sight velocities are maximal. The maximum value of the half-width at bottom (HWB) is 1170 km s⁻¹. The highest possible offset is v_shift = 301 km s⁻¹. We note that in our models the discs are optically thin even when they are viewed edge-on.

Fig. 5 Density and velocity distributions (Σ and vlos, respectively) in model B 1 at three points in time (top row) and the corresponding background-normalised Ca II lines (bottom row). The discs are viewed from below. The sublimation radius is marked by a black dashed line, and the orbital direction of the planet is shown with a black arrow. In the Panels showing the Ca II lines, vertical dashed lines indicate the centre of the line. The retrograde precession of the disc causes a periodic shift in the intensity maximum.

3.3 Disc viscosity

This section examines the evolution of discs and the Ca II line profile while using different α viscosity parameters. Values of α = 0.05, 0.01 and 0.001 (models B3, B1 and B2, respectively) are analysed.

We can estimate the depletion timescale of a disc by performing a log-linear fit on the disc mass data, and extrapolating to the point in time when the disc mass is equal to that of the background. Extrapolation of the data of model B2 yields a depletion timescale of 13.9 years. However, due to its small viscosity (α = 0.001), this model evolves very slowly, and fails to develop a disc structure even after 1.2 years. Furthermore, by the time the disc structure can fully evolve, the disc circularises. This means that the Ca II line profiles will only be asymmetric in the first and second phase of disc formation (see Sect. 3.1). This model has not been considered further.

In the α = 0.01 and 0.05 models, the smaller the viscosity, the weaker the asymmetry: ξ = 0.13 in the α = 0.01 model, and ξ = 0.45 in the α = 0.05 model.

The evolution of the disc mass and the flux are correlated: as the discs start to deplete, the emitted flux also drops by orders of magnitude. This means that the depletion of the disc also corresponds to the end of the observability for the Ca II triplet. As such, measuring the degree of asymmetry and the other line parameters (dip, v_shift and HWB) is only meaningful in the period of best observability, that is when the disc is fully formed and its mass is orders of magnitude higher than that of the background disc.

The depletion timescale also depends significantly on the viscosity. The mass of the discs decreases more slowly if the viscosity is smaller. The depletion timescales are estimated to be 1.2 years (α = 0.05) and 4.3 years (α = 0.01).

This is also reflected in the total emitted flux: while the flux in the α = 0.01 model is saturated at ≃10²⁶ erg s⁻¹ cm⁻² for the whole runtime, the flux in the α = 0.05 model shows a constant decrease, and is only ≃10²² erg s⁻¹ cm⁻² at the end of the simulation.

The same is true for the accretion rate of the WD. In model B 1, it is constant with a value of ≃10⁸ g s⁻¹ throughout the 1.2 years of the simulation. However, in model B3 it drops from 10⁸ g s⁻¹ to 10⁵ g s⁻¹ in just 0.9 years.

The retrograde disc precession causes a periodic intensity change in the line profile. To determine the precession period τ_p, the Fourier transform of the asymmetry time series is analysed. Again, this is only meaningful to calculate when the disc is fully formed, but still has a significant mass (10⁻³ M_a).

Reducing viscosity not only extends asymmetry and disc lifetime but also significantly influences the precession period of the discs. The precession periods of the discs analysed in this section are found to be 27.2 and 10.6 days for models B1 and B3, respectively.

For the models considered, the kinematic viscosity and, therefore, the viscous timescale differ. For geometrically flat discs in the Shakura & Sunyaev α-viscosity prescription, the kinematic viscosity is defined as $$\nu =\alpha h^2{a}^{1/2},$$(7)

which means that the kinematic viscosity in the models considered is ν = 1.94 × 10⁻⁶ m² s⁻¹.

The viscous timescale for each disc at the radius corresponding to the semi-major axis of the asteroid, r_a, is $$t_{\nu }=\frac{r_{\mathrm{a}}^2}{\nu }\frac{\sqrt{M_{\ast }}}{2\pi }.$$(8)

Table 2 displays the viscous timescales and precession periods in each model.

Fig. 6 Comparison of a Keplerian disc and model B1. The upper panel displays the line-of-sight velocity distribution in a Keplerian disc as well as that of model B1 from the hydrodynamic simulations (t = 132.3 days). The sublimation radius is indicated by black dashed lines. The lower panel shows the calculated Ca II line profiles with dashed and solid lines in the two cases. The line centres are marked by vertical grey lines.

Fig. 7 Ca II line profiles for model B1 plotted with different disc inclinations. Fainter lines represent larger inclinations.

3.4 Disc geometry

In this section, model B4 is analysed, where the aspect ratio of the disc is h = 0.02. Its viscosity is α = 0.01.

The model shows significant line asymmetry, with a maximum of ξ = 0.48. The disc precesses at a slower rate than any other model presented in this study, the precession period being 177.4 days.

The decrease in disc mass is slow enough for the Ca II signal to remain observable after 1.2 years in both models. The depletion timescales of the disc is found to be 9.5 years. Comparing these values with model B1 (where h = 0.05 and the depletion time is 4.3 years), it is evident that the aspect ratio of the discs has a strong effect on the disc lifetime. The lower aspect ratio of h = 0.002 results in a disc that persist for about twice as long as their high aspect ratio counterparts.

Regarding the total flux, the model outputs a total flux of ≃5 × 10²⁶ erg s⁻¹ cm⁻². The accretion rate of the WD saturates to ≃5 × 10⁸ g s⁻¹ and is constant during the simulation. This accretion rate is different from the B1 model, where h = 0.05, because in the Shakura & Sunyaev α-prescription, the kinematic viscosity depends on the aspect ratio of the disc (see Eq. (7)).

3.5 Asteroid disruption rate

In this section, we discuss models where the disintegration of the asteroid occurs in 1.2 years. Three such models are compared: B5, B6 and B7.

Disc mass is growing throughout all simulations. This means that the depletion timescale is longer than the simulation time, i.e, all discs persist for at least 1.2 years. The only time we are able to observe persisting asymmetry is if α = 0.1. For the models that use α = 0.01, the asymmetry declines swiftly, and is completely gone in about 0.5 years.

With regard to the α = 0.1 model (B5), the maximum degree of asymmetry is ξ = 0.13. The precession period is found to be 126.7 days. At the end of the simulation, disc mass is saturated to ≃10⁻² M_a, the total emitted flux to ≃10²⁶ erg s⁻¹ cm⁻², and the accretion rate of the WD to ≃10⁷ g s⁻¹.

Models B6 and B7 were run with different asteroid eccentricities. However, they yield the exact same results regarding the monitored parameters. This means that the initial eccentricity of the asteroid has no effect on the long-term morphology of the lines, or the persistence of the gas disc.

4 Discussion

4.1 Comparison to observations

According to observations, the accretion rate of the DZ WDs that host gas discs fall in the range of 10⁵−10¹² g s⁻¹ (Koester & Wilken 2006; Koester et al. 2014). For the models presented here, accretion rates of 10⁵ − 5 × 10⁸ g s⁻¹ are derived from the hydrodynamic simulations during the periods of best observability. It is important to note that the accretion rate is directly proportional to the mass of the asteroid and the above-listed values can be scaled without altering the morphology of the Ca II lines.

Trevascus et al. (2021) assumes that Ca II line asymmetry develops in the second phase of disc formation, when material is in a ring around the WD. While the Ca II lines are indeed asymmetric in this phase, the gas does not reach the atmosphere of the WD until the third phase, when the disc is fully developed (see Sect. 3.1). This means that the DZ nature of the atmosphere and accretion can only be simultaneously observed in the third phase of disc formation.

The observed asymmetric Ca II lines have a HWB of 130–1150 km s⁻¹, as determined by the Doppler shift of the peaks of the emission lines (Dennihy et al. 2020; Gänsicke et al. 2006, 2008; Gentile Fusillo et al. 2021; Manser et al. 2016b; Melis et al. 2010, 2012, 2020). It is important to note that these observations are degenerate with the inclination of the discs. In the simulations presented in this work, assuming the disc to be seen as edge-on, the maximum HWB of the Ca II line is 1170 km s⁻¹. For WDJ0234-0406 (Gentile Fusillo et al. 2021), the HWB of the observed lines is 1150 km s⁻¹, which is only consistent with the simulations if the disc is seen close to edge-on (i ≈ 90°). The HWB measured in each model, along with other line parameters, are presented in Table 3.

As discussed in Sect. 3.2, the centre of the lines shifts when the velocity distribution is asymmetric. The maximum shift, which occurs for i = 90°, is 301 km s⁻¹ for all models. These highest offsets can be observed in model B4, whose aspect ratio is only 0.02, as opposed to the value of 0.05 in other models (see Table 3). The above-mentioned offset values are consistent with those of Melis et al. (2010), who reported line centre displacements of 17–104 km s⁻¹ in three observed systems, and Dennihy et al. (2020), who measured a shift of up to 260 km s⁻¹ of the Ca II emission in the HE 1349-2305 system.

The average central dip of the lines ranges from 0.43 to 0.54 (see Table 3). Regarding previous observations, the central dip of 28 (by eye fit) examined Ca II lines is between 0.13 and 0.8, with about 17% of which lie within the 0.43 to 0.54 range. It should be noted the values of the central dip shown in Table 3 are an average over the period of best observability, and can be 0.1 larger or smaller at a given time. Taking this into account, the models presented here are able to reproduce the central dips seen in 71% of observations.

The central dip of the Ca II triplet lines is both dependent on the inclination and the optical depth of the disc. As shown in Fig. 7, the dip is larger for a larger inclination. For sufficiently large inclinations, the disc may become optically thick (Horne & Marsh 1986). However, this is not the case in the models presented here, as the discs would only become optically thick if their masses were significantly greater (by about 100–1000 times). Therefore, the observed dip values that differ from those presented here can be explained by a different inclination than that adopted in our models (45°).

The line profiles have steep wings in all of our models, and the Ca II flux is zero in the T > 5000 K regions of the discs. This agrees with the assumption of Gänsicke et al. (2006) that the steep wings are not caused by a complete absence of material, but by additional ionisation of Ca II into Ca III. This is supported by the fact that all WDs with discs have a DZ atmosphere: the gas does reach the star’s surface but is no longer present as Ca II.

The observed average extent of the gas disc around WDs is 2.35 × 10⁻³−5.64 × 10⁻³ au (Gänsicke et al. 2008; Farihi et al. 2016). For the simulations in this study, the disc extent is between 0.0005 and 1 × 10⁻² au. These models thus produce discs with larger extents than those observed thus far.

The disruption rate of the asteroid also plays a role in determining the total emitted flux: for models where the asteroid disrupts in one orbit, the flux is about five times higher than if it were disintegrating in 1.2 years. This is a direct consequence of the optically thin LRT model, in which the emitted flux is directly proportional to the disc mass.

The models in this work assume a stellar temperature of 17 000 K. Observations suggest that DZ WDs that host a gas disc have effective temperatures ranging from 13 300 K to 29 000 K (Farihi et al. 2012; Gentile Fusillo et al. 2021). Figure 8 shows the Ca II line profile assuming a stellar temperature of 13 300, 17 000 and 29 000 K. If the stellar temperature is larger, the inner edge of the disc (where T = 5000 K) will be farther from the WD. This implies that the maximal line-of-sight velocities, and thus the HWB of the Ca II lines, will also be reduced. The asymmetry of the lines decreases, and can even disappear completely for substantially large temperatures, while the central dip decreases with stellar temperature. We note that for cooler WDs, the extent of the gas disc is smaller. In such a case it is necessary to explain how the debris created by the disrupting asteroid on the Roche radius can reach the sublimation radius. Conversely, around hot stars, the sublimation radius may extend beyond the Roche limit, causing the dust and gas discs to spatially overlap.

Figure 9 shows the spatial arrangement of three dust and gas discs presented in Melis et al. (2010), the extents of which are very similar. Additionally, the temperature profile investigated in this study is plotted, assuming two stellar temperatures, their corresponding sublimation radii (where T = 1500 K), and the location of the Roche limit. Compared to our models, the dust discs appear to be inside the sublimation radius. The dichotomy presented in Fig. 9 can be resolved if the inclination of the gas disc differs significantly from those obtained by Melis et al. (2010).

The asymmetry of the Ca II lines has been observed for over twelve years (Manser et al. 2016b). To maintain this asymmetry for such a long period, it is necessary for the disc to also persist. According to our simulations, this is possible in two cases: either when the viscosity is α ≤ 0.01 and the asteroid disrupts in one orbit, or if α = 0.1, but the disruption of the asteroid is prolonged to 1.2 years. Table 3 summarises the depletion timescale τ₀ for each model, obtained by fitting a log-linear function to the temporal evolution of the disc mass.

The fading of the Ca II lines has only been observed in the SDSS J1617+1620 system (Wilson et al. 2014). The flux increased significantly in 2008 and then showed a steady decrease until 2014 when it disappeared completely. Similarly, Gentile Fusillo et al. (2021) showed a 50% decrease in the flux of Ca II lines over a three-year period. These observations support the idea that disc depletion can occur on a timescale of a few years, resulting in a decrease in the flux of Ca II lines below what is observable. Extrapolating the simulations of the first 1.2 years shown in this study, the disc in models B1, B3 and B4 are estimated to deplete in 4.3 and 1.3, and 9.5 years, respectively (see Table 3). These models, therefore, show a rough agreement with the two observations mentioned above.

Based on the periodic fluctuation in the asymmetry of the Ca II lines, previous research has shown that the precession timescales of the observed gas discs range from a few hours (123 minutes, Manser et al. 2019), to a few months (120 days, Dennihy et al. 2020), and even several years (1.4–37 years, Dennihy et al. 2018; Manser et al. 2016b; Fortin-Archambault et al. 2020; Goksu et al. 2024; Hartmann et al. 2016). The precession period of the models analysed in this paper ranges from 27.2 days to 177.4 days (all precession periods are listed in Table 3). The precession periods of a few months found by Dennihy et al. (2018, 2020) can thus be explained by these models. However, they cannot account for the two-hour precession periods measured by Manser et al. (2019) or the years- and decades-long precession periods found by Hartmann et al. (2016), Manser et al. (2016b) and Goksu et al. (2024). It should be noted that poor temporal resolution of measurements may result in precession periods appearing to be several years or decades long. Many of the above-mentioned studies had a time gap of months, and in some cases, a year or two, between two measurements. Therefore, if the discs precess at a rate faster than this, the measurements will not be able to detect the rapid variation.

According to the findings presented in Sect. 3.3, a viscosity of α ≤ 0.01 is required to observe asymmetry for a minimum of 1.2 years. This is inconsistent with the assumption of Trevascus et al. (2021), where α = 0.05. However, Trevascus et al. (2021) only followed the evolution of the discs for 8.61 days, so no claims can be made about the persistence of the asymmetry over time. In our case, model B2 (where α = 0.001) suggests that the viscosity cannot be reduced to any arbitrary value. The lower the viscosity, the longer the disc remains in the ring phase; however, if the viscosity is too low, the disc will circularise before its material reaches the atmosphere of the WD. Therefore, it is not possible to simultaneously explain accretion and Ca II line asymmetry in this scenario.

In the work of Trevascus et al. (2021), the asteroid has a mass of 1.6 × 10⁻⁶ Earth masses (0.1 times the mass of Ceres). In this work, the asteroid’s mass is assumed to be 10⁻¹² Earth masses, and the disc is assumed to be optically thin in all simulations. If the asteroid’s mass were to be reduced, the absolute flux of the Ca II lines would decrease proportionally. However, the mass of the asteroid does not affect the strength or timescale of the asymmetry as long as the disc is optically thin. The discs would only become optically thick if the asteroid’s mass were increased above 10⁻¹⁰−10⁻⁹ Earth masses.

The models presented in this work demonstrate that the Ca II emission signal can be detected over multiple precession periods if the asteroid’s disruption occurs on a timescale commensurable to the length of the observations. This is in a rough agreement with the calculations of Manser et al. (2019), which show that a rocky body with a mass of ≤10⁻⁹ Earth mass orbiting at the Roche radius of the WD can survive for 1.5 years before it completely disrupts. Malamud & Perets (2020b,a) also reach a similar conclusion: the disruption of a body at a WD’s Roche limit does not occur in a single orbit. The latter is consistent with the models presented in this work, since the asteroids with eccentric orbits also cross the sublimation radius several times.

The presence of minor planets orbiting DZ WDs (Vanderburg et al. 2015; Manser et al. 2019; Vanderbosch et al. 2020; Guidry et al. 2021) supports the validity of the model presented in this paper. Around DZ WDs, the asteroid and the gas disc, and thus the Ca II emission, are observed simultaneously. This implies that asteroid disintegration and gas production are simultaneous, that is disruption does not precede gas formation.

The presence of a second or third body also raises an interesting question. A giant planet can excite the eccentricity of its companion, causing it to enter the Roche or sublimation radius and triggering the formation of dust or gas. The reason why we do not see these companion giant planets in the DZ WD systems with transiting exoplanets is likely due to the small size of the star. Indeed, a giant planet orbiting a WD at 6 × 10⁻³ au will not transit for i ≥ 0.44°.

Fig. 8 Ca II line profiles assuming different stellar temperatures in model B1. The snapshots depict the largest possible red-side asymmetry.

Fig. 9 Spatial extent of three gas and dust discs that are detected by Melis et al. (2010). The dust discs are marked by grey rectangles, while the gas discs are marked by orange rectangles. To improve visibility, the discs for each system are offset vertically. For the two stellar effective temperatures used in the Melis et al. (2010) study, namely 18 500 K (red lines) and 22 000 K (blue lines), the blackbody temperature profile (light red and blue lines) is shown. A horizontal green line indicates the sublimation temperature (1500 K). The figure displays the sublimation radius for each configuration with vertical dashed lines. The dashed grey line indicates the radius of the WD, while the solid black line represents the Roche limit.

4.2 Caveats and limitations

The hydrodynamic model presented in this paper only follows the formation of gas discs. It does not take into account the asteroid’s path to its orbit or the process of disintegration due to the significant computational effort required to describe the multi-component dust and gas discs and their interaction. Therefore the presented model cannot explain the existence of systems where only dust discs are observed. Nevertheless, several possibilities exist to resolve this: firstly, future measurement tools may be sensitive enough to detect lower-mass gas discs. Secondly, it is possible that a gas disc existed in these systems at an earlier time, but viscous spreading caused it to move far away from the star, rendering them undetectable by the Ca II signal. It is also possible that the depletion of gas discs does not necessarily mean the depletion of dust discs.

If the Roche limit is outside the sublimation radius, the dust produced on the Roche radius can only approach the WD to a limited extent due to its radiative pressure and the Yarkovsky or YORP effect. To observe only a dust disc, there must be no dust present under the sublimation radius. Therefore, a mechanism is required to increase the semi-major axis of the dust particles, which can only be achieved by the Yarkovsky effect assuming a retrograde particle rotation. It is worth noting that the anomalous metal content of WDs can still be observed for up to thousands to even millions of years, even after the discs have been entirely depleted.

The GFARGO2 code calculates the surface density and velocity of discs in two dimensions. However, to gain a better understanding of the time evolution of gas discs, three-dimensional studies are necessary. This would enable the modelling of vertical flows in the disc and the three-dimensional nature of WD accretion (funnel flow, as described in Romanova et al. 2003).

In our model, the internal computational boundary is set at the point where the disc temperature reached 5000 K. In this case, the process of accretion cannot be satisfactorily modelled.

The simulations presented here do not model the stellar magnetic field. Instead, its effect is parametrised using the Shakura & Sunyaev α-viscosity. To accurately describe reality, magneto-hydrodynamic simulations would be necessary. However, their computational complexity exceeds that of the hydrodynamic simulations discussed in this work.

In this investigation, the gas discs around WDs are assumed to be optically thin. However, there are observations in the literature that suggest discs are optically thick (Gänsicke et al. 2006; Jura 2003). Their conclusion is supported by the fact that the centre dip of the Ca II lines is deeper than expected for an optically thin disc (Horne & Marsh 1986). As the Ca II triplet is seen in emission, the emitting region must be optically thin. However, there can still be an underlying, optically thick layer of gas underneath that has a lower temperature compared to the atmosphere.

Our models assume that the discs are in LTE. Establishing a fully non-LTE model means that the level population of the lines must be derived by taking into account that a) excitation temperatures may deviate from gas (kinematic) temperature and b) level populations may be strongly radiatively self-coupled. Furthermore, the temperature dependence of the equations of state must also be taken into account. This requires solving the energy equation in addition to the continuity and Navier-Stokes equations. The necessary numerical methods for this latter, computationally expensive task have been implemented in GFARGO2 (see Tarczay-Nehéz et al. 2020). However, including such complexity is beyond the scope of the current paper.

5 Conclusions

This study examines the evolution of gas discs around WDs and the resulting Ca II emission. Gas discs are assumed to be formed by the disruption of asteroids in eccentric orbits at the sublimation radius, which is situated below the Roche-radius of the WD. The long-term observability of the Ca II line asymmetry and the metal pollution of the WD are monitored simultaneously.

Hydrodynamic simulations in a locally isothermal regime are carried out on a two-dimensional, cylindrically symmetric grid to model the evolution of the density and velocity distribution of gas discs. The results of the hydrodynamic simulations are used to calculate synthetic Ca II emission lines with a simplified line radiative transfer model. The analysis of the results includes the monitoring of the 1.2-year evolution of the disc mass, the degree of the Ca II asymmetry, the emitted flux, the accretion rate of the WD and the precession period and depletion timescale of the discs. The main results of the simulations are summarised below.

1. Asteroids disrupting in eccentric orbits produce globally eccentric, retrogradely precessing discs in three phases: 1) spiral tail; 2) ring; and 3) disc. The asymmetry of the Ca II line is explained by the asymmetric density and velocity distribution of the disc. Gas eventually falls onto the WD, explaining its anomalous metallicity.

2. Our results show that the asymmetric Ca II line profile and the anomalous metallicity of the WD can only be observed simultaneously during the third phase of disc formation. This finding resolves the discrepancy of previous results that assume that the Ca II line asymmetry originates in the second phase of disc formation, when the gas is distributed along a ring and does not reach the surface of the WD.

3. The centre of the Ca II lines is shifted (by up to 301 km s⁻¹) if the disc velocity distribution is non-Keplerian. The shift is a direct indicator of disc eccentricity.

4. For a moderate inclination of 45°, the asymmetry of the Ca II triplet lines ranges from ξ = 0.13 to ξ = 0.45. The width of the lines is approximately 836 km s⁻¹. The central dip of the lines ranges, on average, from 0.43 to 0.54. The magnitude of the asteroid’s eccentricity does not affect the morphology of the lines significantly.

5. There are two ways to explain accretion and asymmetry simultaneously on a timescale of at least 1.2 years: (a) the asteroid disrupts in one orbit, and disc viscosity is 0.001 < α < 0.05; (b) the asteroid disrupts in a time span of 1.2 years, but the disc viscosity high, in the order of α = 0.1

6. Disc precession causes a periodic change of the red and blue peaks of the Ca II lines whose period falls in the range of 27.2 to 177.4 days.

In summary, our work suggests that the persistence of Ca II asymmetry over decades and its periodic change in the peaks can be explained by asteroids on eccentric orbits in two scenarios. In the first case, the asteroid disrupts on a short timescale (couple of orbits), and the gas has a low viscosity range (0.001 < α < 0.05) to maintain the Ca II signal for decades. In the other scenario, the asteroid disrupts on a timescale of a year, and the viscosity of the gas is required to be high, α = 0.1.

The precise origin of WD pollution remains a topic of debate. However, a number of models have been proposed that suggest that small bodies may be scattered or evolve into highly eccentric orbits over time (see, for example, Veras et al. 2020). Moreover, a highly eccentric, comet-like body has been observed to orbit DZ WDs (Vanderbosch et al. 2020). It is therefore important to explore this scenario further in our next study, with a model that considers the disruption of bodies on highly eccentric orbits (e ≳ 0.9). Our preliminary findings indicate that in this case, the gas discs must experience a depletion between the pericentre passages of the comet in order for the disc to maintain its eccentricity over a year-long timescale.

The density distributions provided by the hydrodynamic simulations can be used to calculate the transit light curves of slowly disintegrating asteroids, which can be compared to the asymmetric transit light curves of observed objects. To investigate the nature of the transit asymmetry, implementing the effect of the stellar radiation pressure into GFARGO2 appears promising.

Acknowledgements

We thank the anonymous referee, whose helpful comments improved the quality of the paper. The authors would like to express their gratitude to J. Vinkó for fruitful discussions on the nature of line radiative transfer in the discs under consideration, which discussions have had a significant impact on the quality of our study. VF thanks for the financial support from the ÚNKP-23-2 New National Excellence Program of the Ministry for Culture and Innovation from the source of the National Research, Development and Innovation Fund. VF acknowledges the financial support of the ’SeismoLab’ KKP-137523 Élvonal grant of the Hungarian Research, Development and Innovation Office (NKFIH) and the undergraduate research assistant program of Konkoly Observatory.

Appendix A Simplified line radiative transfer model

The Ca II emission line profiles are calculated via a numerical implementation of a simplified radiative transfer model. The density and velocity distributions obtained from the two-dimensional hydrodynamic simulations are used for the calculations. The line radiative transfer model assumes that the discs are optically thin.

If the medium is both emitting and absorbing along the propagation of radiation s, the basic equation for radiative transfer is: $$\frac{d{I}_{\nu }(s)}{ds}=j_{\nu }(s)-{\chi }_{\nu }(s)I_{\nu }(s),$$(A.1)

where ds is the infinitesimal path length along the propagation of the radiation, I_v(s) is the monochromatic intensity at frequency v, and χ_ν(s) and j_ν(s) are the absorption and emission coefficients of the medium at the given frequency. The monochromatic optical depth along the beam of radiation is τ_ν(s) = χ_ν(s)ds. As such, that of a medium with a total thickness of D is $${\tau }_{\nu }={\int }_0^D{\chi }_{\nu }(s)ds.$$(A.2)

With the above, the radiative transport equation becomes $$\frac{d{I}_{\nu }\left({\tau }_{\nu }\right)}{d{\tau }_{\nu }}=\frac{j_{\nu }\left({\tau }_{\nu }\right)}{{\chi }_{\nu }\left({\tau }_{\nu }\right)}-I_{\nu }\left({\tau }_{\nu }\right)=S_{\nu }\left({\tau }_{\nu }\right)-I_{\nu }\left({\tau }_{\nu }\right),$$(A.3)

where S_ν(τ_ν) is the source function of the emission. The solution of this equation describes the propagation of a radiation with frequency v in a medium with an optical depth of τ_ν. If I_ν(0) is defined as the intensity at the source, it is found in the form of $$I_{\nu }\left({\tau }_{\nu }\right)=I_{\nu }(0)e^{-{\tau }_{\nu }}+{\int }_0^{{\tau }_{\nu }}{S}_{\nu }\left({\tau }^{\mathrm{\prime }}\right)e^{-({\tau }_{\nu }-{\tau }^{\mathrm{\prime }})}d{\tau }^{\mathrm{\prime }}.$$(A.4)

The absorption in a medium with a density of ρ(s) can be given as χ_ν(s) = ρ(s)κ_ν, where κ_ν is the monochromatic mass absorption coefficient of the medium, assumed to be constant along the path of the radiation, namely κ_ν(s) = κ_ν. Now Eq. (A.2) becomes $${\tau }_{\nu }={\int }_0^D\rho (s){\kappa }_{\nu }ds={\kappa }_{\nu }{\int }_0^D\rho (s)={\kappa }_{\nu }\mathrm{\Sigma }\frac{1}{\cos (i)},$$(A.5)

where Σ is the density integrated along the path of the radiation, interpreted as the surface density of the disc in the 2 D model. The inclination of the disc is denoted by i, and i = 90° corresponds to an edge-on disc.

Assuming that the investigated medium radiates as a blackbody, its radiation is described by the Planck function, $$B_{\nu }(T)=\frac{2h{\nu }^3}{c^2}{[\exp \left(\frac{h\nu }{kT}\right)-1]}^{-1},$$(A.6)

where k and h are the Boltzmann, and Planck constants, respectively, and c is the speed of light in vacuum. It should be noted that the blackbody presumption might be questioned in an optically thin disc in the LTE state. Nevertheless, modifying the Planck function in both temperature and frequency domains results in a line morphology that is nearly indistinguishable from the one obtained using the canonical source function for blackbody radiation. For further details, refer to Appendix B. If the medium is assumed to be in thermal equilibrium, the energy that it emits and that it absorbs must be equal. In this case dI_ν(τ_ν)/dτ_ν = 0, and applying Eq. (A.3) we get Kirchhoff’s law, S_ν(τ_ν) = I_ν(τ_ν) = B_ν(T). This means that the source function in Eq. (A.4) can be replaced by the Planck function, and as such, the intensity is $$I_{\nu }=I_{\nu }(0)\exp (-{\kappa }_{\nu }\mathrm{\Sigma })+B_{\nu }(T)(1-\exp (-{\kappa }_{\nu }\mathrm{\Sigma })).$$(A.7)

This equation describes that the originally emitted intensity is increased due to the emission of the medium (second term), while it is decreased due to the absorption of the medium (third term).

To solve Eq. (A.7), we need to know the monochromatic mass absorption coefficient, κ_ν. For a simple atomic model with two energy levels, the energy of the lower level is ϵ₁, and that of the upper level is ϵ_u = ϵ₁ + hν₀. This atom can partake in three different radiation processes. Firstly, it can spontaneously emit a photon with a frequency of ν₀. The transition probability for spontaneous emission per unit time is given by the Einstein-coefficient A_ul. Next, it can absorb a photon with a frequency of ν₀. The probability of this process is proportional to the mean intensity at ν₀, $\overline{J}$, defined as $$\overline{J}=\frac{1}{4\pi }{\int }_0^{\mathrm{\infty }}{\int }_0^{2\pi }{\int }_0^{\pi }{I}_{\nu }\sin (\theta )d\theta d\phi {\varphi }_{\nu }d\nu .$$(A.8)

The probability of the transition per unit time is $B_{\mathrm{l}\mathrm{u}}\overline{J}$, where B_lu is another Einstein-coefficient. The line profile function Φ_ν in Eq. (A.8) describes the relative effectiveness of neighbouring frequencies of ν₀ for causing a transition. The third possible radiation process is stimulated emission. In this case, the transition probability is again proportional to the mean intensity and is given by $B_{\mathrm{u}\mathrm{l}}\overline{J}$. The record of the three radiation processes in LTE, assuming that the number of atoms in the lower state and in the upper state are n_l and n_u, respectively, is $$n_{\mathrm{u}}{A}_{\mathrm{u}\mathrm{l}}=n_{\mathrm{u}}{B}_{\mathrm{l}\mathrm{u}}\overline{J}-n_1{B}_{\mathrm{u}\mathrm{l}}\overline{J}.$$(A.9)

The ratio of the number of species in the upper and lower level states can be given by $$\frac{n_{\mathrm{u}}}{n_{\text{l}}}=\frac{g_{\text{l}}\exp (-{\epsilon }_{\text{l}}/kT)}{g_{\mathrm{u}}\exp (-({\epsilon }_{\text{l}}+h{v}_0)/kT)}=\frac{g_{\mathrm{u}}}{g_{\text{l}}}\exp (-\frac{{\epsilon }_{\mathrm{u}}-{\epsilon }_{\text{l}}}{kT}),$$(A.10)

where g_l and g_u are the degeneracy of the lower and upper states, respectively. Using Eqs. (A.9) and (A.10), and due to the fact that in LTE $\overline{J}=B_{\nu }(T)$ must be valid for all temperatures, the Einstein coefficients (A_ul, B_ul and B_lu) are coupled by the following relations: $$B_{\mathrm{u}\mathrm{l}}=\frac{c^2}{2h{\nu }_0^3}{A}_{\mathrm{u}\mathrm{l}},B_{\mathrm{l}\mathrm{u}}=\frac{g_{\mathrm{u}}}{g_{\mathrm{l}}}{B}_{\mathrm{u}\mathrm{l}}.$$(A.11)

The energy absorbed in frequency range dν, solid angle dΩ unit time dt by a surface element dA is $$dE=\frac{h{\nu }_0}{4\pi }{n}_{\text{l}}{B}_{\mathrm{l}\mathrm{u}}{\mathrm{\Phi }}_{\nu }{I}_{\nu }dAdsdtd\mathrm{\Omega }d\nu .$$(A.12)

We can now apply the relation between energy and intensity: dE = I_νdA dt dΩ dν. Using the above two relations, the absorption coefficient of the gas is $${\chi }_{\nu }(s)=\frac{h\nu }{4\pi }{n}_{\text{l}}{B}_{\mathrm{l}\mathrm{u}}{\mathrm{\Phi }}_{\nu }.$$(A.13)

We have seen that stimulated emission is proportional to the mean intensity. Furthermore, we assume that it only affects the photons along the given beam. In this case we can view stimulated emission as negative absorption. As a consequence, the absorption coefficient of the gas corrected for stimulated emission is $${\chi }_{\nu }(s)=\frac{h\nu }{4\pi }(n_l{B}_{lu}-n_u{B}_{ul}){\mathrm{\Phi }}_{\nu }.$$(A.14)

Gas is assumed to be a locally homogeneous medium, which means that ρ = nm_mol, where m_mol is the mass density of the ions. Both are constant along the propagation of the radiation. Combining all our previous assumptions, the mass absorption coefficient of the gas can be written as $${\kappa }_{\nu }=\frac{1}{n{m}_{\mathrm{m}\mathrm{o}\mathrm{l}}}\frac{h\nu }{4\pi }(n_{\text{l}}{B}_{\mathrm{l}\mathrm{u}}-n_{\mathrm{u}}{B}_{\mathrm{u}\mathrm{l}}){\mathrm{\Phi }}_{\nu }.$$(A.15)

Assuming only thermal excitation, the total number of the molecules, namely the number density, is the sum of molecules in the ith state, $$n={\displaystyle \sum _{\mathrm{i}}}{n}_{\mathrm{i}}=\frac{n_0}{g_0}{\displaystyle \sum _{\mathrm{i}}}{g}_{\mathrm{i}}\exp (-\frac{{\epsilon }_{\mathrm{i}}-{\epsilon }_0}{kT}).$$(A.16)

In the above expression, n₀ is the number of ions in the ground state, and n_i is the number of ions in the i th excited state, while g₀ and g_i are the degeneracies of the ground and excited states, respectively. If the energy of the ground state is ε₀, and that of the excited states are ε_i, we can introduce Q_T, the partition function, which can be approximated as follows, $$Q_{\mathrm{T}}={\displaystyle \sum _i}{g}_{\mathrm{i}}\exp (-\frac{{\epsilon }_{\mathrm{i}}-{\epsilon }_0}{kT}).$$(A.17)

The level population of the i^th excited state can be expressed by $$n_{\mathrm{i}}=n\frac{g_{\mathrm{i}}}{Q_{\mathrm{T}}}\exp (-\frac{{\epsilon }_{\mathrm{i}}-{\epsilon }_0}{kT}).$$(A.18)

It is expedient to find an approximate formula for Q_T as a function of temperature T in a polynomial form, given by $$\ln (Q_{\mathrm{T}}-g_0)={\displaystyle \sum _{\mathrm{i}}}{a}_{\mathrm{i}}\ln {\left(\frac{5040}{T}\right)}^{\mathrm{i}}.$$(A.19)

The coefficients a_i for the Ca II ion can be calculated by the transitional parameters given by Bolton (1970), and can be found in Table A.1.

Combining Eqs. (A.11), (A.15) and (A.18), the monochromatic mass absorption coefficient of the gas at temperature T for a given transition at frequency ν₀ can be expressed as $${\kappa }_{\nu }=\frac{1}{8\pi }\frac{1}{Q_{\mathrm{T}}}\frac{A_{\mathrm{u}\mathrm{l}}{g}_{\mathrm{u}}}{m_{\mathrm{m}\mathrm{o}\mathrm{l}}}\frac{c^2}{{\nu }_0^2}[\exp (-\frac{{\epsilon }_1}{kT})-\exp (-\frac{{\epsilon }_{\mathrm{u}}}{kT})]{\mathrm{\Phi }}_{\nu }.$$(A.20)

Regarding the calculation of the intrinsic line profile function Φ_ν, we have to take into consideration the Doppler line broadening caused by the thermal motion of emitting molecules. The frequency of the emitted radiation of the molecule in its own frame is ν₀, which corresponds to a frequency of ν₀ ± Δν for an observer. The line emitted by a certain amount of gas is broadened without change in the net flux. Assuming a non-relativistic line-of-sight velocity (v_los/c ≪ 1), the Doppler shift causes a change of ν − ν₀ = ν₀(v_los/c) in the observed frequency. The line-of-sight velocity of the gas in the observer’s frame, v_los, can be calculated from the velocity components obtained from the hydrodynamical simulations: v_los = [v_r sin(ϕ) + v_t cos(ϕ)] cos(i), or one can assume it to be Keplerian. The line-of-sight velocity is responsible for the width of the lines, and the position of the line centre. Assuming the Maxwell speed distribution for the ions, the number of molecules emitting in the frequency range of ν₀ + dν can be written as $$n(\nu )=\exp [-{\left(\frac{c}{{\nu }_0}\right)}^2\frac{m_{\mathrm{m}\mathrm{o}\mathrm{l}}}{2kT}{(\nu -{\nu }_0)}^2].$$(A.21)

Since the strength of the emission (i.e. Φ_ν) is proportional to the number of emitting ions, the thermally broadened line profile is described by a Gaussian function, $${\mathrm{\Phi }}_{\nu }=\frac{1}{\sigma \sqrt{\pi }}\exp [-{\left(\frac{\nu -{\nu }_0}{\sigma }\right)}^2].$$(A.22)

The broadening of the lines is caused both by thermal excitation and the turbulent movement of the ions. Thermal broadening is denoted by σ_term and, from Eq. (A.21) can be given as $${\sigma }_{\text{term\hspace{0.17em}}}=\frac{{\nu }_0}{c}\sqrt{\frac{2kT}{m_{\mathrm{m}\mathrm{o}\mathrm{l}}}}.$$(A.23)

The turbulent line broadening, σ_turb, is $${\sigma }_{\text{turb\hspace{0.17em}}}=\eta \sqrt{\frac{\mathrm{\Gamma }{k}_{\mathrm{B}}T}{m_{\mathrm{C}\mathrm{a}}}},$$(A.24)

where the adiabatic index of the gas is assumed to be Γ = 1.4, and m_Ca is the mass of the Ca atom. The turbulent velocity is assumed to be η times the sound speed in the gas. However, in this study we use the Shakura & Sunyaev α-parameter, so the turbulent broadening becomes $${\sigma }_{\text{turb\hspace{0.17em}}}=\alpha \sqrt{\frac{\mathrm{\Gamma }{m}_{\mathrm{m}\mathrm{o}\mathrm{l}}}{m_{\mathrm{C}\mathrm{a}}}}{\sigma }_{\text{term}}.$$(A.25)

The total line broadening is then $\sigma =\sqrt{{\sigma }_{\text{term\hspace{0.17em}}}^2+{\sigma }_{\text{turb\hspace{0.17em}}}^2}$.

Due to the gas orbiting the central star, the line profile is shifted by Δν compared to the ν₀ frequency. Taking into account the Doppler-shift, the line profile at a given point of the disc (r, ϕ) can be given as $${\mathrm{\Phi }}_{\nu }=\frac{1}{\sigma \sqrt{\pi }}\exp [-{\left(\frac{\nu -{\nu }_0+\mathrm{\Delta }\nu }{\sigma }\right)}^2].$$(A.26)

To calculate the final line profile one has to integrate the Φ_ν profiles in each cell of the hydrodynamic simulation grid. Taking into consideration the inclination of the disc, the final line profile is $$\mathrm{\Phi }=\int {\mathrm{\Phi }}_{\nu }dA\cos (i),$$(A.27)

where dA is the area of the emitting grid cell.

The data of all three lines of the Ca II triplet can be found in Table A.2 with the line relevant for this work highlighted in bold. In this study, the 8542 Å wavelength Ca II line is calculated.

Appendix B Effects of a perturbed source function for blackbody radiation on line morphology

As previously stated in Appendix A, the canonical form of the Planck function described by Eq. A.6 may not be able describe an optically thin disc in LTE. We undertook an analysis of two perturbed Planck functions. In the first scenario, it is subject to a power-law perturbation in temperature space: $$B_{\nu }(T)=\frac{2h{\nu }^3}{c^2}{[\exp \left(\frac{h\nu }{kT}\right)-1]}^{-1}{\left(\frac{T}{T_{\text{in\hspace{0.17em}}}}\right)}^{\beta },$$(B.1)

where T_in is the temperature measured at the inner edge of the disc (5000 K). In the second case, the function is perturbed in the frequency domain as follows, $$B_{\nu }(T)=\frac{2h{\nu }^3{\nu }^{\ast }}{c^2}{[\exp \left(\frac{h\nu {\nu }^{\ast }}{kT}\right)-1]}^{-1},$$(B.2) $${\nu }^{\ast }={\left(\frac{\nu }{{\nu }_0}\right)}^{\alpha },$$(B.3)

where ν₀ is the centre of the λ = 8542 Å line.

The top panel of Fig. B.1 illustrates B_ν as a function of frequency, with varying values of β. The bottom panel depicts the calculated λ = 8542 Å emission line for each case, assuming a Keplerian disc that is unperturbed in both density and velocity space. As can be observed, the morphology of the background-normalised lines does not exhibit significant variation from the unperturbed case. However, the absolute flux of the lines undergoes a slight alteration in response to the perturbation on a range from 7 × 10¹³ erg s⁻¹ cm⁻² Hz⁻¹ at β = 2 to 2.25 × 10¹⁴ erg s⁻¹ cm⁻² Hz⁻¹ at β = −2.

Figure B.2 depicts B_ν as a function of frequency, with varying values of α. The lines perturbed in frequency space are not shown, as they are practically indistinguishable from the unperturbed line profiles, both in terms of morphology and absolute flux. Furthermore, selecting a value of ν₀ that is either half or twice the initial value (the value corresponding to λ = 8542 Å) does not result in a notable change in the background-normalised line morphology.

Based on our analysis, neither the temperature nor the frequency space perturbations cause a difference in the relative strength of the Ca II triplet lines. This is due to the fact that the applied perturbation on the Planck function results in a very small variation in the close vicinity of the Ca II triplet. As a summary, we conclude that the applied perturbations have little to no effect on the line morphologies, and we proceed with assuming an canonical source function for blackbody radiation.

Fig. B.1 Effect of temperature domain perturbations on the Ca II line profile. Top: Planck functions as a function of frequency. Different colours correspond to different perturbations in the temperature domain (see Eq. B.1). The unperturbed model is denoted by a black line. The inlet shows the close proximity of the Ca II triplet, marked by a dashed rectangle. Bottom: background-normalised λ = 8542 Å Ca II line profiles calculated using the perturbed Bν functions.

Fig. B.2 Same as the top panel of Fig. B.1, but with the perturbation occurring in frequency space (see Eq. B.2). The unperturbed model is shown by the solid grey line.

Appendix C Spectrograms

Figure C.1 displays spectrograms of all models, showing a series of Ca II line profiles as a function of time in velocity and wavelength space. The decrease of line intensity can clearly be observed (see, e.g. model B3), which corresponds to the decrease in flux and disc mass and the eventual depletion of the disc.

Fig. C.1 Spectrograms of all models. Colours code the absolute flux of the Ca II lines.

References

All Tables

All Figures

	Fig. 1 Summary of the parameters describing line morphology: asymmetry ξ, central dip, shift of the line centre in velocity space (vshift) and half width at base (HWB).
In the text

	Fig. 2 Process of disc formation. When the asteroid passes under the Roche limit, it loses material, which initially exists as dust until the asteroid passes the sublimation radius (T = 1500 K), within which gas is produced. The gas then spreads out viscously and forms a disc. The inner limit of Ca II emission is T = 5000 K. We note that the figure is not to scale.
In the text

	Fig. 3 Arrangement of asteroid orbits (marked by red lines) at different eccentricities relative to the sublimation radius (white dashed lines). We note that eccentric asteroids spend only part of their orbits below the sublimation radius.
In the text

	Fig. 4 Phases of disc formation: snapshots of the density distribution in model B1. Grey dashed lines mark the sublimation radius. A black arrow shows the orbital direction of the asteroid.
In the text

Fig. 5 Density and velocity distributions (Σ and vlos, respectively) in model B 1 at three points in time (top row) and the corresponding background-normalised Ca II lines (bottom row). The discs are viewed from below. The sublimation radius is marked by a black dashed line, and the orbital direction of the planet is shown with a black arrow. In the Panels showing the Ca II lines, vertical dashed lines indicate the centre of the line. The retrograde precession of the disc causes a periodic shift in the intensity maximum.

In the text

Fig. 6 Comparison of a Keplerian disc and model B1. The upper panel displays the line-of-sight velocity distribution in a Keplerian disc as well as that of model B1 from the hydrodynamic simulations (t = 132.3 days). The sublimation radius is indicated by black dashed lines. The lower panel shows the calculated Ca II line profiles with dashed and solid lines in the two cases. The line centres are marked by vertical grey lines.

In the text

	Fig. 7 Ca II line profiles for model B1 plotted with different disc inclinations. Fainter lines represent larger inclinations.
In the text

	Fig. 8 Ca II line profiles assuming different stellar temperatures in model B1. The snapshots depict the largest possible red-side asymmetry.
In the text

Fig. 9 Spatial extent of three gas and dust discs that are detected by Melis et al. (2010). The dust discs are marked by grey rectangles, while the gas discs are marked by orange rectangles. To improve visibility, the discs for each system are offset vertically. For the two stellar effective temperatures used in the Melis et al. (2010) study, namely 18 500 K (red lines) and 22 000 K (blue lines), the blackbody temperature profile (light red and blue lines) is shown. A horizontal green line indicates the sublimation temperature (1500 K). The figure displays the sublimation radius for each configuration with vertical dashed lines. The dashed grey line indicates the radius of the WD, while the solid black line represents the Roche limit.

In the text

Fig. B.1 Effect of temperature domain perturbations on the Ca II line profile. Top: Planck functions as a function of frequency. Different colours correspond to different perturbations in the temperature domain (see Eq. B.1). The unperturbed model is denoted by a black line. The inlet shows the close proximity of the Ca II triplet, marked by a dashed rectangle. Bottom: background-normalised λ = 8542 Å Ca II line profiles calculated using the perturbed Bν functions.

In the text

	Fig. B.2 Same as the top panel of Fig. B.1, but with the perturbation occurring in frequency space (see Eq. B.2). The unperturbed model is shown by the solid grey line.
In the text

	Fig. C.1 Spectrograms of all models. Colours code the absolute flux of the Ca II lines.
In the text