Polluted White Dwarfs: Mixing Regions and Diffusion Timescales

TL;DR

This paper models mixing processes in polluted white dwarf atmospheres to better understand accretion rates and the influence of various mixing mechanisms, providing new diffusion timescales and implications for planetary debris accretion.

Contribution

It introduces detailed MESA white dwarf models that explore convection, sedimentation, overshoot, and thermohaline mixing, offering new diffusion timescales and accretion rate estimates.

Findings

Accretion rates up to 10^13 g/s explain observed pollution in hot DA WDs.

Thermohaline mixing dominates sedimentation in certain regimes.

Polluted DB WDs experience less thermohaline mixing, especially below 18,000 K.

Abstract

Many isolated white dwarfs (WDs) show spectral evidence of atmospheric metal pollution. Since heavy element sedimentation timescales are short, this most likely indicates ongoing accretion. Accreted metals encounter a variety of mixing processes at the WD surface: convection, gravitational sedimentation, overshoot, and thermohaline instability. We present MESA WD models that explore each of these processes and their implications for inferred accretion rates. We provide diffusion timescales for many individual metals, and we quantify the regimes in which thermohaline mixing dominates over gravitational sedimentation in setting the effective settling rate of the heavy elements. We build upon and confirm earlier work finding that accretion rates as high as $10^{13} \, \rm g \, s^{-1}$ are needed to explain observed pollution in DA WDs for $T_{\rm eff} > 15,000 \, \rm K$, and we provide tabulated results from our models that enable accretion rate inferences from observations of polluted DA WDs. If these rates are representative of young WDs, we estimate that the total mass of planetesimal material accreted over a WD lifetime may be as high as $10^{28} \, \rm g$, though this estimate is susceptible to potential selection biases and uncertainties about the nature of disk processes that supply accretion to the WD surface. We also find that polluted DB WDs experience much less thermohaline mixing than DA WDs, and we do not expect thermohaline instability to be active for polluted DB WDs with $T_{\rm eff} < 18,000 \, \rm K$.

Tables

Table 1. Table 1: Comparison of MESA  and Koester ( 2009 ) results for the mass of the surface convection zone and diffusion timescales for Ca 40 superscript Ca 40 {}^{40}{\rm Ca} on a 0.6 ​ M ⊙ 0.6 subscript 𝑀 direct-product 0.6\,M_{\odot} WD. The models of Koester ( 2009 ) adopt a fixed value of log ⁡ g = 8.0 𝑔 8.0 \log g=8.0 , so we include a column for log ⁡ g 𝑔 \log g from the MESA  model to note how it varies slightly about this value as the WD cools. Columns labeled with D Paq subscript 𝐷 Paq D_{\rm Paq} refer to diffusion calculations using the coefficients of Paquette et al. ( 1986a ) , and those labeled D SM subscript 𝐷 SM D_{\rm SM} refer to calculations using the coefficients of Stanton & Murillo ( 2016 ) . Columns labeled with Z 2 + 1 subscript 𝑍 2 1 Z_{2}+1 or Z 2 − 1 subscript 𝑍 2 1 Z_{2}-1 refer to calculations for which the charge of Ca 40 superscript Ca 40 {}^{40}{\rm Ca} is taken to be 1 larger or smaller than the value given by the ionization calculations based on Paquette et al. ( 1986b ) .

$T_{\mathrm{eff}}$ [K]	$\log \left(M_{\mathrm{cvz}}/M\right)$		$\log g$	$\log \left({\tau }_{\mathrm{diff}}/\mathrm{yr}\right)$ for ${}_{}{}^{40}\mathrm{Ca}$
	Koester	MESA	MESA	Koester	MESA		MESA ($Z_2+1$)		MESA ($Z_2-1$)
				($D_{\mathrm{Paq}}$)	($D_{\mathrm{Paq}}$)	($D_{\mathrm{SM}}$)	($D_{\mathrm{Paq}}$)	($D_{\mathrm{SM}}$)	($D_{\mathrm{Paq}}$)	($D_{\mathrm{SM}}$)
6000	-7.722	-7.8094	8.0342	4.2924	4.2449	4.13	4.2977	4.1827	4.1835	4.0684
7000	-8.432	-8.5222	8.0306	3.7125	3.7107	3.6139	3.7662	3.6695	3.6473	3.5513
8000	-8.954	-8.9849	8.0272	3.3303	3.4113	3.3372	3.4674	3.3923	3.3616	3.2908
9000	-9.607	-9.517	8.0238	2.8725	3.0408	2.9957	3.1009	3.0531	2.9715	2.9304
10000	-10.738	-10.251	8.0202	1.9997	2.476	2.4679	2.5493	2.5371	2.3884	2.3852
11000	-12.715	-11.872	8.0164	0.4845	1.1984	1.2236	1.3214	1.3478	1.0337	1.0566
12000	-15.618	-14.698	8.0127	-1.6941	-1.0767	-1.071	-0.81118	-0.80264	-1.5571	-1.5573
13000	-16.408	-16.103	8.0094	-2.359	-1.9677	-1.9629	-1.7151	-1.7073	-2.4523	-2.4534
14000	-16.672	-16.292	8.006	-2.6305	-2.0968	-2.0931	-1.8252	-1.8185	-2.5713	-2.5734
15000	-16.698	-16.43	8.0026	-2.6277	-2.1953	-2.1926	-1.9216	-1.9159	-2.6857	-2.6887
16000	-16.744	-16.622	7.9991	-2.622	-2.3333	-2.3318	-2.0573	-2.0526	-2.8272	-2.8312
17000	-16.634	-16.836	7.9953	-2.4688	-2.4941	-2.4939	-2.2153	-2.2119	-2.9901	-2.9952
18000	-16.586	-16.787	7.9914	-2.4213	-2.469	-2.4694	-2.1889	-2.1862	-2.9668	-2.9724
19000	-16.538	-16.703	7.9872	-2.3804	-2.4182	-2.4191	-2.1374	-2.1352	-2.9171	-2.9231
20000	-16.439	-16.644	7.983	-2.3077	-2.3847	-2.3862	-2.1033	-2.1015	-2.8849	-2.8913

Table 2. Table 2: MESA  diffusion timescales for the 0.6 ​ M ⊙ 0.6 subscript 𝑀 direct-product 0.6\,M_{\odot} WD model calculated using the coefficients of Stanton & Murillo ( 2016 ) . Supplemental tables for other WD masses are available at https://doi.org/10.5281/zenodo.2541235 (Bauer, 2019 ) .

$T_{\mathrm{eff}}$ [K]	$\log \left(M_{\mathrm{cvz}}/M\right)$	$\log \left({\tau }_{\mathrm{diff}}/\mathrm{yr}\right)$
		${}_{}{}^{12}C$	${}_{}{}^{16}O$	${}_{}{}^{23}\mathrm{Na}$	${}_{}{}^{24}\mathrm{Mg}$	${}_{}{}^{27}\mathrm{Al}$	${}_{}{}^{28}\mathrm{Si}$	${}_{}{}^{40}\mathrm{Ca}$	${}_{}{}^{48}\mathrm{Ti}$	${}_{}{}^{52}\mathrm{Cr}$	${}_{}{}^{56}\mathrm{Fe}$
6000	-7.8094	4.2573	4.3303	4.0949	4.0564	4.001	3.8702	4.13	4.0304	3.931	3.89
7000	-8.5222	3.8338	3.9012	3.646	3.6249	3.4721	3.4464	3.6139	3.5067	3.4046	3.4617
8000	-8.9849	3.5523	3.6214	3.3551	3.3105	3.1875	3.1688	3.3372	3.1957	3.0915	3.185
9000	-9.517	3.2304	3.2723	2.9774	2.935	2.8664	2.8485	2.9957	2.8552	2.7713	2.8289
10000	-10.251	2.7752	2.7551	2.4842	2.4591	2.2849	2.3975	2.4679	2.3434	2.3276	2.3132
10500	-10.928	2.3481	2.2194	1.9873	1.8792	1.8257	1.973	1.9514	1.8086	1.9119	1.8567
11000	-11.872	1.7424	1.552	1.2795	1.0427	1.2254	1.3768	1.2236	1.1432	1.233	1.184
11500	-13.147	0.74402	0.57918	0.21533	0.17977	0.38994	0.54935	0.22032	0.31781	0.28456	0.21358
12000	-14.698	-0.55196	-0.70814	-1.3172	-0.85012	-0.63472	-0.65315	-1.071	-0.88082	-0.9153	-0.94716
12500	-15.953	-1.6971	-1.9663	-2.124	-1.6439	-1.6901	-1.7108	-1.8661	-1.9454	-1.9789	-2.0115
13000	-16.103	-1.9121	-2.0535	-2.2121	-1.74	-1.8271	-1.8073	-1.9629	-2.0423	-2.0772	-2.1094
13500	-16.203	-1.9954	-2.1254	-2.2838	-1.809	-1.9515	-1.8764	-2.0319	-2.1114	-2.1462	-2.1785
14000	-16.292	-2.0585	-2.1858	-2.3321	-1.87	-2.0932	-1.9374	-2.0931	-2.1725	-2.2074	-2.2272
15000	-16.43	-2.1621	-2.2885	-2.4473	-1.9693	-2.2683	-2.0368	-2.1926	-2.2721	-2.3071	-2.3393
16000	-16.622	-2.3042	-2.4307	-2.5896	-2.1081	-2.6023	-2.1757	-2.3318	-2.4114	-2.4464	-2.4786
17000	-16.836	-2.4658	-2.5932	-2.7529	-2.2694	-2.8193	-2.3373	-2.4939	-2.5737	-2.6088	-2.6411
18000	-16.787	-2.4429	-2.5704	-2.7301	-2.2449	-2.6484	-2.3128	-2.4694	-2.5493	-2.5843	-2.6167
19000	-16.703	-2.3906	-2.5212	-2.6809	-2.1946	-2.2487	-2.2625	-2.4191	-2.4989	-2.5339	-2.5663
20000	-16.644	-2.261	-2.4894	-2.6491	-2.1617	-2.2136	-2.2296	-2.3862	-2.466	-2.5009	-2.5333

Table 3. Table 3: Critical accretion rates in our 0.6 ​ M ⊙ 0.6 subscript 𝑀 direct-product 0.6\,M_{\odot} MESA  WD models above which thermohaline mixing modifies the equilibrium surface mass fractions from the prediction of Equation ( 2 ). This assumes accretion of material with bulk earth composition (McDonough, 2001 ) .

$T_{\mathrm{eff}}$ [K]	6,000	7,000	8,000	9,000
${\dot{M}}_{\mathrm{crit}}$ [$g{s}^{-1}$]	${10}^{12}$	${10}^{10}$	${10}^9$	${10}^8$
$T_{\mathrm{eff}}$ [K]	10,000	11,000	$>12$,000
${\dot{M}}_{\mathrm{crit}}$ [$g{s}^{-1}$]	${10}^7$	${10}^6$

Equations

$${\tau }_{\mathrm{diff},i}\equiv \frac{M_{\mathrm{cvz}}}{4\pi r^2\rho v_{\mathrm{diff},i}},$$

$$X_{\mathrm{eq},i}=\frac{{\dot{M}}_i}{M_{\mathrm{cvz}}}{\tau }_{\mathrm{diff},i}.$$

$$\displaystyle v_{\rm diff}=D_{12}\bigg{[}$$

$$\displaystyle+\left(\frac{Z_{2}}{Z_{1}}-1\right)\frac{\partial\ln p_{i}}{\partial r}+\alpha_{T}\frac{\partial\ln T}{\partial r}\bigg{]}~{},$$

$${\nabla }_T-{\nabla }_{\mathrm{ad}}<\frac{\phi }{\delta }{\nabla }_{\mu }<0,$$

$$\frac{(\phi /\delta ){\nabla }_{\mu }}{{\nabla }_T-{\nabla }_{\mathrm{ad}}}>\frac{{\kappa }_{\mu }}{{\kappa }_T}.$$

$${\kappa }_T=\frac{4ac{T}^3}{3\kappa {\rho }^2{C}_P},$$

$$D_{\mathrm{th}}={\alpha }_{\mathrm{th}}{\kappa }_T\frac{3}{2}\frac{(\phi /\delta ){\nabla }_{\mu }}{{\nabla }_T-{\nabla }_{\mathrm{ad}}},$$

$$ϵ\equiv 1-\frac{{\kappa }_{\mu }/{\kappa }_T}{(\phi /\delta ){\nabla }_{\mu }/({\nabla }_T-{\nabla }_{\mathrm{ad}})},$$

$$\mathrm{Ro}\sim \sqrt{\frac{N^2}{4{\Omega }^2}\frac{(\phi /\delta ){\nabla }_{\mu }}{{\nabla }_T-{\nabla }_{\mathrm{ad}}}},$$

$$\mathrm{Ro}\sim 10^{-2}\sqrt{\frac{R}{H}}\left(\frac{{\Omega }_{\mathrm{crit}}}{\Omega }\right).$$

$$D_{\mathrm{overshoot}}(r)=D_0\exp \left(-\frac{2\left|r-r_0\right|}{H_P}\right),$$

$${\tau }_{\mathrm{ovr},i}\equiv \frac{M_{\mathrm{ovr}}}{4\pi r^2\rho v_{\mathrm{ovr},i}},$$

$$X_{\mathrm{eq},i}=\frac{{\dot{M}}_i}{M_{\mathrm{ovr}}}{\tau }_{\mathrm{ovr},i}$$

$${\Sigma }_{st}^{(lj)}=\frac{4\pi }{{\pi }^{3/2}}{\int }_0^{\infty }dv\exp \left(\frac{-v^2}{{\alpha }_{st}^2}\right)\frac{v^{2j+3}}{{\alpha }_{st}^{2j+4}}{S}_{st}^{(l)},$$

$$S_{st}^{(l)}=2\pi {\int }_0^{\infty }(1-{\cos }^l{\chi }_{st})bdb$$

$$S_{st}^{(1)}=2.41\pi Ze{\left(\frac{{\alpha }_N}{{\mu }_{st}}\right)}^{1/2}\frac{1}{v},$$

$${\Sigma }_{st}^{(11)}=3.62\pi \frac{Ze}{{\alpha }_{st}}{\left(\frac{{\alpha }_N}{{\mu }_{st}}\right)}^{1/2}.$$

$$K_{st}=\frac{2}{3}{n}_s{n}_t{\mu }_{st}{\alpha }_{st}{\Sigma }_{st}^{(11)}=2.41\pi n_s{n}_{t}Ze{\left({\mu }_{st}{\alpha }_N\right)}^{1/2}.$$

$$K_{st}=\frac{16\sqrt{\pi }}{3}\frac{n_s{n}_t{Z}_s^2{Z}_t^2{e}^4}{{\mu }_{st}{\alpha }_{st}^3}\ln \Lambda ,$$

$$\frac{K_{st}^{(\mathrm{dipole})}}{K_{st}^{(\mathrm{coulomb})}}=\frac{3.62(2k_{B}T{)}^{3/2}{\left(\pi {\alpha }_N\right)}^{1/2}}{8e^3\ln \Lambda }.$$

$$\frac{K_{st}^{(\mathrm{dipole})}}{K_{st}^{(\mathrm{coulomb})}}=\frac{0.054}{\ln \Lambda }.$$

Full text

accretion , accretion disks – diffusion – instabilities – minor planets, asteroids: general – planetary systems – white dwarfs

††software: MESA (Paxton et al., 2011, 2013, 2015, 2018), Astropy (Astropy Collaboration et al., 2013, 2018), Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011), SciPy (Jones et al., 2001–), MesaScript (Wolf et al., 2017)

1 Introduction

A large fraction ($25\%$–$50\%$) of isolated white dwarf (WD) atmospheres show signatures of polluting metals (Zuckerman et al., 2003; Koester et al., 2014). Heavy element sedimentation timescales are short (Schatzman, 1945, 1948), and this implies recent or ongoing accretion of observed heavy elements (Vauclair et al., 1979; Koester et al., 2014). Polluted WD spectra are often accompanied by infrared emission from a dust disk (Koester et al., 1997; Farihi et al., 2009; Girven et al., 2012; Farihi, 2016), and the predominant model for the origin of this dust is debris from disrupted planetesimals (Jura, 2003; Jura & Young, 2014; Vanderburg et al., 2015). Models for WD surface mixing allow inferences of the composition of these planetesimals and the rates at which WDs accrete this material (Koester & Wilken, 2006; Koester, 2009; Dufour et al., 2010, 2012; Koester et al., 2011; Farihi et al., 2013; Raddi et al., 2015).

While many have relied on elemental sedimentation timescales to make inferences about polluted WD accretion, recent work has revealed that thermohaline instability is active and substantially modifies the inferred accretion rate $\dot{M}_{\rm acc}$ (Deal et al., 2013; Wachlin et al., 2017). Our work in Bauer & Bildsten (2018, Paper I) greatly expanded the range of $T_{\rm eff}$ explored for polluted WD models accounting for thermohaline mixing. In Paper I, we found that some DA WDs must accrete bulk earth composition at rates as large as $\dot{M}_{\rm acc}=10^{13}\,\rm g\,s^{-1}$ for our models to match observed surface metal abundances.

In this work, we build on the results of Paper I with further examination of the surface mixing processes relevant for heavy element pollution. We construct models that include these processes using the stellar evolution code MESA (Paxton et al., 2011, 2013, 2015, 2018). We use MESA version 10398 unless otherwise stated. For models that include thermohaline mixing, we use MESA version 11191 (see Section 4.1 for a description of the relevant changes to the code that motivate using this version). MESA inlists and other input files necessary to reproduce our MESA models are available online at https://doi.org/10.5281/zenodo.2541235. In Section 2, we describe the hydrogen-dominated surface convection zones that metals first encounter when they accrete onto DA WDs. In Section 3, we quantify the individual element diffusion timescales for sedimentation beneath the convection zones in our MESA models, and provide tabulated results. In Section 4, we examine other forms of mixing that can be relevant in non-convective regions. These include a greatly refined and expanded treatment of thermohaline mixing (Section 4.1), as well as convective overshoot (Section 4.2). Our results confirm the findings of Paper I that earth composition accretion rates approaching $\dot{M}_{\rm acc}=10^{13}\,\rm g\,s^{-1}$ are necessary to match observed calcium abundances in DA WDs with $T_{\rm eff}\gtrsim 15,000\,\rm K$. We close with discussion and conclusions in Sections 5 and 6.

2 Surface Convection Zones in Pure Hydrogen

Polluting metals quickly mix into the WD surface convection zone, so that the base of this fully mixed region is where the gravitational sedimentation rate must be evaluated (Vauclair et al., 1979; Koester, 2009). Here we quantify the total mass, $M_{\rm cvz}$, in the surface convective layers of our MESA WD models with pure hydrogen atmospheres and compare to previous work.

To facilitate comparison to the work of Koester (2009, 2010), we adopt the ML2 convection prescription (Bohm & Cassinelli, 1971) with $\alpha_{\rm MLT}=0.8$. This value of $\alpha_{\rm MLT}$ is similar to those calibrated against 3D convection by Tremblay et al. (2013, 2015), but it should be noted that the calibrated values show some variation with $T_{\rm eff}$. The depth of the surface convection zone is also sensitive to the surface boundary condition of the model. We find that the gray iterative atmosphere procedure in MESA (Paxton et al., 2013) provides values of $M_{\rm cvz}$ in good agreement with Koester (2009) for $T_{\rm eff}\gtrsim 9,000\,\rm K$. At lower effective temperatures, we switch to the pre-computed DA WD atmosphere tables adapted from Rohrmann et al. (2012). When using these tables, MESA models attach to boundary conditions given at the optical depth $\tau_{\rm Ross}=25$, so the tables are not suitable for WDs with very shallow convection zones that do not extend to $\tau_{\rm Ross}>25$. Fortunately, the gray iterative atmosphere procedure is adequate for all cases of shallow convection zones, and it is only necessary to switch to the tables for cooler WDs with large convection zones. For the remainder of this work, we use models that switch from gray atmosphere boundary conditions to the WD atmosphere tables below $T_{\rm eff}=9,000\,\rm K$.

Figure 1 shows a comparison of the mass of the surface convection zone between MESA WD models and the DA models of Koester (2009, 2010).111 Most recent tables found at http://www1.astrophysik.uni-kiel.de/~koester/astrophysics/astrophysics.html. For hotter WDs where no surface convection is present ($T_{\rm eff}\gtrsim 15,000\,\rm K$), the mass exterior to the photosphere is the relevant parameter for pollution mixing calculations, so we show this value as well. The tables of Koester (2009) give whichever is larger of mass in the surface convection zone and mass exterior to the photosphere. The hydrogen ionization transition that drives convection results in a steep change in the convection zone mass around $10,000\,{\rm K}\lesssim T_{\rm eff}\lesssim 13,000\,\rm K$. We see small disagreements in the exact location of this feature, and otherwise are in excellent agreement with Koester (2009). Although the disagreement in mass at fixed $T_{\rm eff}$ can be up to an order of magnitude, the steep slope of the curve in this region means that small variations of $T_{\rm eff}$ within typical observational uncertainties can bring these values into agreement.

3 Diffusion Timescales

With the structure of MESA WD models and convection zone masses established, diffusion timescales can now be calculated for all trace heavy elements. These timescales are essential to inferring accretion rates and compositions from observations of trace elements in the photosphere. Section 2 of Paper I shows the equations relating these diffusion timescales to accretion rates and observable surface abundances. Here we only repeat a few key definitions for convenience. When no other mixing occurs beneath the surface convection zone, the sedimentation timescale for trace element $i$ is

[TABLE]

where $r$ is the radius, $\rho$ is the density, and $v_{{\rm diff},i}$ is the sedimentation velocity of element $i$ evaluated at the base of the surface convection zone where it sinks away from the fully mixed surface region. An approximate expression for $v_{{\rm diff},i}$ for a trace diffusing element is given later in Equation (3), but in general our MESA models calculate diffusion velocities from a full solution of the Burgers equations (Burgers, 1969) as described in Paxton et al. (2015, 2018). For a constant accretion rate $\dot{M}_{i}$ of species $i$ over timescales much longer than $\tau_{{\rm diff},i}$, the surface mass fraction approaches the equilibrium value

[TABLE]

For observational inferences, it is assumed that this equilibrium state has been reached, so that the elemental accretion rate can be derived from the observed mass fraction as $\dot{M}_{i}=X_{{\rm obs},i}M_{\rm cvz}/\tau_{{\rm diff},i}$. We denote the total accretion rate as $\dot{M}_{\rm acc}\equiv\sum_{i}\dot{M}_{i}$.

3.1 MESA Diffusion Results

The diffusion velocities necessary to calculate diffusion timescales using Equation (1) are readily available from MESA models. We obtain these by simply introducing a polluting metal (e.g. ${}^{40}{\rm Ca}$) accreting at a rate of $\dot{M}_{\rm acc}=10^{7}\,\rm g\,s^{-1}$. After accretion takes place for many diffusion timescales, so that the abundance in the surface convection zone has reached equilibrium, we calculate the diffusion timescale using Equation (1) along with the diffusion velocity reported by MESA from the solution of the Burgers equations. These diffusion calculations include thermal diffusion and properly account for any degree of electron degeneracy as described in Paxton et al. (2018).

Diffusion calculations according to the Burgers equations rely on coefficients calculated using a binary scattering formalism. The well-established coefficients of Paquette et al. (1986a) are based on a screened Coulomb potential treatment for calculating binary Coulomb collision cross sections. Recent updates to MESA (Paxton et al., 2018) have included options for using the coefficients of Stanton & Murillo (2016), who provide an improvement upon this method with a more sophisticated treatment of screening. Table 1 shows some comparisons for diffusion timescales in a $0.6\,M_{\odot}$ WD including calculations using the coefficients of Paquette et al. (1986a). In general, both sets of coefficients give similar results except for the deepest convection zones, where the increased Coulomb screening due to electrons in the calculations of Stanton & Murillo (2016) allows for faster diffusion.

Table 1 also shows comparisons to the diffusion timescales given by Koester (2009, see link in Footnote 1 for the most up-to-date diffusion timescale results), which employ the coefficients of Paquette et al. (1986a). When using these same coefficients, the MESA timescale results agree well as long as the convection zone depth is comparable. For $T_{\rm eff}\gtrsim 11,000\,\rm K$, the convection zone depths differ by up to an order of magnitude between MESA and Koester (2009), and the diffusion timescales disagree accordingly. Table 2 gives MESA diffusion timescales for ten commonly observed elements, using the coefficients of Stanton & Murillo (2016).

3.2 Approaching Equilibrium

Figure 2 shows surface mass fractions for several accreting elements in a MESA model, first approaching equilibrium after accretion turns on and continues for many diffusion timescales, then sinking away after accretion shuts off. For comparison, this figure also shows the analytic solution described in Paper I  for this constant accretion rate for ${}^{16}\rm O$ with a diffusion timescale of $\log(\tau_{\rm diff}/\rm yr)=3.62$. This verifies that the metals approach the equilibrium surface mass fraction predicted by Equation (2) for the diffusion timescales given in Table 2. The accretion episode shown in Figure 2 has all elements accreting at equal rates ($10^{7}\,\rm g\,s^{-1}$ for each element) to illustrate the effects of the hierarchy of diffusion timescales. This manifests as a clear ordering of abundances, where those with the longest diffusion timescales appear as the most abundant over all phases. In contrast, Figure 3 shows a more realistic accretion scenario, where the elements accrete at the total rate $\dot{M}_{\rm acc}=10^{8}\,\rm g\,s^{-1}$, but with the bulk earth abundance ratios of McDonough (2001). In this case, both the relative accreted mass fraction and diffusion timescale for each element play a role in establishing the final hierarchy of observed surface abundances. Neither of these calculations include thermohaline mixing (see Section 4.1).

3.3 Ionization States for Trace Metals

The partial ionization of metals in the surface layers relevant for pollution has important effects for the diffusion timescales. If we denote the background material in which diffusion takes place by the index $1$ (hydrogen in the case of a DA WD atmosphere), and denote the pollutant by index $2$, then in the limit of a trace pollutant ($n_{2}\ll n_{1}$) its diffusion velocity can be expressed as (cf. Pelletier et al. 1986; Dupuis et al. 1992; Koester et al. 2014)

[TABLE]

where $c_{2}\equiv n_{2}/(n_{1}+n_{2})$ is the concentration of the pollutant, $p_{i}\equiv p_{1}+p_{2}$ is the ion pressure, $\alpha_{T}$ is the thermal diffusion coefficient, and $A$ and $Z$ refer to the mass and charge of each species respectively. Note that this equation is appropriate for any degree of electron degeneracy (Pelletier et al., 1986), and it agrees with our MESA diffusion treatment based on the Burgers equations in the limit of trace particles diffusing in a hydrogen background. The charge of each species influences the diffusion velocity in two important ways: the direct influence on the forcing terms felt by each ion seen in Equation (3), and the influence of the charge of each particle on the Coulomb scattering that results in the diffusion coefficient $D_{12}$. The diffusion coefficient is related to the resistance coefficients used for MESA diffusion calculations described in Paxton et al. (2015, 2018) by $D_{ij}\propto K_{ij}^{-1}$. For Coulomb collisions, the resistance coefficients described in Paxton et al. (2015) scale with the charge approximately as $K_{ij}\propto(Z_{i}Z_{j})^{2}$, and hence diffusion calculations can be very sensitive to the ionization treatment adopted for the partially ionized surface regions of WD models.

Formally, each ionization state of a given element may be treated as a separate species with its own integer charge $Z_{i}$ for purposes of diffusion calculations. In order to simplify the problem, MESA calculations instead adopt an average state $\bar{Z}$ for each element as described in Paxton et al. (2015, 2018) so that each isotope corresponds to only one diffusion species. We use the ionization treatment of Paquette et al. (1986b) to find an average charge state for each diffusion species everywhere in the MESA model.222We note that the expression in Paquette et al. (1986b) for the depression of the continuum for ionization potentials contains a typo in Equation (21), where a factor of $\rho^{1/3}$ is missing from the last line. The MESA ionization routine instead follows Equation (3) of Dupuis et al. (1992), which correctly includes this factor. Our ionization treatment is very similar to that of Koester (2009), who also notes correcting the missing factor of $\rho^{1/3}$ for the most recent calculations hosted on his website (see link in Footnote 1).

Figure 4 displays some of the charges used as input for diffusion calculations reported in MESA WD models. Since the ionization procedure based on Paquette et al. (1986b) involves comparing ionization potentials to an effective threshold potential, it always selects an integer value for the average charge. This results in the stair-stepped profiles seen in Figure 4, which have been smoothed slightly to improve the numerical stability of diffusion calculations. The last columns of Table 1 present results from diffusion calculations for which the charge $Z_{2}$ is taken to be one larger or smaller than the value obtained from the Paquette et al. (1986b) routine. Comparison of these timescales quantifies the rough uncertainty associated with the average ionization calculations here.

Our diffusion calculations assume that every species is at least singly ionized. Diffusion coefficients for a neutral species require collision integrals for dipole scattering, which result in significantly smaller collision cross sections and correspondingly faster diffusion timescales (Appendix A). Options for such diffusion coefficients are not currently available in MESA. Since diffusion fluxes for neutral elements can be much faster than those for singly ionized elements, even a small fraction of neutral particles in the relevant layer can significantly modify overall sedimentation timescales, and it is no longer appropriate to treat ionization with an average charge $\bar{Z}<1$. Diffusion timescales presented in this work are only accurate for models where surface temperatures are hot enough or surface convection zones reach depths sufficient for at least single ionization of pollutants. Due to a thin or absent surface convection zone, these conditions fail to be satisfied around $T_{\rm eff}\approx 15,000\,\rm K$, and corresponding disagreement is evident between our results and those of Koester (2009) in Table 1 for this regime.

4 Other Mixing

We now explore additional mixing other than element diffusion beneath the convective layer. We focus on two fluid processes that can cause additional mixing: the thermohaline instability and convective overshoot.

4.1 Thermohaline Mixing

In the context of WD pollution, Deal et al. (2013) were the first to explore the possibility that accreted metals in WD atmospheres may lead to thermohaline instability. Subsequent work by Wachlin et al. (2017) confirmed the importance of the resulting mixing. In Paper I, we extended the parameter space for polluted WDs where thermohaline instability may occur, finding that thermohaline mixing significantly modifies inferred accretion rates in hydrogen-atmosphere WDs with $T_{\rm eff}\gtrsim 10,000\,\rm K$, with some rates reaching ${\dot{M}_{\rm acc}\approx 10^{13}\,\rm g\,s^{-1}}$ for $T_{\rm eff}\gtrsim 15,000\,\rm K$. Our exploration in Paper I was limited to WD models of mass $M=0.6\,M_{\odot}$ ($\log g=8.0$). We now expand upon that work with models of other masses to allow interpolation in $\log g$. We also adopt a refined treatment of thermohaline mixing based upon the work of Brown et al. (2013), which is calibrated against 3D simulations.

Two criteria must be satisfied for the thermohaline instability to be active. First, there must be an inverted molecular weight gradient in a region that is stable to convection:

[TABLE]

where $\nabla_{T}=(\partial\ln T/\partial\ln P)$ is the temperature gradient in the fluid, $\nabla_{\rm ad}=(\partial\ln T/\partial\ln P)_{s}$ is the adiabatic temperature gradient, ${\nabla_{\mu}=(\partial\ln\mu/\partial\ln P)}$ is the mean molecular weight gradient, ${\varphi=(\partial\ln\rho/\partial\ln\mu)_{P,T}}$, and ${\delta=-(\partial\ln\rho/\partial\ln T)_{P,\mu}}$. The instability is then driven by thermal exchange of perturbed fluid elements with their surroundings, whereupon a density contrast due to $\nabla_{\mu}$ leads to further mixing. Thus, the second criterion for thermohaline instability is that the magnitude of the molecular weight gradient $\nabla_{\mu}$ and the thermal diffusivity $\kappa_{T}$ must be large enough to excite the instability before particle diffusivity $\kappa_{\mu}$ within a perturbed fluid element can adjust its composition (Baines & Gill, 1969; Garaud, 2018):

[TABLE]

Assuming that heat transport is radiative, the thermal diffusivity is

[TABLE]

where $\kappa$ is the opacity and $C_{P}$ is the heat capacity. The particle diffusivity is derived from the diffusion coefficients of the various polluting metals that determine the molecular weight of a fluid element.

The mixing that results from thermohaline instability can be approximated with a coefficient that scales with thermal diffusivity and the molecular weight gradient (Kippenhahn et al., 1980):

[TABLE]

where $\alpha_{\rm th}$ is a dimensionless efficiency parameter. In Paper I, we explored inferences for polluted WD accretion rates using the mixing treatment of Equation (7) with $\alpha_{\rm th}=1$. Note that this mixing treatment does not explicitly check the criterion for instability given in Equation (5), but in Paper I we verified that it is satisfied for regions of interest in polluted WDs where thermohaline mixing may occur.

MESA also offers a thermohaline mixing treatment based on the work of Brown et al. (2013), which is calibrated against their 3D hydrodynamic simulations. This treatment explicitly accounts for the criterion in Equation (5) and produces a mixing coefficient designed to scale smoothly to zero as conditions approach the limit defined there. Figure 5 shows the surface Ca mass fraction in polluted $0.6\,M_{\odot}$ WD models after accreting for many diffusion timescales, with thermohaline mixing according to either Equation (7) or Brown et al. (2013). These MESA models also include element diffusion at all times. Figure 5 shows that results for polluted WD models using thermohaline mixing based on Brown et al. (2013) are qualitatively similar to those using Equation (7) with $1\lesssim\alpha_{\rm th}\lesssim 10$.

For small accretion rates, thermohaline mixing is not active, and the equilibrium surface mass fractions shown in Figure 5 match the prediction of Equation (2) for the diffusion timescales given in Table 2. Above a critical accretion rate $\dot{M}_{\rm crit}$, the metal concentration at the surface builds up a sufficient magnitude of $\nabla_{\mu}$ to excite thermohaline instability, and MESA models including thermohaline mixing diverge from the prediction of Equation (2). For a ${T_{\rm eff}=11,000\,\rm K}$ WD, Figure 5 shows that this critical rate is around $\dot{M}_{\rm crit}\approx 10^{6}\,\rm g\,s^{-1}$. Table 3 gives values of this critical rate for MESA WD models over a range of $T_{\rm eff}$. For models with ${T_{\rm eff}>12,000\,\rm K}$, the surface convection zones are so small that thermohaline mixing is active for all accretion rates in the range that we explored ($\dot{M}_{\rm acc}>10^{4}\,\rm g\,s^{-1}$).

The curve shown for the Brown et al. (2013) prescription in Figure 5 varies slightly from the similar plot shown in Figure 3 of Paper I. This is due to a small correction to the MESA implementation of this routine that affects the mixing coefficient in the regime near the limit of thermohaline instability. This correction was introduced after MESA release version 10398, which was used for Paper I, but it is present in MESA version 11191, which we use for all models that include thermohaline mixing in this paper. The asymptotic analysis regimes presented in Appendix B of Brown et al. (2013) form the basis of the 1D mixing treatment. In particular, their Appendix B.3 addresses the regime in which the fluid is near the limit imposed by Equation (5). The method relies on an expansion in the parameter

[TABLE]

which is assumed to be small. The implementation for this regime in MESA version 11191 ensures that this parameter is sufficiently small whenever applying the method of Brown et al. (2013) Appendix B.3, yielding more consistent results than version 10398. With these corrections, the MESA implementation shows more mixing near the boundary of thermohaline instability defined by Equation (5). Hence, models employing the Brown et al. (2013) routine in MESA version 11191 diverge from the prediction of diffusion alone at the lower accretion rates seen in Figure 5.

Using this updated thermohaline mixing treatment, we construct a large grid of accreting DA WD models as in Paper I. Effective temperatures of the models span the range ${6,000\,{\rm K}<T_{\rm eff}<20,500\,\rm K}$, and accretion rates for each temperature span $10^{4}\,{\rm g\,s^{-1}}<\dot{M}_{\rm acc}<10^{12}\,\rm g\,s^{-1}$. All models accrete bulk earth material (McDonough, 2001). We tabulate values of $X_{\rm Ca}$ present at the surface of each model after 100 diffusion timescales as defined by Table 2. We then interpolate on these tables to map observed values of $X_{\rm Ca}$ to total inferred accretion rates $\dot{M}_{\rm acc}$. We also expand upon the results of Paper I by providing these tables for models with three different WD masses to allow interpolation in $\log g$: $M_{\rm WD}/M_{\odot}=0.38,0.60,0.90$ (${\log g\approx 7.5,8.0,8.5}$). These tables are available along with simple python interpolation routines at https://doi.org/10.5281/zenodo.2541235 (Bauer, 2019).

Figure 6 shows inferred accretion rates based on these tables for the same sample of polluted DA WDs (Koester & Wilken, 2006) that was discussed in Paper I. In general, accretion rates are similar to the inferences made in Paper I, though the very highest inferences are slightly lower than the previous highest values. A few WDs also show adjustments due to observed values of $\log g$ different from the value of $8.0$ assumed in Paper I, especially for ${9,000{\,\rm K}<T_{\rm eff}<13,000\,\rm K}$, where the surface convection zone masses are especially sensitive to $\log g$ (see Figure 1). However, the overall qualitative picture remains the same. Thermohaline mixing causes inferred accretion rates to increase by several orders of magnitude for WDs with ${T_{\rm eff}>10,000\,\rm K}$!

4.1.1 Non-constant Accretion Rates

The previous section shows results when accretion occurs in a steady state for many diffusion timescales, allowing the surface metal pollution to approach equilibrium abundances. However, if the source of accretion supplied to the surface varies with time, this can introduce complexities in the $\nabla_{\mu}$ profile that sets the conditions for thermohaline instability according to Equations (4) and (5). In particular, heavy elements must be continually supplied to the surface to maintain $\nabla_{\mu}<0$ in the mixing region relevant to observable pollution. If the accretion rate decreases significantly, the gradient necessary for thermohaline instability can disappear, halting thermohaline mixing.

As a simple illustration, we show in Figure 7 a MESA model including thermohaline mixing that accretes at a constant rate until it approaches equilibrium, followed by a cessation of accretion after $10\,{\rm yr}\approx 100\tau_{\rm diff}$. Due to the sudden disappearance of an inverted $\nabla_{\mu}$ in the surface region governing observable metal pollution, thermohaline mixing is no longer relevant. Instead, the diffusion timescales of Table 2 dictate the fast exponential decay of metal pollution at the photosphere. These results contrast with the diffusion-only MESA model shown in Figure 3, where the same diffusion timescale governs both the approach to equilibrium and exponential decay after accretion ceases.

Figure 7 also demonstrates important features involving differentiation of the accreted composition. When thermohaline mixing is active near the surface during the constant accretion phase, no composition differentiation occurs because fluid elements that dominate the mixing transport all elements together. However, once thermohaline mixing ceases near the surface, individual particle diffusion dominates, and significant differentiation quickly occurs within a few diffusion timescales. The middle panel of Figure 7 shows that the deeper layers where thermohaline mixing is still active reflect the accreted bulk earth composition (McDonough, 2001), but separate diffusion timescales for each element quickly rearrange the surface composition. Elements with the shortest diffusion timescales such as ${}^{56}{\rm Fe}$ disappear from the surface much sooner, even when they were previously among the most abundant due to the accreted composition.

4.1.2 Helium-dominated Atmospheres

WDs with helium-dominated atmospheres do not experience the same corrections due to thermohaline mixing that hydrogen-dominated atmospheres do. Two effects conspire to greatly reduce the potential for a $\nabla_{\mu}$ large enough to excite thermohaline instability. First, the mean molecular weight of the dominant background material (He) is more than double that in the case of a hydrogen atmosphere, so the contrast with accreting metals is not as severe. Second, surface convection zones for helium atmospheres contain much more mass than hydrogen at a given temperature (Koester, 2009). This dilutes accreted metals and prevents the buildup of a significant $\nabla_{\mu}$ below the convection zone.

For example, we constructed a $0.59\,M_{\odot}$ MESA WD model with $T_{\rm eff}=18,000\,\rm K$ and a pure He atmosphere. We found that the surface convection zone mass of this model was $M_{\rm cvz}=8\times 10^{-8}\,M_{\odot}$, and diffusion timescales for accreted metals were on the order of $10^{5}$ years. These values agree with the tables of Koester (2009) for $\log g=8.0$ DB WDs. We explored MESA runs for this WD model accreting bulk earth composition at rates in the range $10^{4}\,{\rm g\,s^{-1}}<\dot{M}_{\rm acc}<10^{12}\,\rm g\,s^{-1}$. We included thermohaline mixing in the runs using the treatment of Equation (7) with $\alpha_{\rm th}=10$ (the MESA treatment based on Brown et al. (2013) is not applicable here because it assumes a hydrogen-dominated background). Even for the highest accretion rates, we find adjustments of at most one order of magnitude to inferred accretion rates compared to calculations that assume no thermohaline mixing (Figure 8). Figure 6 shows that typical accretion rates inferred for DB WDs in this temperature range are $10^{8}$–$10^{10}\,\rm g\,s^{-1}$, and our MESA models show negligible corrections due to thermohaline mixing in this regime. The surface convection zone grows up to three orders of magnitude larger for cooler WDs (Koester, 2009), and the largest rates inferred for DB WDs only reach $10^{11}\,\rm g\,s^{-1}$, so thermohaline mixing will be inconsequential for He-atmosphere WDs with $T_{\rm eff}\lesssim 18,000\,\rm K$.

Below $T_{\rm eff}\approx 16,000\,\rm K$, the cool, dense, neutral helium at the surface of the WD falls outside the regime covered by opacity tables currently available in MESA (Paxton et al., 2011). The code is therefore not able to set a physical outer boundary condition for the model below this temperature. Tabulated outer boundary conditions such as those used in the case of hydrogen-dominated atmospheres (Section 2) have been implemented in other WD codes (e.g., Camisassa et al., 2017), but no such option is currently available in MESA. In the context of polluted WDs, it is also unclear whether atmosphere conditions tabulated for pure helium would be sufficient, since opacity may be sensitive to contaminating metals through effects such as He- free-free absorption. Without the ability to set an appropriate outer boundary condition, MESA models cannot give reliable structures for the outer layers and depths of surface convection zones. A more thorough investigation of polluted WDs with helium-dominated atmospheres in MESA awaits extensions to atmosphere capabilities that can account for these issues.

4.1.3 Rotation

Rotational mixing and its interplay with other fluid processes can be important in stars (Sengupta & Garaud, 2018). This potential impact is quantified with the Rossby number ${\rm Ro}=U/2\Omega L$, where $\Omega$ is the rotational frequency, and $U$ and $L$ are the characteristic velocity and length scale for the relevant fluid process. Large values of the Rossby number indicate that rotation is not expected to have a strong influence, while ${\rm Ro}\lesssim 1$ indicates potential for significant modifications. Sengupta & Garaud (2018) studied the effect of rotation on thermohaline mixing in stellar interiors, where they derived the Rossby number in an actively mixing region as

[TABLE]

where $N$ is the Brunt-Väisälä frequency. In a non-degenerate WD atmosphere, this frequency is of order $N^{2}\sim g/H$, where ${H=k_{\rm B}T/m_{\rm p}g}$ is the local pressure scale height. For our polluted WD models experiencing moderate amounts of thermohaline mixing, we estimate $(\varphi/\delta)\nabla_{\mu}/(\nabla_{T}-\nabla_{\rm ad})\sim 10^{-4}$ (cf. Paper I). We can therefore rewrite Equation (9) in terms of the critical rotation rate $\Omega_{\rm crit}=\sqrt{GM/R^{3}}$ as

[TABLE]

This requires $\Omega/\Omega_{\rm crit}\gtrsim 10^{-2}$ for rotation to be important (${\rm Ro}\lesssim 1$). However, typical isolated WD rotation periods are around one day (Hermes et al., 2017), while the critical rotation period is on the order of a few seconds, so we do not expect rotation to influence the thermohaline mixing in typical polluted WDs.

Thermohaline mixing has also been discussed as a mechanism for explaining observed surface abundances in low-mass giant stars (Charbonnel & Zahn, 2007; Denissenkov & Pinsonneault, 2008; Cantiello & Langer, 2010), but this may require an implausibly large mixing efficiency $\alpha_{\rm th}>100$ for implementations such as Equation (7). This mixing efficiency appears to be inconsistent with the mixing found in our models based on Brown et al. (2013). Sengupta & Garaud (2018) suggested that the interplay of rotation with thermohaline instability may enhance mixing near the cores of giant stars. Since we estimate that rotation would be irrelevant for thermohaline mixing in polluted WDs, this may alleviate the apparent tension between thermohaline mixing efficiency inferred in these different regimes.

4.2 Overshoot

Convective overshoot may cause fluid motions that can keep composition thoroughly mixed well below the formal boundary for convective instability according to the Ledoux criterion (Freytag et al., 1996; Koester, 2009; Tremblay et al., 2015), even in the absence of thermohaline instability. This will lead to a larger effectively mixed region and longer diffusion timescales for a given $T_{\rm eff}$ (Brassard & Fontaine, 2015; Tremblay et al., 2017).

To estimate mixing due to overshoot beneath the convective zone, we follow the results of Tremblay et al. (2015) and use a diffusion coefficient that decays exponentially with pressure scale height:

[TABLE]

where $r_{0}$ is the radial coordinate of the base of the convection zone, $H_{P}$ is the pressure scale height there, and $D_{0}$ is the mixing coefficient from MLT near that location. Figure 9 shows the resulting diffusion coefficient profiles for two MESA models. In the absence of thermohaline mixing, this will lead to a new mass of the fully mixed surface region ($M_{\rm ovr}$) defined by the location where the overshoot mixing decays to where element diffusion takes over ($D_{\rm overshoot}<D_{\rm diffusion}$, see Figure 9). There is then a corresponding new diffusion timescale for each element

[TABLE]

where $\rho$, $r$, and $v_{{\rm ovr},i}$ are all evaluated at the base of the new mixing region defined by $M_{\rm ovr}$. The equilibrium observable abundance of an accreted element will then be

[TABLE]

instead of the analogous value given in Equation (2).

Figure 10 shows how observable abundances of accreting metals change when including this form of overshoot in our MESA models. For models at $T_{\rm eff}=11,000\,\rm K$, the new diffusion timescale for Ca is $\tau_{\rm ovr}=1300\,\rm yr$, almost $100$ times larger than the timescale without overshoot. The larger mixing region means that accreted metals are more diluted for a given accretion rate, and so larger accretion rates are needed for thermohaline mixing to cause the observable abundances to diverge from the prediction of Equation (13). Still, for accretion rates of $\dot{M}_{\rm acc}>10^{7}\,\rm g\,s^{-1}$, thermohaline mixing begins to dominate the final observed abundance, and overshoot causes only small adjustments when thermohaline mixing is active (see also the left panel of Figure 9). For the case of a $T_{\rm eff}=14,000\,\rm K$ WD shown in Figure 11, overshoot extends the small surface mixing region to $M_{\rm ovr}\approx 8\times 10^{-16}\,M_{\odot}$ (see right panel in Figure 9), but this is still so thin that thermohaline mixing dominates even for modest accretion rates.

While the results shown in this section may serve as a useful qualitative description of effects that can be expected from overshoot, it is likely that the overshoot mixing prescription given in Equation (11) is too simplistic for WD pollution applications. Simulations are beginning to probe regimes specific to convective overshoot in WDs (Montgomery & Kupka, 2004; Tremblay et al., 2015; Kupka et al., 2018), and they appear to show that a simple exponential decay in the diffusion coefficient is only accurate within a few scale heights of the convective boundary. Extrapolation down to the much smaller diffusion coefficients relevant for particle diffusion is likely inaccurate. Simulations by Lecoanet et al. (2016) found overshoot mixing that decays with a Gaussian profile ($D_{\rm overshoot}\propto\exp[-(r-r_{0})^{2}/2H^{2}]$ for some scale height $H$) rather than the exponential of Equation (11). A few more recent results appear to confirm this Gaussian overshoot profile in other contexts (Jones et al., 2017; Korre et al., 2018). This faster decay of the diffusion coefficient would imply that the extra extent of overshoot mixing is smaller than what is shown in Figure 9. We therefore refrain from a complete exploration of MESA models including overshoot until simulations can provide better constraints on overshoot mixing well below convective boundaries.

Our results are sufficient to conclude that overshoot will have negligible effects on most accretion rate inferences for $T_{\rm eff}\gtrsim 12,000\,\rm K$, where thin surface mixing regions result in strong concentrations of metals that make thermohaline mixing dominant. For lower temperatures, Figure 10 suggests that overshoot may cause significant adjustments to accretion rate inferences in cases where thermohaline mixing is not active. Even at higher temperatures, the new timescales due to overshoot may be important for decay phases where supply of fresh accreted material has ended and there is nothing to maintain the $\nabla_{\mu}$ needed to drive thermohaline instability near the surface mixing region. In this case, the $M_{\rm ovr}$ and $\tau_{\rm ovr}$ parameters will govern the exponential decay of observable surface abundances.

5 Discussion

A significant fraction of WDs show evidence of pollution (Koester et al., 2014), and if this fraction represents the fraction of the lifetime of each individual WD that it is polluted, then Figure 6 may be taken as approximately showing a complete history of accretion rates experienced over a WD lifetime. In this case, the total mass of planetesimal material accreted over a WD lifetime would be dominated by the high rates experienced by young WDs, yielding a high estimate of ${M_{\rm tot}\sim(3\times 10^{8}\,{\rm yr})(10^{12}\,{\rm g\,s^{-1}})\sim 10^{28}\,\rm g}$. However, Koester et al. (2014) point out that the Ca based sample used to construct Figure 6 may be biased toward objects that are especially heavily polluted, since the optical Ca lines used to select this sample require higher Ca abundances to be detectable for $T_{\rm eff}\gtrsim 15,000\,\rm K$ as total WD flux moves primarily into the UV. While this is unlikely to change the accretion rates inferred for the objects shown in Figure 6, it could hide a much larger intrinsic scatter in the accretion rates for young DA WDs. Therefore, $10^{28}\,\rm g$ could be an overestimate of the total mass accreted over a WD lifetime.

This sample of polluted WDs may reveal that some young WDs are undergoing short timescale bursts of accretion such as those suggested by Rafikov (2011b); Metzger et al. (2012). This may help explain the discrepancy with DBZ WDs for $T_{\rm eff}\gtrsim 15,000\,\rm K$. The rates here are too high to be explained by Poynting-Robertson drag (Rafikov, 2011a), but rare runaway bursts would leave very different observational signatures for DA and DB WDs (Farihi et al., 2012). DA WDs approach a quasi-equilibrium surface abundance within days or years in this temperature range even for our MESA models including thermohaline mixing. On the other hand, the diffusion timescales in DB WDs are of order $10^{5}$-$10^{6}$ years, and bursts lasting less than $10^{4}$ years would never approach an equilibrium surface pollution level suggesting a high rate. Instead, DB WD surfaces may represent a more accurate estimate of accretion rates averaged over their much longer diffusion timescales. A more conservative mass estimate for total planetesimal material may then be ${M_{\rm tot}\sim(3\times 10^{8}\,{\rm yr})(10^{10}\,{\rm g\,s^{-1}})\sim 10^{26}\,\rm g}$.

Alternative processes could supply polluting material for longer timescales at rates higher than the limits of Poynting-Robertson drag, e.g., collisional cascades (Kenyon & Bromley, 2017a, b) or viscous evolution of earth-mass dust disks (van Lieshout et al., 2018). Hence, short bursts are not strictly necessary to explain the rates shown in Figure 6, but longer timescale processes may then require that the planetesimal environments form with significantly different amounts of mass around DA and DB WDs. Wyatt et al. (2014) found that stochastic accretion of a distribution of planetesimal sizes may be able to explain some discrepancies in inferred accretion rates for DA and DB WDs without the need to appeal to large bursts, but this analysis assumed accretion rates inferred without accounting for thermohaline mixing.

Finally, we note that some authors have pointed out trends of inferred accretion rates that decline with WD age over timescales of Gyr (e.g., Hollands et al., 2018; Chen et al., 2018), consistent with slow depletion of the planetesimal reservoirs that can obtain highly eccentric orbits on which they will eventually be tidally disrupted (e.g., Debes et al., 2012; Mustill et al., 2018). Our results appear to suggest that this decline may be more dramatic during the first Gyr of evolution when thermohaline mixing is accounted for. In particular, the broken power-law for accretion rates over time used by Chen et al. (2018) may not be necessary for rates inferred using our MESA models. Instead, a single power-law may work for all WD ages, consistent with the rate at which asteroids dynamically encounter the WD in the model of Chen et al. (2018).

6 Conclusions

We have confirmed the result of Paper I that thermohaline mixing in polluted DA WDs with $T_{\rm eff}\gtrsim 10,000\,\rm K$ requires accretion rates several orders of magnitude larger than calculations assuming only gravitational sedimentation. We have provided results from an expanded grid of models to allow interpolation in $\log g$ as well as $T_{\rm eff}$ (https://doi.org/10.5281/zenodo.2541235, Bauer, 2019). We also find that thermohaline mixing is inconsequential in polluted DB WDs with $T_{\rm eff}\lesssim 18,000\,\rm K$ due to much more massive surface convection zones. Polluted DA WDs experience a regime of accretion rates low enough that thermohaline mixing is not active (Table 3), and so Table 2 provides diffusion timescales based on our MESA models. These timescales are also relevant for WDs where accretion is no longer ongoing, as they govern the exponential decay of metals sinking away from the surface where thermohaline mixing is no longer active, even when it was active during accretion. Finally, we have also provided a qualitative description of the effects of convective overshoot, though we refrain from a full exploration of its effects due to quantitative uncertainty in the overall extent of overshoot. However, we note that for WDs with thin surface convection zones ($T_{\rm eff}\gtrsim 12,000\,\rm K$), thermohaline mixing dominates down to layers deeper than overshoot can extend, and hence we do not expect significant modifications to inferred accretion rates in this regime.

Acknowledgments: We are grateful to the anonymous referee for comments that improved the clarity of the paper. We are grateful to Matteo Cantiello, Tim Cunningham, Jay Farihi, Gilles Fontaine, Boris Gänsicke, Pascale Garaud, Marc Pinsonneault, Josiah Schwab, Sutirtha Sengupta, Andrew Swan, and Pier-Emmanuel Tremblay for inspiring and encouraging discussions relating to many topics in this paper. We thank Bill Paxton for continuous efforts in support of broad MESA usage. This work was supported by the National Science Foundation through grants PHY 17-148958 and ACI 16-63688.

Appendix A Diffusion Coefficients for Neutral Atoms

When diffusion occurs near the surface of a star, even a very small fraction of neutral particles for a given species can have a dramatic effect on the net diffusion flux for the element. The Coulomb collision formalism for diffusion coefficients no longer applies for neutral atoms, and induced dipole scattering of the neutral atoms with background ions becomes the relevant physical process for diffusion (Vennes et al., 2011a, b). MESA does not currently offer options for diffusion coefficients based on dipole scattering, but Section 9 of Paxton et al. (2015) describes much of the formalism necessary to construct them. Referring to Equations (84)–(86) of Paxton et al. (2015), we see that the relevant coefficients for the Burgers (1969) diffusion equations are expressed in terms of

[TABLE]

where $\alpha_{st}=2k_{B}T/\mu_{st}$, $\mu_{st}=m_{s}m_{t}/(m_{s}+m_{t})$ is the reduced mass of the scattering particles, and

[TABLE]

are the traditional scattering cross section integrals. The scattering angle $\chi_{st}$ is a function of both impact parameter $b$ and relative velocity $v$, depending on the physics of the scattering process between particles $s$ and $t$. For dipole scattering, Chapter 2 in Draine (2011) gives the scattering cross section of an ion with a polarizable neutral atom as

[TABLE]

where $Z$ is the charge of the ion and $\alpha_{N}$ is the polarizability of the neutral atom. The result for $\Sigma^{(11)}_{st}$ is then

[TABLE]

Hence the result for the resistance coefficient for ions scattering with induced dipoles of neutral atoms is

[TABLE]

For comparison, the resistance coefficient given by Burgers (1969) for diffusion of ions with other ions is

[TABLE]

where $\ln\Lambda$ is the Coulomb logarithm. Diffusion velocities scale inversely with the resistance coefficients $K_{st}$, so we can use these expressions to estimate the relative speeds of neutral and ionized metals.

Consider the case of a mixture of singly ionized oxygen and neutral oxygen diffusing in an ionized Hydrogen background. The above equations give that the ratio of the ion-neutral coefficient to the ion-ion coefficient is

[TABLE]

The polarizability of neutral oxygen is given in Table 2.1 of Draine (2011) as $\alpha_{N}=5.326a_{0}^{3}$, where $a_{0}$ is the Bohr radius. For $T=15,000\,\rm K$, we find that the resistance coefficient ratio is

[TABLE]

For conditions near the surface of a WD, we estimate the Coulomb logarithm as $\ln\Lambda\approx 5$ (Spitzer, 1962, Table 5.1). The final result for the ratio of coefficients is $K_{st}^{(\rm dipole)}/K_{st}^{(\rm coulomb)}\sim 10^{-2}$. This means that the neutral oxygen particles will have diffusion velocities approximately $100$ times faster than the singly ionized particles. A mere $1\%$ of particles being neutral for a particular element can therefore significantly modify the net diffusion flux for that element. This effect would be most noticeable for metals with relatively high first ionization potentials.

Since the physics of scattering is fundamentally different for charged and neutral particles, the diffusion and resistance coefficients do not scale in a meaningful way as $Z\to 0$. It is hence meaningless to adopt an average charge for an element for purposes of diffusion calculations in the case of $\bar{Z}<1$. This is why the diffusion implementation in MESA currently assumes that all diffusing particles are at least singly ionized. Extending the implementation to account for neutral particles would require two improvements: a) the ability to separate neutral particles off as distinct diffusion classes, and b) incorporating tables of atomic polarizabilities to use with Equation (A5) for the resistance coefficients to use in the Burgers equations.