Horizontal spreading of planetary debris accreted by white dwarfs

White dwarfs with metal-polluted atmospheres have been studied widely in the context of the accretion of rocky debris from evolved planetary systems. One open question is the geometry of accretion and how material arrives and mixes in the white dwarf surface layers. Using the 3D radiation-hydrodynamics code CO$^5$BOLD, we present the first transport coefficients in degenerate star atmospheres which describe the advection-diffusion of a passive scalar across the surface-plane. We couple newly derived horizontal diffusion coefficients with previously published vertical diffusion coefficients to provide theoretical constraints on surface spreading of metals in white dwarfs. Our grid of 3D simulations probes the vast majority of the parameter space of convective white dwarfs, with pure-hydrogen atmospheres in the effective temperature range 6000-18000 K and pure-helium atmospheres in the range 12000-34000 K. Our results suggest that warm hydrogen-rich atmospheres (DA; $\gtrsim$13000 K) and helium-rich atmospheres (DB, DBA; $\gtrsim$30000 K) are unable to efficiently spread the accreted metals across their surface, regardless of the time dependence of accretion. This result may be at odds with the current non-detection of surface abundance variations at white dwarfs with debris discs. For cooler hydrogen- and helium-rich atmospheres, we predict a largely homogeneous distribution of metals across the surface within a vertical diffusion timescale. This is typically less than 0.1 per cent of disc lifetime estimates, a quantity which is revisited in this paper using the overshoot results. These results have relevance for studies of the bulk composition of evolved planetary systems and models of accretion disc physics.


INTRODUCTION
It is now accepted that 25-50 per cent of all white dwarfs show clear evidence of metal pollution in their surface layers (Zuckerman et al. 2010;Koester et al. 2014;Farihi et al. 2016). For most of these objects the unequivocal origin of this material is the evolved planetary system hosted by the white dwarf (Jura 2003;Gänsicke et al. 2012;Vanderburg et al. 2015;Gänsicke et al. 2019;Manser et al. 2019;Vanderbosch et al. 2020;Harrison et al. 2018;Bonsor et al. 2020;Xu et al. 2019). The canonical model is that a circumstellar disc borne of one or many tidally disrupted planetesimals is able to constantly feed planetary material onto the white dwarf surface (Jura 2003). These metals disappear from the observable layers at rates defined by the sinking timescale (Koester 2009). In accretiondiffusion equilibrium, the atmospheric abundances remain constant and decrease once accretion stops (Dupuis et al. 1992). Other models suggest that direct impacts of asteroids on star-grazing orbits may ★ E-mail: timothy.cunningham@warwick.ac.uk † STFC Ernest Rutherford Fellow provide a small percentage (∼1-2 per cent) of photospheric metal pollution (Wyatt et al. 2014;Brown et al. 2017). In all cases, surviving planets are required to gravitationally scatter asteroids and minor planets onto the central regions (Debes & Sigurdsson 2002;Debes et al. 2012;Veras 2016).
There have been significant theoretical efforts focused on accretion onto white dwarfs. Of particular importance, Wyatt et al. (2014) provided a theoretical constraint on the size of particles reaching the white dwarf surface and Turner & Wyatt (2020) found that the singlebody tidal disruption assumption is invalid in over 20 per cent of cases. Regarding the geometric configuration of accretion, Metzger et al. (2012) provided a comprehensive discussion regarding models of discs, finding that magnetic fields as weak as 0.1-1 kG at the white dwarf surface can affect the accreted material, potentially directing it toward the poles. A predicted outcome of this pole-on accretion is the production of X-rays. This is a detection yet to be made for metal-polluted white dwarfs (Farihi et al. 2018).
Accretion in cataclysmic variables (CVs) is a more mature field than that of metal-polluted white dwarfs. For weakly magnetic CVs there are models which favour near-spherical deposition (Patterson & Raymond 1985) for low accretion rates ( < 10 16 g s −1 ) and equatorial deposition (Piro & Bildsten 2004) for high accretion rates ( > 10 16 g s −1 ). The regime of metal-polluted white dwarfs, where the accretion rates are many orders of magnitude lower ( ∼ 10 9 g s −1 , see, e.g., Farihi et al. 2016), has however not been widely explored.
The sinking timescales of trace elements through white dwarf atmospheres has received, and continues to receive, considerable attention (Schatzman 1945;Paquette et al. 1986a,b;Dupuis et al. 1992;Koester 2009;Bauer & Bildsten 2019;Cunningham et al. 2019;Koester et al. 2020;Heinonen et al. 2020). White dwarfs develop convective atmospheres and envelopes below 18 000 K (hydrogenrich; Cunningham et al. 2019) or 40-60 000 K (helium-rich; Bergeron et al. 2011;Cukanovaite et al. 2019). In the accretion-diffusion model (Dupuis et al. 1992;Koester 2009) it is normally assumed that for white dwarfs with a convective atmosphere, the observed spectroscopic abundance is distributed homogeneously throughout the mixed region, in the vertical and horizontal directions.
Rapid convective mixing in the vertical (or radial) direction is welljustified given the shallow convection zones restricted to the hydrogen or helium envelopes (depth of 0.1-100 km) and short turnover timescales in the range c = 0.1-1000 s (Van Grootel et al. 2012). On the other hand, the geometric extent of convection zones in the horizontal (or surface-plane, or azimuthal) directions can be many orders of magnitude larger; with a typical white dwarf circumference of 2 WD ∼ 50 000 km. Thus the assumption that convection is able to evenly distribute metals over the entire stellar surface is less obvious, especially if the metals sink vertically before they are able to spread. A deviation from this assumption could impact the measured photospheric abundances, and hence the accretion rates and accreted masses at these systems. A heterogeneous distribution of surface metals may also produce variations in the spectroscopic line strengths and equivalent widths which may be detectable on the rotational period of the white dwarf.
A theoretical investigation of metal transport across the surface of a white dwarf was made by Montgomery et al. (2008) who sought to constrain the homogeneity of metals distributed on the surface of the ZZ Ceti G29-38. The authors argued that sinking was faster than lateral spreading and suggested that any observed variability could provide a constraint on the distribution of metals at the surface. However, no variations in metal line strength have been robustly identified in G29-38 (Debes & López-Morales 2008;Reach et al. 2009), nor any other white dwarf so far .
We revisit the question of sinking versus lateral spreading with 3D radiation-hydrodynamics (RHD) which has recently led to a better understanding of the convective atmospheres of white dwarfs (Tremblay et al. 2013a;Kupka et al. 2018) and enhanced numerical techniques for characterising the advective transport of material (Cunningham et al. 2019). In particular, we have recently shown that convective overshoot can increase the sinking times in hydrogenatmosphere white dwarfs by two orders of magnitude. This will directly impact any conclusion about the balance of sinking times and lateral spreading times.
In Section 2 we describe our numerical approach to the problem. In Section 3 we present results from our 3D model atmosphere simulations. Sections 4 and 5 focus on the implications of our results on the white dwarf population and finally we provide conclusions in Section 6.

NUMERICAL SETUP
Our approach is to model the transport of trace metals across the white dwarf surface as an effective diffusion process. This requires the local diffusion coefficient describing mixing in the radial and angular directions. With a vertical extent of less than ∼10 km, a white dwarf convection zone is typically much less than a per cent of the stellar radius, and thus we can simplify the problem to mixing in the radial and horizontal, or surface-plane, directions. In the next two sections we describe first the derivation of the effective diffusion coefficients and secondly the incorporation of these diffusion coefficients into a model describing the surface spreading.

Diffusion coefficients
We present the first 3D RHD CO 5 BOLD (Freytag et al. 2012) simulations designed to directly study the horizontal spreading of trace material in the surface layers of white dwarf stars. Our simulations were developed by taking well-relaxed models from two existing grids, adding tracer densities to quantify the horizontal advectiondiffusion and running these as a new set of simulations. A grid of deep, RHD simulations for the surface hydrogen convection zones of white dwarfs in the effective temperature range 11 400-18 000 K were presented in Cunningham et al. (2019, hereafter CU19). Earlier, Tremblay et al. (2013a) developed a grid of RHD simulations for the surface convection zones of cooler DA white dwarfs with effective temperatures between 6000-13 500 K. Here we combine both grids such that the full parameter space of convective DA white dwarfs is probed: from when convective instabilities initially arise at 18 000 K down to when the convection zone is sufficiently deep that it couples with the degenerate core. For comparison, we show some key physical properties of both grids in Fig. 1. The model grid from CU19 is shown in blue and the model grid from Tremblay et al. (2013a,c) is shown in orange. The horizontal diffusion coefficients derived from the evolution of a tracer-density are used to confirm that the diffusion coefficient can be computed from mean 3D parameters -an output of the simulations -providing a computational advantage. The tracer-density experiments were limited to pure-hydrogen models only, but the parametrisation using mean 3D properties allowed us to exploit the grid of 3D DB (pure-helium) convection zone models from Cukanovaite et al. (2019) -spanning effective temperatures of 12 000-34 000 K and surface gravities (log ) of 7.5-9.0 -to derive horizontal diffusion coefficients for pure-helium atmospheres.
The CO 5 BOLD simulations presented here, and those previously published, were performed on a Cartesian grid using the box-in-a-star setup with periodic side boundary conditions, an open top boundary and bottom boundary that is either closed or open to convective flows but always open to radiative flux. The horizontal diffusion coefficients were derived using a tracer density to study the local flows in 3D RHD simulations and via mean 3D parameters from relaxed simulations. The length scale over which the horizontal diffusion coefficient is calculated is typically an order of magnitude less than the extent of the simulation box and thus the presence of periodic horizontal boundaries is not expected to influence the determination of the horizontal diffusion coefficient. To test this hypothesis, in Section 3.1 we perform a convergence test with a simulation domain four times larger in area.
Extending the working of CU19, we trace the effective diffusion in the surface layers of white dwarfs by following the local advective flows in relaxed 3D simulations of surface convection zones. A simulation is deemed relaxed when the mean properties (i.e., effective temperature, horizontally-averaged velocities) reach a quasi-steady state. This was previously confirmed for all simulations used in this study. This ensures that the inferred quantities (i.e., horizontal diffusion coefficient) are derived for a simulation which best represents the physical characteristics of a white dwarf at a given effective temperature and surface gravity.
The horizontal diffusion coefficient is calculated by adding vertical slabs of a passive scalar tracer density throughout the simulation domain. In order to probe diffusion in a given direction, a slab must be oriented such that the surface normal is parallel to said direction, i.e., to probe the diffusion in the -direction, a slab must be placed in the − plane. Multiple slabs (≈24-28) are initialised at evenly spaced locations in the -direction. This allows to quantify the uncertainty on diffusion in the -direction. The spatial extent of these ensembles, if undergoing a diffusion process, should evolve proportionally to the square-root of time ( 1/2 ), such that the diffusion coefficient in the -dimension, , can be calculated as 2 log 10 ( * ( , )) = log 10 ( ) + log 10 (2 ) , where * is the density-weighted standard deviation of the tracers defined as where is the number of cells in the dimension, t ( ) is the tracer density mean-averaged across the and dimensions, and * is the density-weighted mean -position of the tracer density given by The initial condition for each independent slab of tracer density added at = is defined as As the tracer density, , is a passive scalar the actual number densities only bare relevance to the precision of the numerical reconstruction schemes and rounding errors. The initial positions, , are chosen to evenly sample the full horizontal extent of the simulation box in the -dimension. This procedure is repeated in the -dimension such that an alternative component to the horizontal diffusion coefficient, , is found by placing slabs of a passive scalar density in the − plane. As it is not expected that the horizontal diffusion coefficient should vary across the -and -dimensions the final diffusion coefficient is computed as a mean across all slabs, such that where and are the number of slabs added in the -anddimensions, respectively. The standard deviation around this mean quantity (Eq. 5) is chosen to represent the error attributed to each horizontal diffusion coefficient. The results of these computations are considered in Section 3.1 and compared to the well-known 1D estimation of the diffusion coefficient. We note that this procedure ultimately provides a single diffusion coefficient describing transport in the surface plane. This omits any radial dependence of the diffusion coefficient, and precludes any differential rotation effects.
In 1D mixing-length theory the local approximation of the diffusion coefficient is typically given (e.g., Zahn 1991;Ventura et al. 1998 where d is a characteristic convective length scale and conv is the mean convective velocity. Pressure scale height, P , has often been proposed to be a suitable choice of characteristic mixing length.
Here we explore the validity of this approximation by comparing this expression to results from tracer-density evolution in the convection zone of white dwarfs at nine effective temperatures spanning the full range of convective DA white dwarfs: 6000-18 000 K. The vertical diffusion coefficients for 11 400-18 000 K were presented in CU19. In Section 3 we present the results of these experiments, providing the first horizontal diffusion coefficients derived from multidimensional simulations in white dwarfs.

Simulating surface transport
Here we introduce the method by which we model the global transport of metals across the surface of the white dwarf. We begin by taking the limiting case of a delta function accretion, both in the temporal and spatial domains. Physically, this effectively represents the collision of an asteroid with the surface of the white dwarf. Modelling the surface of the white dwarf as a ring allows us to exploit the wellknown diffusion equation: Here, ( , ) is the local concentration of some scalar which, in our case, represents the accreted heavy elements, and is the diffusion coefficient. As a preliminary scenario we consider the situation that metals are not permitted to sink out of the observable layers. This assumption allows for a simple analytic form of position dependent abundance to be written down. We also assume throughout this study that the diffusion coefficient is independent of the position across the surface. In this case, with only surface (horizontal) transport and a delta function accretion the local concentration can be written as where corresponds to the position along the circumference of the white dwarf and surf is the diffusion coefficient describing the effectiveness of material transport in the plane of the surface. This diffusion coefficient is assumed to have no depth dependence for the time being. For spatial domain, , contained in interval [− , ] and with initial and boundary conditions where denotes the derivative of with respect to , we can write the analytic solution to the heat equation as This study is focused on deriving theoretical constraints on the timescales over which one may or may not expect to observe variability at a single spatially unresolved white dwarf. As a measure of the homogeneity at the white dwarf surface we define the contrast, Δ , between the origin ( = 0) and an arbitrary point on the surface ( = ) to be The circumferential positions which will be used to quantify the timescales of interest throughout this study will be the opposite side of the white dwarf ( = WD ) and a quarter of a circumference ( = WD /2). Fig. 2 shows the solution to Eq. (8) for the evolution of a delta function scalar on a circular ring of radius WD . This serves as a demonstration of the simple global transport model we employ to estimate the spreading timescales.
We now go on to present the diffusion coefficients derived from our grid of 3D models. We then implement these diffusion coefficients into the simple model outlined previously for simulating surface transport. This is followed by a presentation of the results from a model which accounts for the sinking of material and depth dependence of the diffusion coefficient, though for narrative clarity we introduce that second model as part of the discussion in Section 3.2.2.

Diffusion coefficients
The horizontal diffusion coefficient, surf , derived from tracerdensity evolution on nine simulations between 6000 and 15 000 K, are shown in Fig. 3. The tracer experiments were carried out with slabs of tracer density initially placed in the − and − planes, allowing to follow the diffusion process in the horizontal and directions, respectively. Each slab probes the entire vertical extent of the simulation box and does not allow for an explicit depth dependence of surf , instead providing a vertically-averaged horizontal diffusion coefficient. The final horizontal diffusion coefficient (and associated error) is computed as the mean (and standard deviation) of the horizontal diffusion coefficient derived for each slab.
The pink points in Fig. 3 represent the mean horizontal diffusion coefficient throughout the entire box. For the closed bottom simulations (≥ 12 000 K) the mean horizontal diffusion coefficient includes a contribution from the non-mixing layers beneath the convective and overshoot region. We account for this contribution by limiting the maximum amplitude of the tracer density distributions (see Fig. 4), effectively removing the vertical layers with no or very little mixing. The purple circles in Fig. 3 show the horizontal diffusion coefficient when only accounting for the regions where the tracer density is 10 per cent or lower than the maximum initial value. The arbitrary value of 10 per cent is used to crudely separate the parts of the tracer slab in the deep, non-diffusing layers from those in the shallow, rapidlydiffusing layers. A more robust, but computationally-demanding, approach would be to follow the trajectories of individual tracer particles, preserving the depth dependence of the horizontal diffusion coefficient. A major limitation with this approach however is the large range of advective time scales spanned by a single simulation box, where the vertical diffusion coefficient was found in CU19 to vary by up seven orders of magnitude. This crude clipping procedure is thus adopted to disentangle the contributions of horizontal diffusion at different depths, without doing the detailed path integration. The distributions resulting from a clipping value of 10 percent are shown in purple in Fig. 4.
To investigate whether this clipping procedure negatively impacts our derived diffusion coefficients we consider its behaviour on the open-bottom and closed-bottom simulations, respectively. The clipping procedure should not alter the open-bottom simulations because the entire simulation domain is contained within the fast-moving photospheric layers. The closed-bottom simulations are those which possess the deep, slowly-diffusing layers, the contribution of which we seek to remove from the derived horizontal diffusion coefficients. We see in Fig. 3 that the open-bottom simulations ( eff < 11 000 K) are barely affected by this clipping value of 10 per cent. On the contrary, the closed-bottom simulations ( eff > 11 000 K) exhibit an increase in the diffusion coefficient of 0.2-0.5 dex.
By extending this approach it is possible to provide an upper limit on the horizontal diffusion coefficient. We probe the most rapidly diffusing layers by considering only the lowest 0.01 per cent of each tracer density distribution, corresponding to the wings of the tracer densities (see Fig. 4). Given that convective velocities peak just below the photosphere (Tremblay et al. 2015a), it is likely that these are the layers being described by this diffusion coefficient (cyan on Fig. 3). In fact, this case is well represented by as shown by a dashed black curve in Fig. 3, where the characteristic length is half the characteristic granulation size, gran , defined from x p = πR WD x p = πR WD /2 Figure 2. Left: Evolution of a passive scalar, , according to Eq. 8 with a horizontal diffusion coefficient of surf = 2 × 10 9 cm 2 s −1 , characteristic of a white dwarf with eff = 12 000 K and log = 8.0, and a typical white dwarf radius of WD = 0.01288 . To make clear the evolution, every tenth (logarithmic) time step is indicated with a colour; from yellow to purple. Right: Examination of the time-dependent contrast between points of interest along the white dwarf circumference. Here we define contrast (Δ ) according to Eq. 13 and show the evolution of Δ at half and a quarter of the circumference. For the closed bottom simulations, the trace density in the stable lower layers serves to reduce the effective diffusion coefficient by averaging over layers outside of the convectively mixed region. To mitigate this effect we employ a clipping procedure at <100, <10 and <0.01 per cent of the initial maximum tracer density in pink (i.e., no clipping), purple and cyan, respectively (see Fig. 4). A lower clipping value corresponds to probing increasingly fast and rare local diffusion events within the resolved simulation. The cyan line, representing the fastest diffusing layers or resolved diffusion events, can be considered a robust upper limit on the horizontal diffusion coefficient for passive scalars. The purple line corresponds to our best estimate of a mean diffusion coefficient across the accessible convectively mixed layers. The black solid line shows the analytic approximation of the local diffusion coefficient given by Eq. (15). The black dashed line shows an alternative analytic prescription for the diffusion coefficient where the length scale is chosen to be the characteristic granule size, rather than the pressure scale height. the emergent intensity. The granule size corresponds most directly to the characteristic size of convective flows in the photosphere (Tremblay et al. 2013b), suggesting the equation characterises the mixing capabilities of the uppermost surface layers. We could not clearly establish the reason for the factor of 1 2 to the granulation size compared to the 1D diffusion equation of Eq. (6), but this could be because granulation size is a 2D representation of a 3D process, or that in the definition of granule size a single granule is composed of both a downdraft and an updraft.  To remove the contribution of the (non-diffusing) stable lower layers from the mean tracer densities, the distributions are clipped at <10% (purple) and <0.01% (cyan) of the initial maximum density as described in the text. These colours correspond to the lines shown in Fig. 3. The boundary conditions in the x-dimension are parallel such that the tracer density ( tr ) satisfies tr (0) = tr ( b ), for b = 7.5 km. The apparent gap at = 3.5 km is a consequence of plotting only a subset of the 28 tracer density slabs included in this simulation. This omission is for clarity only, and the horizontal diffusion coefficient of each slab is not expected to vary significantly.
The photosphere only contains a small amount of mass compared to the full convection zone, hence Eq. (14) is unlikely to be a good representation of the overall horizontal mixing of accreted debris in the convection zone. To represent this process, we use instead, as describe above, the clipped estimate where the tracer density is 10 per cent or lower than the maximum initial value. By construction of our models (Tremblay et al. 2013a), this corresponds to mixing in surface layers above Rosseland optical depth R 1000. This case is best described by which is shown by a solid black curve in Fig. 3. The characteristic length real P is the actual pressure scale height that corresponds to the geometrical distance between the photosphere ( R = 1) and the layer below where the pressure drops by a factor . It is not clear whether the factor of real P /2 in Eq. (15), when compared to Eq. (6), has a fundamental physical justification, as here the diffusion process corresponds to an average over a range of convective layers.
For open bottom simulations with eff < 12 000 K, Eq. (15) may still represent an upper limit on the efficiency of the overall diffusion of accreted material across the stellar surface. The reason is that the lower portions of the convection zones, which are not included in open bottom simulations, have lower convective velocities and are possibly less efficient to mix material. While this would not necessarily delay the process of achieving a fully homogeneous surface composition across the stellar disc, the actual absolute metal abundances could vary on longer timescales, as the bottom of the convective zone is still adjusting. A lower limit on horizontal diffusion could possibly rely on using the pressure scale height at the bottom of the convection zone in Eq. (6).
We provide the results of the tracer-density evolution experiments in Table 1 where one can find the mean horizontal diffusion coefficient in the mixed layers (purple circles in Fig. 3) as log surf . We also provide the upper limits (corresponding to the photosphere) on the diffusion coefficient at each effective temperature as log max surf . To test the convergence of our results we also performed a direct tracer experiment on a simulation with twice the spatial extent in the x and y dimensions. From CU19 we adapted the simulation C1-2 which was identical in setup to simulation C1 from the same study -with eff = 13, 500 K, log = 8.0 and box size of 150 3 -except that the x and y dimensions each extended to 15 km, rather than 7.5 km. The results of direct tracer experiments for this simulation are shown in Fig. 3 with open triangles. We find the results of this simulation (C1-2) to be in reasonable agreement with the one with a lesser spatial extent (C1). The minimum diffusion coefficient (pink), which includes the slowest diffusing layers, is in excellent agreement, whilst the most physical estimate for the diffusion coefficient ( surf ; purple), after the removal of the contribution from the overshoot layers, agrees to within a standard deviation.
By way of comparison to previous studies, Montgomery et al. (2008) derived a horizontal diffusion coefficient of 1.5 × 10 10 cm 2 s −1 for the metal-polluted white dwarf G29-38 using white dwarf evolutionary models and stellar parameters constrained by observations. The derived parameters for this star are eff = 11 400 K and log = 8.02 (Gentile Fusillo et al. 2019) and using these we can compare the horizontal diffusion coefficient of Montgomery et al. (2008) to those shown in Fig. 3. We see that it is in reasonable agreement with our predicted maximum value max surf (blue) corresponding to the photosphere. However, it is about an order of magnitude larger than our adopted coefficient (purple) considering mixing over the entire convection zone. Table 2 contains the diffusion coefficients derived analytically using Eqs. (14) and (15) from mean 3D properties for the grid of the warmer (>11 400 K) convective pure-hydrogen atmospheres at log = 8.0 presented in CU19. These quantities correspond to the black solid and dashed line shown in Fig. 3. In the same table and Fig. 5 we show the same quantities for the grid of cooler (<13 000 K) convective DAs from Tremblay et al. (2013a) for log = 7.0-9.0. The two grids have an overlap in effective temperature between 11 400-13 000 K at log = 8.0. In this region we include only the values for the warmer grid (i.e., that from CU19). We find the diffusion coefficients where the grids overlap agree to within 0.06 dex, well within the estimated uncertainty of 0.1 dex. We also present the diffusion coefficients across the grid of helium-rich models with log (H/He) = −10.0 published by Cukanovaite et al. (2019). These are given in Table 3 and Fig. 6. Table 1. Diffusion coefficients derived from direct diffusion experiments using nine simulations which were adapted from CU19 and Tremblay et al. (2013a) for pure-hydrogen atmosphere white dwarfs at log = 8.0. The horizontal diffusion coefficient is given as surf which represents the mean diffusion coefficient across the convectively mixed layers and is shown explicitly in Eq. (5). The associated error is the standard error on the mean. The upper limits on the diffusion coefficient (corresponding to the photospheric layers) are also given. The simulations used to derive these results included a surface gravity of log = 8.0, grey opacities and boxes with 150 grid points in all three (x, y, and z) dimensions.  Figure 5. Horizontal diffusion coefficients for pure-hydrogen white dwarf models defined using mean 3D values derived from the grid of models presented in Tremblay et al. (2013a). Diffusion coefficients are derived using Eq. (15) (circles, solid lines) and Eq. (14) (plusses, dotted lines). The latter represents photospheric diffusion which serves as an upper limit on the horizontal diffusion coefficient, max surf . The mean error on the diffusion coefficients derived from direct diffusion experiments (Table 1) was found to be (log surf ) = 0.1 which we adopt as an estimate for the error on points shown in this figure. The surface gravity (log ) of each line pair is shown in the corresponding colour. All values shown here are displayed in Table 2.

Surface spreading
In the following we consider two idealised geometries of accretion. Firstly we examine instantaneous accretion onto a point source. In nature this could arise from, for instance, the direct impact of an asteroid. Probabilistic arguments suggest that this scenario could provide on the order of 1-2 per cent of metals in white dwarfs. Nonetheless it allows for a highly idealised problem setup and bears Table 2. Diffusion coefficients defined by Eq. (15) for pure-hydrogen DA white dwarf atmospheres. Computed using mean 3D quantities from a grid of cool DAs ( 13 000 K) from Tremblay et al. (2013a), except for those with 11 600 K eff 18 000 K and log = 8.0, which were computed using the grid of deep simulations from CU19. The simulations leading to the results in this table utilised a non-grey opacity scheme for the radiative transfer. We also show the upper limit on the diffusion coefficient (corresponding to the photospheric layers) defined by Eq. (14). The real effective temperature of a simulation can stabilise at a slightly different effective temperature to the desired value. To represent this we include the quantity Δ( ) such that the true effective temperature is given by eff + Δ( ). An explicit error is not available for these results, but we note that typical statistical errors from the direct diffusion experiments (Table 1) were on the order of (log surf ) ≈ 0.1. An additional source of uncertainty arises in Fig. 3 from the imperfect fit of the directly derived diffusion coefficients to those computed from physical parameters. We find at the warmest end of the grid the difference is as large as 0.2 and 0.4 dex for surf and max surf , respectively. relevance for the decreasing phase at the end of any steady-state accretion.
Secondly, we consider the more likely scenario of continuous accretion onto a point. This would correspond to the canonical model of debris from a disc constantly feeding the white dwarf surface. Both models can be scaled to account for equatorial accretion, rather than accretion onto a spot. If the white dwarf possesses a sufficiently strong magnetic field, material may be channelled toward the poles (Metzger et al. 2012). This scenario is well described by our second model, where material is constantly fed to two spots, rather than one. We revisit the importance of magnetic fields regarding the delivery of material to the surface in Section 4.3 and the impact of magnetic fields on convective instabilities in Section 4.1.1. In the following we focus our modelling on the transport of material due to convective motions after it arrives at the surface.

Instantaneous accretion
We begin by considering the evolution of a system whereby a large quantity of a passive scalar is instantaneously delivered either to a point on the white dwarf surface or equatorially, e.g. from an accretion disc.
In the first model, to account for the large range of timescales involved, we consider the spreading and sinking to be governed by two independent, one-dimensional equations. One governs the horizontal spreading of material and assumes the passive scalar concentration is restricted to the convection zone. The other equation relaxes the massless assumption, allowing the accreted material to sink via microscopic diffusion processes; i.e., gravitational settling. Clearly this problem setup is a strong simplification but it allows for a simple analytic formulation of the problem. The analytic approach refers to the problem setup outlined in Section 2.2.  Fig. 7 shows the time required, cover , for the metal abundance at the white dwarf surface to homogenize. From black to grey this is shown for increasing levels of homogeneity; namely contrasts of Δ = 0.001, 0.1, 0.37, 0.99. The vertical diffusion (sinking) times from Koester et al. (2020) are shown in dotted purple for eight elements across the full effective temperature range. The mean sinking time arising from the model of convective overshoot (CU19) is shown in green for effective temperatures between 11 400-18 000 K. For comparison our updated estimate of a typical disc lifetime is shown in dashed blue, with the associated error (see Section 4.1.2). We  Table 3.
find moderate dependence on whether the accretion is spot-like or equatorial from this model. The expected behaviour that equatorial accretion homogenizes faster is recovered, but the difference is no more than a factor of three. We introduce in Fig. 8 the ratio of the spreading time, cover , and the sinking time, sink . This dimensionless quantity is a metric by which we will assess the capability of a white dwarf envelope to efficiently spread material. One can consider that when this ratio is below unity the surface layers are able to spread material out quicker than it sinks. In this figure we show the cover / sink ratio for cover corresponding to Δ = 1/ ≈ 0.37, i.e., chemical abundance variations of less than a factor ≈3.

Continuous Accretion
So far we have considered the scenario that the accreted metals arrive in a small region on the surface at one instant in time. This would correspond to a direct impact type scenario and is considered one of the least likely sources of accretion at white dwarfs (Wyatt et al. 2014). We have also required that the metals remain confined to the surface when modelling the spreading. We now relax these criteria and go on to consider the scenario in which a white dwarf accretes for a prolonged period with constant accretion rate. This would represent the stable accretion of material from a debris disc.
To model the behaviour of the metals at the surface of the white dwarf in such a system we adopt a two-zone model. In this model the horizontal diffusion coefficient, surf , is confined to some upper region and a lower adjacent region has a drift velocity in the vertical direction, drift . This describes the physical scenario of a trace concentration being removed from a mixed region at the surface with a drift velocity in the vertical direction defined at the base of said region.
The equation describing the evolution of this scenario is as follows where = ( , , , ) is the concentration of a passive scalar de- Koester (2020) 10 8 10 9 Cooling Age [yrs] The overshoot correction in the range 18 000-11 400 K is taken from (CU19) (thick) and is computed as the mean correction over elements Fe, Ca and Mg. The overshoot correction for cooler models is estimated to be an additional one pressure scale height of mixing below the lower Schwarszchild boundary (thin). The dashed cyan line depicts an estimate of a disc lifetime and the expected uncertainty (see Table 4). This figure highlights the strong, effective temperature dependence of the characteristic time scales of horizontal spreading and vertical sinking. The minimum in cover around 12 000 K corresponds to the maximum in the horizontal convective velocity at that temperature (see pending on depth, surface position and time. Once this system has reached a steady state, the equation reduces to For a constant flux delivered to = = 0, this equation has the analytic solution of the form the base of the surface zone and is the horizontal extent of the surface. We are interested in the global spreading properties which for a spot would require = = WD and for an equatorial belt would require = = WD /2. Thus in the case of a white dwarf with equatorial accretion, for convection to be efficient at spreading before material sinks requires where cvz is the depth of the convection zone, drift represents a characteristic velocity for the settling of a trace amount of metal out of the convection zone. This can also be written where the transport timescales are estimated as For the drift velocity we adopt the vertical diffusion velocity at the base of the mixed region. Fig. 9 shows Cooling Age [yrs] Figure 10. Comparison of two models, instantaneous accretion (purple lines) and continuous accretion (grey circles), in the comparison of spreading and sinking timescales for H-rich atmospheres at log = 8.0. The instantaneous accretion model is shown where the vertical sinking velocity, drift , is either given by the 1D mixing length theory results of Koester et al. (2020) or from the 3D overshoot results of CU19 (see Fig. 9) in dotted or solid, respectively. The spreading timescale is defined for a maximum surface abundance contrast of Δ = 0.37 from equatorial accretion (reproduced from lower panel; Fig. 8 (pink-dotted) and at the base of the overshoot region (purple-solid) taken from the diffusion calculations of Koester et al. (2020) and overshoot determination of CU19, respectively. As was shown by the latter, convective instabilities arise at 18 000 K. In the radiative regime ( eff > 18 000 K) diffusion velocities are given at an optical depth of R ∼ 1. We caution that none of these models include the effects of radiative levitation and thermohaline mixing which may become important in the radiative regime (Chayer et al. 1995a,b;Bauer & Bildsten 2018). We plot the left-hand side of Eq. (19) in Fig. 10 in grey circles and we include the results with (solid) and without (open) overshoot. We also reproduce the results for instantaneous accretion (see Fig. 8) in purple, again with (solid) and without (dotted) overshoot. We find the two models make almost entirely overlapping predictions on the ability of a white dwarf to homogenize its surface within a diffusion timescale.
The models we have presented highlight the horizontal mixing capabilities across the full range of convective H-atmosphere white dwarfs. In Figs. 3 & 5 and Table 2 we provide the horizontal diffusion coefficients for all convective pure-H atmosphere white dwarfs with surface gravities log = 7.0-9.0. For brevity, in this study we have only compared those at log = 8.0 with vertical diffusion coefficients (Figs. 7,8 & 10), as the results are similar for other surface gravities but with a small shift in effective temperature. The effect of trace metals in H-rich model atmosphere calculations is very small or negligible, hence our results directly apply to DAZ white dwarfs. We now go on to consider the surface spreading behaviour in heliumatmosphere white dwarfs.

Surface spreading in He-atmosphere white dwarfs
Approximately half of the known polluted white dwarfs have heliumrich atmospheres (Farihi et al. 2016). In this section we repeat some of the procedures detailed in the previous sections to provide a picture of the mixing capabilities of He-atmosphere white dwarfs as they cool. The lower panel of Fig. 11 shows the ratio of the spreading and sinking times for our 3D grids of pure-He (pink) and pure-H (green) model atmospheres. A 3D correction for convective overshoot is not readily accessible for the DBs -requiring an overshoot study similar to CU19 for the early convective DBs. This is outside the scope of this work, where we instead adopt an arbitrary overshoot correction of one pressure scale height of additional mixing below the lower Schwarzschild boundary. A comparison of the solid and dotted pink lines in Fig. 11 shows that this assumption does not significantly impact the results of the pure-He models. Given that the majority of metal-polluted DBs have effective temperatures less than 25 000 K, this figure suggests that metal-polluted DBs are far more likely to be homogeneously mixed within a diffusion timescale. While there are no 3D simulations of helium-rich atmospheres cooler than 12 000 K (DC, DZ and DQ white dwarfs), diffusion timescales calculated from 1D model atmospheres suggest that spread / sink plateaus at values in the range 10 −2 -10 −3 , depending on trace hydrogen and metal abundances.

APPLICATION TO WHITE DWARF POPULATION
For any white dwarf, the numerical solution of Eq. (16), using the horizontal diffusion coefficients presented in this paper as well as vertical diffusion coefficients from CU19 or the literature, would yield a reasonable description of the distribution of metals across the white dwarf surface. This could provide the basis for the generation of synthetic spectra under different accretion conditions (geometries and time-dependence) and rotational periods. In Section 5.1.1 we generate synthetic spectra for the DAZ white dwarf SDSS J104341.53+085558.2, assuming the observed metals are confined to varying spot sizes on the visible surface. In the following, we firstly explore a heuristic discussion of these systems across a wide parameter space. The results we have shown so far focus mostly on the competition between horizontal spreading and vertical sinking. This led to a separation, based on atmospheric parameters, between metal-rich white dwarfs which are expected to have a homogeneous surface composition and those where metals must be concentrated in and near regions of accretion. When it comes to predict how likely it is to observe photometric or spectral variability in a white dwarf, we also need to consider the amplitude, periodicity and total duration of a predicted contrast, which depend on further physical variables. In addition to (1) the white dwarf atmospheric parameters (spreading and sinking times), we must also consider (2) the geometry of accretion, (3) the instantaneous accretion rate, (4) the duration of the event (e.g. the disc lifetime combined with the sinking time), and (5) the white dwarf rotation period and inclination with respect to the line of sight. Fortunately, some of these quantities are observationally or theoretically constrained and limiting cases are sometimes sufficient when considering the white dwarf population as a whole.

Accretion geometry
To predict the time evolution of surface inhomogeneities on white dwarfs, one should consider the initial distribution of metals onto the white dwarf surface, whether that be a spot or multiple spots as seen in the Shoemaker-Levy 9 disruption and accretion of an asteroid onto Jupiter (Brown et al. 2017), a stable accretion disc funneling material either to an equatorial belt or to the poles (Wyatt et al. 2014;Farihi et al. 2018), or a near-homogeneous (spherical) deposition. We emphasise that we obtain a trivial result if accretion is fully spherical, in which case the surface is homogeneous and diffusion is only of relevance in the radial direction, where 3D convective effects nevertheless still play a role. We will consider the relevance and motivation for each of these scenarios in the following.
Direct impacts: To understand how prevalent direct impacts are as a source of observed metal pollution one must consider the likelihood of an asteroid colliding directly with the white dwarf surface. A key input to this problem is whether an asteroid can survive the tidal forces as it approaches the white dwarf. This can be quantified with the Roche radius which describes the distance at which a substellar body, held together by its own self-gravity, will succumb to the tidal forces of the white dwarf and tidally disrupt.
where is the density of the substellar body and = 1.3-2.9 accounts for additional variables not explicitly included here such as strength and shape. If a given substellar body is assumed to be a strengthless rubble pile, the Roche radius is typically Roche ∼ 100 WD . This implies that direct impacts are unlikely to contribute more than 1-2 per cent of metal pollution at white dwarfs (Wyatt et al. 2014). If the substellar body has internal strength, however, the Roche radius is decreased (Veras 2016) and such high-strength objects can survive many passes within the rubble pile Roche radius of the white dwarf. For instance, the planetesimal observed at SDSS J1228+1040 ) lies at 0.73 , which is 25 per cent within what is typically considered the white dwarf's rubble pile Roche radius. The existence of such a system provides evidence that solid bodies can exist within the Roche radius, and thus direct impacts cannot be ruled out as unphysical. However, the frequency of metal-pollution in white dwarfs with short sinking times (i.e., hydrogen-rich DAZ white dwarfs) strongly supports a model of (continuous) accretion from a long-lasting source, such as an accretion disc.
Accretion discs: If a substellar body comes within a Roche radius and does not result in a direct impact it most likely suffered tidal disruption, forming a disc of dust and gas. There is significant evidence that the bulk of metal pollution at white dwarfs arrives this way. For one, warm DAZs have sinking time scales on the order of days to years, so the likelihood of observing them is low, unless the source of material has longevity. Furthermore, the observed frequency of IR excess due to dust at metal-polluted white dwarfs confirms the importance of this scenario (Rocchetto et al. 2015). Once a disc is formed it is likely to accrete (quasi-)continuously driven by drag forces (Rafikov 2011;Metzger et al. 2012). This provides the ongoing source of material at white dwarfs, but constraints on where the material arrives on the stellar surface are inexistent.
It has been shown that magnetic fields are likely to play an important role in the trajectory of material between the disc and the white dwarf surface. The Alfvén radius describes the distance at which the magnetospheric and disc stresses equalize, allowing a surface magnetic field to dynamically alter the trajectory of in-falling material. If the surface magnetic field is dipolar, the Alfvén radius for a gas disc accreting onto a white dwarf can be written (Ghosh & Lamb 1978) as where WD is the surface field strength and the accretion rate. For a white dwarf with mass WD = 0.6 and radius WD = 9.0 × 10 8 cm this becomes If the Alfvén radius is smaller than the white dwarf radius ( A < WD ) then accretion from a disc is most likely to land in an equatorial belt on the white dwarf surface. On the other hand, if A > WD , then any material within the Alfvén radius, and external to the white dwarf, is likely to be magnetically driven, i.e., towards the magnetic poles. Metzger et al. (2012) showed that fields as weak as 0.1-1 kG were sufficient to affect the flow of gas near the sublimation radius around white dwarfs. The authors pointed out two caveats to this prediction. Firstly, current observations constrain the total surface magnetic field, of which the dipole may only be a minor component. And secondly, debris in actively accreting systems could be diamagnetically screened, essentially exposed to systematically lower magnetic fields.
Another geometric possibility is accretion through multiple misaligned discs/rings, only a fraction of which are within our line-ofsight. Asteroids scattered to the Roche radius are highly unlikely to all be in the same plane. In fact, they could easily vary in inclination. Mustill et al. (2018) showed that the orbital inclinations of bodies as they cross the Roche limit are roughly isotropic. Whilst we have only modest constraints on the accretion geometry at metal-polluted white dwarfs, a comparison can be drawn with a better studied class of accreting white dwarf system; cataclysmic variables (CVs).

Analogy to accretion on CVs:
The hypothesis of magnetically-driven accretion is observationally confirmed for the highly-magnetic CVs (polars) such as AM Herculis (Gänsicke et al. 2001;Schwope et al. 2020). But although they provide the closest, well-studied analogy to debris-accreting white dwarfs, it is important to understand that the range of accretion rates of 10 7 − 10 11 g s −1 estimated for metalpolluted white dwarfs (Fig. 12) is 5-10 orders of magnitude lower than that observed in CVs.
In the equatorial accretion belt scenario, material is slowed down from Keplerian velocities in the inner disc to the rotation velocity of the white dwarf. Since the free-fall velocity for a white dwarf is ∼1000 km s −1 , compared to rotational velocities of the order of 1-10 km s −1 , that must involve significant shear, turbulence, shocks and mixing above the white dwarf atmosphere. For CVs, Patterson & Raymond (1985) found that if the accretion rate is sufficiently high ( > 10 16 g s −1 ) the accretion disc may survive all the way down to the white dwarf surface. For lower accretion rates ( < 10 16 g s −1 ), still orders of magnitude above the range observed for metal-polluted white dwarfs, they found that the disc was likely to evaporate before reaching the white dwarf surface, forming a more diffusive equatorial belt, or even a spherical shroud, covering much larger latitudes. A theoretical study by Meyer & Meyer-Hofmeister (1994) considered the accretion disc around dwarf novae. They found that a siphon flow is capable of evaporating some of the inner disc, leading to a hole between the white dwarf surface and the inner disc. These studies support the idea of close-to-spherical accretion for CVs with low accretion rates ( < 10 16 g s −1 ). However, Piro & Bildsten (2004) found the mixing region to be relatively small and restricted to less than one degree of the equator.

Modelling disc lifetimes
With enhanced mixing from convective overshoot (CU19), our new diffusion models can be used to revisit and update a known issue with metal accretion rates and duty cycles across the population of all polluted white dwarfs, especially the distinct behaviour of Hversus He-rich atmospheres. It has been known for some time that, all else being equal, if identical calculations are made for instantaneous accretion rates onto known DAZ stars and known DBZ or similar He-rich stars, there is a sizable gap between the two inferred rates for each sub-population ). Based on the new results presented in this work and in CU19, these accretion rate inferences, and their implications are updated.
H-atmosphere white dwarfs permit instantaneous accretion rate inferences based on their relatively short diffusion timescales, which reduce the full evolutionary accretion behaviour to a simple equation (Jura et al. 2009, Equation 1). This is not the case for He-dominated atmospheres because their sinking timescales are sufficiently long that the full differential equation is not immediately solvable without a knowledge of the full history. Nevertheless the same, simple equation can be applied to better understand the population as a whole, and to put both H-and He-atmosphere stars onto the same footing . Previous work has shown that when this is done carefully with the most up-to-date models at the time, and assuming the same fractional representation of metal abundance by primary atmospheric constituent, there is a notable difference where some of the He atmosphere stars exhibit accretion rates that are orders of magnitude larger than the largest H-atmosphere accretion rates .
Assuming that the accretion rate is constant and that the accretion and diffusion (atmospheric settling) processes have reached a steady state, the rate of accretion can be estimated following Koester (2009) where is the atmospheric abundance of element , cvz is the mass of the convectively mixed surface layers and is the diffusion timescale at the base of the mixed region. It must be kept in mind that these are not true accretion rates, but an average over the past sinking timescale, where the calculation simply takes the mass of metals currently in the mixing layer and divides by a characteristic sinking timescale. Nevertheless, they are informative in at least two ways. The first is in the numerator of this calculation, which represents the total mass of metals currently within the outer, visible and well-mixed layers of the star. The difference between the H-and He-atmosphere stars indicates that metals continue to accumulate far longer than a typical sinking timescale for H-atmosphere white dwarfs, and this difference has been used to calculate a typical disc lifetime . Specifically, if one takes the larger mass of metals that can be inferred within He-atmosphere stars, and divides by the instantaneous accretion rate inferences for H-atmosphere stars, the outcome has dimensions of time, which can be interpreted as a disc lifetime, if the accretion rate is constant over a typical He-atmosphere sinking timescale.
The second difference is in the accretion rates themselves. When the ongoing accretion rates for H-atmosphere stars are compared with the time-averaged accretion rates for He-atmosphere stars, there is a similar, stark difference, and the only hypothesis in the literature suggests this is likely due to highly variable accretion rates ). The idea is straightforward: only modest accretion rates are seen to be ongoing currently via H-rich white dwarfs, but over a typical, Myr timescale for metal diffusion in He-rich white dwarfs, there have been high-rate bursts that are not yet witnessed as ongoing. This result implies the empirical disc lifetime estimates made from these calculations are uncertain at best, but it does provide a constraint on the duty cycle of high-and modest-rate accretion for the population of white dwarfs as a whole, from which further tests might be made.
In Fig. 12, the new model calculations provide updated accretion rates and accreted masses for white dwarfs with confirmed IR excess. It is most noteworthy that they corroborate the previous findings discussed above. That is, there remains a gap between the time-averaged accretion rates for He-atmosphere stars, and the instantaneous accretion rates for H-atmosphere stars. Table4 provides the mean accretion rates for DAZ stars and accreted masses for DBZ stars with confirmed IR excess, values which we will use in the following discussion regarding disc lifetimes. Girven et al. (2012) derived a crude estimate for the lifetime of a disc of disrupted planetesimals around a white dwarf to be disc ∼ DBZ/DZ,disc where DAZ,disc is the mean inferred accretion rate at warm polluted hydrogen-rich DAZ stars and DBZ/DZ,disc is the mean mass of metals enclosed in He-rich (DBZ and DZ) white dwarf convection zones, with both quantities derived from a sample of white dwarfs with confirmed IR excess. As discussed before, this calculation assumed a Jura-like disc model and a constant accretion rate over the full sinking timescale. In these calculations it is estimated that Ca makes up 1.6 per cent of the accreted material, a value characteristic of the bulk Earth abundance. Girven et al. (2012) found for the thirteen DAZ white dwarfs with IR excess confirmed with Spitzer the mean accretion rate to be DAZ,disc = 9.7 × 10 8 g s −1 . The average total accreted mass for the eight DBZ/DZ stars with confirmed IR excess was found to be DBZ/DZ,disc = 4.1 × 10 22 g. The accretion rates and accreted masses span orders of magnitude in range (see Fig. 12), and thus the authors favoured a logarithmic mean, such that the mean accretion rate in the DAZ stars was found to be log DAZ,disc = 8.8 [g s −1 ] and mean total accreted mass in the DBZ/DZ stars log DBZ/DZ,disc = 21.9 [g]. This led to a disc lifetime estimate of log( disc ) = 5.6 ± 1.1 [yr].
We revisit this disc lifetime estimate using the same set of assumptions as in the Girven et al. (2012) study and repeat the same analysis with the now slightly larger sample of metal-polluted white dwarfs which exhibit IR excess (Farihi et al. 2016). This comprises nineteen hydrogen-rich DAZ stars and eleven DBZ/DZ stars, with the accreted masses and accretion rates shown in Fig. 12. We consider two definitions of the DAZ accretion rates. One with the base of the mixed region defined at the lower Schwarzschild boundary, akin to standard 1D accretion-diffusion calculations. Secondly, we include the overshoot correction of CU19.
In the no-overshoot picture we find the mean accretion rate in the DAZ stars to be log DAZ,disc = 8.6 [g s −1 ] and mean total accreted mass in the DBZ/DZ stars log DBZ/DZ,disc = 22.4 [g]. This gives rise to a disc lifetime estimate of log( disc ) = 6.3±1.4 [yr]. We find that the new observed sample with the updated sinking times from Koester et al. (2014) supports a longer disc lifetime (Myr) than the earlier estimate of Girven et al. (2012). Without overshoot the disc lifetime estimate has increased by almost an order of magnitude, though due to the large uncertainties in this calculation the two estimates are still consistent.
The inclusion of convective overshoot in the analysis yields a disc lifetime estimate which is just 0.2 dex shorter than its non-overshoot counterpart. This is a direct consequence of the fact that including an overshoot correction to the DAZ stars increases their mean accretion rate by 0.2 dex to log DAZ,disc = 8.8 [g s −1 ]. This relatively small 3D correction is in part because a third of the DAZ stars in the sample have effective temperatures above 18 000 K for which there is no overshoot correction. For white dwarfs with eff < 11 400 K we use the correction applied to the coolest end of the overshoot grid. The mean total accreted mass of the DBZ/DZ stars is also expected to increase with the inclusion of convective overshoot but detailed modelling of this mixing process in helium-atmosphere white dwarfs is beyond the scope of this work. We therefore propose a disc lifetime estimate of log( disc ) = 6.1 ± 1.4 [yr].
Whether or not the accretion-diffusion calculations include convective overshoot, both sets of models support the hypothesis of highly variable accretion rates over the Myr timescales that are typical for metal diffusion in He-atmosphere stars. Although the difference is now smaller with the CU19 overshoot models, the likelihood of catching a high-rate burst in the future, in real time, may be possible with an increasingly large sample of DAZ spectra via large and ongoing surveys. Furthermore, these high-rate bursts of accretion should be accompanied by higher rates of X-ray production and might be detectable ) by surveys such as eRosita.

Predicting variability
Identifying the parameter space of metal-polluted white dwarfs most likely to exhibit spectral variability would be of great advantage for observational campaigns. The parameter space is defined by many independent variables including, but not limited to: vertical sinking time scale, horizontal spreading time scale, disc lifetime and accretion geometry. Given the number of independent variables, we split our discussion based on white dwarf temperature and spectral type, which are readily available observational quantities.

Cool H-atmosphere white dwarfs
Below 10 000 K, H-atmosphere white dwarfs have relatively large convection zones and long sinking timescales. For instance metals will stay visible for at least 10 kyr at 7000 K, but it could be 2-3 orders of magnitude longer when considering the disc lifetime and phase of decreasing abundance. Fig. 10 shows that convection is efficient at spreading metals across the surface in this regime, under any accretion geometry or whether it is instaneous or continuous. Indeed, the time to homogenize the surface such that global variations are less Table 4. Mean parameters for white dwarfs with inferred accretion rates and confirmed IR excess. Accretion rates and accreted masses without overshoot were taken from Farihi et al. (2016) for a sample of 19 H-rich (DAZ) and 11 He-rich (DBZ/DZ) metal-polluted white dwarfs with a detected IR excess in Spitzer observations. For the the DAZ white dwarfs, these quantities were also updated with the overshoot correction from CU19. The accretion rates and accreted masses used to derive the values in this table, with and without the overshoot correction, are shown in Fig. 12. We include the inferred logarithmic mean DAZ accretion rate, log DAZ,disc , and the logarithmic mean DBZ/DZ accreted mass, log DBZ/DZ,disc . The number of objects in the DAZ and DBZ/DZ samples is given by DAZ and DBZ , respectively. We also indicate whether the modelling included convective overshoot.  than 50% takes only ∼1000 yr. Thus surface variations could only be observed during the increasing, non-steady phase of accretion, which is less than 0.1 percent of the disc lifetime. Thus the probability of finding a heterogeneous surface is very low.

Warm H-atmosphere white dwarfs
Above 13 000 K, a DA white dwarf will have a convection zone which is relatively small in radius and mass (log cvz / WD −14) with maximum convective velocities on the order of 1 km s −1 . Fig. 10 shows that such a star is unable to effectively homogenize its surface, regardless of whether accretion is instantaneous or continuous over a disc lifetime. The reason for this behaviour is that any metals accreted do not have enough time to travel across the stellar surface before they diffuse out of the convection zone.
For temperatures above 18 000 K, no convective instabilities are expected to develop in DA white dwarf atmospheres, i.e. those are socalled radiative atmospheres. Thus our improved convective model is not directly relevant. In this regime thermohaline instabilities may be the dominant form of surface mixing (Bauer & Bildsten 2019). These instabilities have similar sizes and timescales compared to convection in the vertical direction, but it is unlikely they contribute significant spreading in the horizontal direction (Garaud 2018), hence we expect spreading times in warm radiative white dwarfs to be similar to those of warm convective objects. As a consequence, it is unlikely that thermohaline mixing can lead to significantly more homogeneous stellar surfaces than predicted in this work.
If the accretion rate is variable in time, additional variations are expected over the sinking time, which -in the case of white dwarfs above 13 000 K -is of the order of days to months (see Fig. 7). Radiative levitation must also be considered when predicting such temporal abundance variations. Overall, we predict that all DAZ white dwarfs above 13 000 K are expected to show surface inhomogeneity. The only known physical mechanism that would be able to prevent this outcome is near-spherical accretion geometry.

H-atmospheres at intermediate temperatures
We see in Fig. 10 that the inclusion of convective overshoot makes a considerable impact to the homogenizing efficacy of the surface layers. Because overshoot has no impact on the horizontal diffusion coefficient one can explain this as a direct consequence of the longer sinking times that result from including overshoot in the modelling (because the base of the mixed region is in deeper, slower diffusing layers). Furthermore, solving Eq. (16) for continuous accretion should in principle include the depth dependence of the horizontal diffusion coefficient (Fig. 5), which is likely to result in a relatively complex surface abundance profile with extended low abundance tails corresponding to the faster photosphere spreading. Therefore, it is difficult to predict a firm outcome for hydrogen-dominated DA white dwarfs in the range 10 000-13 000 K.
We note that this temperature range corresponds to the quickest spreading time, with variations expected to smear out within 100 yr (see Fig. 7). This is a direct consequence of the diffusion coefficient experiencing a maximum at this temperature (i.e., the log = 8.0 track of Fig. 5), which is directly tied to this being the peak in convective velocities. However, this is somewhat balanced out by the relatively short sinking timescales on the order of 1-100 yr.

Helium-rich atmospheres
Convective instabilities in white dwarf surface layers are driven by the opacity of the dominant atmospheric constituents. For DAs the hydrogen opacity is maximal around 12 000 K, whereas for DBs the helium opacity peaks around 28 000 K. One might thus expect a similarity between the two evolutionary tracks, offset by ≈16 000 K. An examination of Fig. 11 reveals such behaviour. We see that DB white dwarfs with effective temperatures cooler than ≈30 000 K are very efficient at moving accreted metals across the surface and that they are likely to be fully mixed within one diffusion timescale, similarly to cool H-atmosphere white dwarfs below 10 000 K as discussed above. Since very few DBZ white dwarfs are known to have temperatures above ≈30 000 K (Kepler et al. 2019), this scenario applies for the vast majority of known DBZ and DZ white dwarfs.

Magnetic white dwarfs
We have discussed previously that global magnetic fields may influence the accretion of planetary material from a debris disc at metalpolluted white dwarfs. We now consider the behaviour of magnetic fields in relation to the local white dwarf atmospheric layers. The impact of magnetic fields on the stellar structure can be estimated from the plasma-parameter where is the thermal pressure and WD is the average magnetic field strength at the surface. It was found from 3D magneto-hydrodynamic (MHD) simulations (Tremblay et al. 2015b) and observations (Gentile Fusillo et al. 2018) that a value of < 1 will inhibit convective energy transfer in the white dwarf atmosphere. Setting plasma-= 1 can therefore provide an indication of the critical magnetic field required to heavily damp convective flux in the atmosphere. Fig. 13 shows the critical magnetic field for white dwarfs with a surface gravity of log = 8.0. Magnetic fields larger than this critical value, corresponding to 2-100 kG in DA and DB white dwarfs, and 0.1-1 MG in DC/DZ stars, will therefore inhibit convective energy transfer in the atmosphere. As a consequence both surface spreading and sinking of metals are likely to have different timescales. We emphasise that the lack of significant convective energy transfer does not necessarily imply that mixing is inhibited. In fact, Tremblay et al. (2015b) demonstrated that for a pure-H atmosphere white dwarf at 10 000 K with an average field of 5 kG, convective velocities were of similar amplitude to the non-magnetic case even though convective energy transfer was largely suppressed. Furthermore, Alfvén waves may contribute to additional mixing processes. In many of the limiting cases we have discussed so far, such as cool DAZ and DZ white dwarfs, changing the diffusion coefficients by even one order of magnitude would not necessarily make surface spreading slower than sinking times, nor result in any inhomogeneity. Magnetic fields are also likely to influence debris disc formation and geometry of accretion. It has been shown that the Alfvén radius, , of a white dwarf with surface field strengths of > 1 kG lies near or outside the typical Roche limit for accretion rates of < 10 10 g s −1 (Farihi et al. 2018). In such a configuration it is highly likely that accreted gas is redirected along magnetic field lines, toward the poles. It was also shown that magnetic fields as weak as = 1 G have an Alfvén radius larger than the white dwarf radius for accretion rates up to < 10 16 g s −1 . There is however no observational evidence of white dwarfs with magnetic fields weaker than ≈ 1 kG, or any prediction that global magnetic fields in the range 1-1000 G could remain stable throughout white dwarf evolution including after the onset of convection.
Given that magnetic fields can affect both accretion flows and atmospheric convection, we briefly consider four regimes of interest. As a reference, Fig. 14, shows Alfvén radius, A , as per Eq. (25), as function of gas accretion rate, , for various fixed strengths of dipolar surface magnetic fields, including the critical field crit as per Eq. (29).
(i) A < WD and WD < crit -surface magnetic fields have no impact on accretion flows or convection. Accretion is likely to be equatorial or spherical with convective surface spreading unimpeded. This could be the case for weakly magnetic white dwarfs or systems with high accretion rates. This is currently the default scenario for most metal-polluted white dwarfs with no measurable magnetic fields.
(ii) A > WD and WD < crit -surface magnetic fields affect gas out to the Alfvén radius, but atmospheric convection is unimpeded. Accreted gas is likely to be funnelled towards the poles and then transported via convective spreading, as described in this work, upon reaching the surface. This scenario is plausible if weak global magnetic fields, below current observational limits, are ubiquitous in white dwarfs.
(iii) A < WD and WD > crit -From Fig. 13, the lowest critical magnetic field in a white dwarf is crit ≈ 2 kG. Thus, an accretion rate in excess of > 10 16 g s −1 would be required to . The magnetic field corresponding to plasma-= 1 in the photosphere is shown in dashed blue. Thus, the blue shaded region shows where magnetic fields are unlikely to impede surface convection and the orange shaded region indicates where convective flux is likely to be heavily damped by magnetic fields. enable this scenario. For metal-polluted white dwarfs this scenario is therefore ruled out by observations.
(iv) A > WD and WD > crit -surface magnetic fields affect gas out to the Alfvén radius and significantly dampen convective flux. Thus accretion is likely to be pole-on with convective surface spreading potentially altered, but not necessarily inhibited. Landstreet & Bagnulo (2019) found that, within twenty parsec, 20 ± 5 per cent of DA white dwarfs have fields in excess of > 30 kG, which is above the critical field according to Fig. 13. Assuming this incidence does not change as a function of spectral type nor does it influence the presence of a planetary system, a significant fraction of evolved planetary systems would fall within this category.

DISCUSSION
One feasible way to constrain heterogeneities in the distribution of metals on a white dwarf surface is to look for photometric or spectroscopic variations on the rotational period. A spectroscopic and photometric analysis of 27 pulsating ZZ Ceti white dwarfs allowed Hermes et al. (2017a) to constrain the rotation of single white dwarfs with masses between 0.51-0.73 . They found that a mean period of 35 ± 28 hr could explain the rotational separation of pulsation frequencies.
Low level (few percent) photometric variations have been detected in non-pulsating white dwarfs and attributed to surface inhomogeneities that rotate in and out of view (Brinkworth et al. 2004;Kilic et al. 2015;Maoz et al. 2015;Hermes et al. 2017b). By combining Kepler photometry with Hubble Space Telescope (HST) ultraviolet (UV) spectroscopy for seven white dwarfs, Hallakoun et al. (2018) tested the hypothesis that low level variations in the photometry could be attributed to an accretion hot spot or inhomogeneous distribution of metals across the white dwarf surface. No such correlation could be claimed although their observations did not allow to rule out the hypothesis either. In contrast, photometric variability has been linked to magnetic white dwarfs Brinkworth et al. 2013;Reding et al. 2020;Gänsicke et al. 2020). In those cases the light intensity is thought to be inhomogeneous across the stellar surface, e.g. due to an uneven temperature distribution from a yet unexplained physical mechanism related to stellar magnetism (Tremblay et al. 2015b). As a consequence, photometric variability cannot yet be uniquely linked to evolved planetary systems.
The metal-rich DAZ white dwarf GD 394 was observed to exhibit a 25 per cent modulation in its EUV flux on a period of 1.15 d (Dupuis et al. 2000). The hypothesis that this variability was due to an accretion spot was tested by Wilson et al. (2019) with time-resolved HST UV spectroscopy. They found no evidence for variations in flux, inferred accretion rate and radial velocity over the observed period. Wilson et al. (2020) reported a 0.12% flux variation in optical photometric observations from the Transiting Exoplanet Survey Satellite (TESS), with a period consistent with that observed by Dupuis et al. (2000). They concluded that the observed variations may be explainable by a spot with an enhanced metal abundance, interaction with an orbiting planet or a magnetically-induced hot spot, but that more observations would be needed to confirm any of these three hypotheses.
The variability of accretion rates onto white dwarfs have also been studied. It was claimed that G29-38 showed evidence of a variable accretion rate based on a variation in the calcium line strength (von Hippel & Thompson 2007). However, by revisiting the observations Debes & López-Morales (2008) found no evidence for variation on timescales of days to years. Vanderbosch et al. (2020) have recently detected marginal calcium line variability at white dwarf ZTF J013906.17+524536.89 with transiting debris. However, they conclude that the spectroscopic variability is likely related to circumstellar debris transiting along the line of sight rather than photospheric variations. All in all, there is no robust evidence of surface metal heterogeneities at white dwarfs. This is in contrast with magnetic chemically peculiar main sequence stars where inhomogeneities are routinely detected (see, e.g., Sikora et al. 2019).
Although variations due to surface heterogeneities at white dwarfs have not yet been observed, variability of gas discs around white dwarfs is well established (Manser et al. 2020). In the majority of cases, the gaseous components of debris discs exhibit variations on the timescale of years to decades (Wilson et al. 2015;Manser et al. 2016a,b;Dennihy et al. 2018). Studies of dusty debris disc variability at white dwarfs have also been performed (Swan et al. 2019;Rogers et al. 2020). For a sample of white dwarfs with IR excess due to dust, Rogers et al. (2020) found no evidence of near-infrared -band variability above the 10% level for 32 objects with optical magnitudes smaller than 18. The upper limit on variation from the inferred dust disc was found to be 1.3% for objects brighter than magnitude 16. In contrast, Swan et al. (2019) used archival data from the Wide-field Infrared Survey Explorer (WISE) at redder near-infrared wavelengths to show that 69% of the 35 white dwarfs with IR excess due to dust show infrared variations more significant than 3 . This provides important information on disc evolution and replenishment Rafikov 2011;Metzger et al. 2012;Kenyon & Bromley 2017a,b), although there is no confirmed causal connection between disc variability and accretion rate variability and geometry.

Spectral line shapes
For any white dwarf with heterogeneous surface metal abundances, the spectral intensity must be calculated at each point of the surface and then integrated over the stellar disc, taking into account center-tolimb variations. All regions without metals will contribute nothing to the equivalent width of a given metal line. Therefore, the maximum depth of a metal line with respect to the continuum flux is roughly proportional to the fraction of the surface covered by metals. In other words, if metals are fully confined to, e.g., 10% of the surface, the flux decrement in any metal line should not be much larger than 10%, with the exact maximum value depending on the region of accretion, inclination and center-to-limb coefficient. Therefore the presence of deep or saturated photospheric metal lines in any white dwarf spectrum could be sufficient to rule out any significant surface inhomogeneities. We test this assumption using one example in the following section.

SDSS 1043+085
We consider in detail the expected manifestation of heterogeneity in surface metal abundances with an analysis of the metal-polluted DAZ white dwarf SDSS J104341.53+085558.2 (hereafter SDSS 1043+085). This star has an observed IR excess due to dust and it has been shown that disc-borne gas emission lines of the Ca II triplet are variable on a time scale of nine years (Gänsicke et al. 2007;Brinkworth et al. 2012).
For this experiment we use high-resolution UV spectroscopic observations of SDSS 1043+085. These were taken over two orbits with the HST Cosmic Origins Spectrograph (COS; Green et al. 2012) on 27 April 2015 as part of the GO program 13700 (Melis & Dufour 2017). These observations utilised the G130M grating, with central wavelength 1291Å, for a total exposure time of 5105 s. More details regarding the observational setup can be found in Melis & Dufour (2017). The authors found evidence of C, O, Al, Si, P, S, Fe and Ni in the COS spectrum. For the purposes of our experiment, we focus only on the strongest absorption features; the Si II resonance doublet at 1260 and 1264 Å.
We deduce the best-fitting parameters and Si abundance for this object by performing a two-step spectroscopic fit on the HST/COS observed spectrum shown in Fig. 15. We use a grid of white dwarf models computed with the atmosphere code of Koester (2010) with parameters log = 7.0-8.50 in steps of 0.1 dex and eff = 9000-30 000 K in steps of 200 K. The grid assumes a mixing length for convection of ML2/ = 0.8 and a pure-H composition (He/H = 0). This allows us to first update the white dwarf parameters by fitting the Ly-feature centred at 1216 Å. The strongest metal absorption lines are excluded during this step of the fit, and the masked wavelengths are shown in shaded grey in the figure. We use a Markov-Chain Monte Carlo (MCMC) approach with 100 walkers that run for 5000 steps. The likelihood was set by the chi-squared value to be −0.5 2 . We applied a Gaussian prior on the parallax based on the value from Gaia-DR2. We applied flat priors on the reddening, based on the value provided by STructuring by Inversion the Local Interstellar Medium (Stilism; Capitanio et al. 2017) for the distance to our system. We also applied flat priors to the effective temperature and surface gravity defined by the range of the model grid. We find the best-fitting parameters to be an effective temperature of eff = 16 536 +76 −70 K, surface gravity of log = 7.90 +0.04 −0.03 , parallax of Π = 5.98 +0.02 −0.05 mas and reddening of ( − ) = 0.01 ± 0.016 mag. These parameters are available in Table 5.
The second step of our fitting procedure was to deduce the silicon −0.05 mas and reddening ( − ) = 0.026 ± 0.003 mag. This fit was made around the Ly-feature, centred at 1216 Å, using an MCMC to find best-fit parameters for effective temperature, surface gravity, parallax and reddening. The 1-, 2-and 3-sigma tolerances on the derived parameters are shown in Fig. A1. The vertical grey shaded regions indicate wavelengths that were masked during the fitting procedure.  Fig. 16 shows the spectroscopic data from HST/COS (blue) and the synthetic spectrum (black) with parameters adopted from Table 5 and the best-fitting Si abundance. In order to test the hypothesis that metals may be heterogeneously distributed, we also show various synthetic spectra assuming the majority of metals are concentrated in a concentric circle on the visible stellar surface. The radius, , of this circle is defined in the 2D plane projection and thus must lie in the range WD > 0. To conserve the total amount of metals, the metallicity of the central spot is defined as where = 2 is the surface area of the spot and is the background metallicity, i.e., that of the remaining stellar surface. A value of = 1 corresponds to the original, best-fit metal abundance, and = 0 would correspond to a zero metal abundance. Fig. 16 shows synthetic spectra for four background abundance values; = 0.0001, 0.01, 0.1 & 0.5 in panels A-D, respectively. The total wavelength-dependent flux is computed by averaging the emission over the stellar surface using an eight-point Gaussian integration. This approach implicitly accounts for centre-to-limb darkening and should provide a reasonable estimate of the spectrum for SDSS 1043+085 under various levels of heterogeneity in its metal distribution. The radii of the central circles with enhanced metals are written in each of the lower four panels of the plot. The spot metallicity, , is also given alongside in brackets.
For a white dwarf at 16 500 K, our model predicts the sinking time to be at least 4 dex shorter than the lateral spreading time (see Fig. 8). Equations (18-20) predict an approximate background abundance in such a case to be vanishingly small (on the order of ∼ exp(−10 4 )). Thus, in Fig. 16 we favour the physical parameters of Panel A (a background abundance of = 0.0001) as being the most appropriate for this object. The lower three panels (B-D) are included as a demonstration of the possible spectroscopic appearance of a white dwarf with various metal-abundance configurations, although our model does not provide any mechanism to produce those higher background abundances for SDSS 1043+085.
Based on the synthetic spectra of Panel A, we find that, as predicted, the ratio of flux in the line cores to that in the continuum decreases with radius (and area). The background abundance contributes so little to the absorption lines that reproducing the deep flux cores is impossible, unless the spot covers most ( 65 per cent) of the surface (see the indigo line; = 0.8). We note that the non-linear dependence of the flux core depths on radius (or area) is a direct consequence of the centre-to-limb darkening. We conclude that the observed metals must be fairly homogeneously distributed across the surface, a finding that somewhat contradicts the predictions of our convectively-driven lateral spreading model. This may necessitate an accretion geometry favouring a relatively homogeneous deposition pattern (i.e., close to spherical), or additional mixing processes such as differential rotation. While we have selected a simple geometry, with the spot at the center of the disc with respect to the observer, the conclusion would be the same with a misaligned spot or an equatorial belt, with only small changes from center-to-limb effects.
We note that several other warm DAZ white dwarfs have deep or saturated metal lines in the UV. For instance, three DAZ white dwarfs observed with HST/COS , PG 0843+516, SDSS 1228+1040 and GALEX 1931+0117, have a deep Si II resonance doublet at 1265 Å with near-zero flux in the line core. Based on model atmosphere fits and radial velocities, it was confirmed that the origin of these absorption lines was the white dwarf photosphere. This suggests that nearly homogeneous metal coverage may be the norm rather than the exception, even at effective temperatures in the range 19 000-23 000 K, where only thermohaline mixing is expected to cause instabilities. While it is tempting to conclude that accretion must take place near-spherically in these objects, we can not rule out that we are missing another physical ingredient, e.g. another source of mixing in addition to convection and thermohaline instabilities.

Constraining accretion with X-rays
The geometry of accretion has so far received little attention in the context of the low accretion rates of metal polluted white dwarfs. Further work is therefore necessary to understand if the accretion process itself can be a source of mixing. Nevertheless, large enough accretion rates may lead to observational signatures such as X-ray fluxes, which is a promising method to test the accretion geometry in white dwarfs. Farihi et al. (2018) presented X-ray observations of several metalpolluted white dwarfs, including WD 1145+017 and G29-38. A reanalysis of XMM-Newton observations of G29-38 led to a new upper limit on the accretion rate based on the non-detection. Their model of X-ray production invoked a magnetically driven, pole-on accretion geometry. From optical spectropolarimetry they constrained the upper limit on the magnetic field of G29-38 to be = 0.7 ± 0.5 kG.
Based on the non-detection, their X-ray accretion model found an upper limit on the accretion rate of G29-38 to be = 1-6 × 10 9 g s −1 , consistent with the accretion rate of = 6.5 × 10 8 g s −1 inferred from spectroscopic abundances (Xu et al. 2014). The effective temperature of G29-38 is ≈11 400 K and the critical magnetic field for heavily damped convective flux in a DA white dwarf at this temperature is crit ≈ 3.5 kG. These parameters place G29-38 securely in scenario (ii) from Section 4.3; where A > WD and WD < crit , implying that the lateral spreading due to convection should be well described by the results presented in this work. An examination of Fig. 10 shows that the spreading and sinking time scales are comparable for a DA white dwarf at 11 400 K and thus local metal abundances may be expected to vary by a factor ∼3 across the surface.
An upper limit on the total X-ray flux at Earth was calculated, from G29-38, to be = 9.7 × 10 −15 erg cm −2 s −1 (Farihi et al. 2018). The plasma parameters adopted for their calculation included a plasma temperature of = 4 keV and a composition equivalent to that of the Sun. Accounting for interstellar absorption and the passband of XMM-Newton they found an upper limit on the flux to be = 8.0 × 10 −15 erg cm −2 s −1 . In Fig. 17 we show an estimate for the at-Earth X-ray flux from all DAZs in the Spitzer sample. This crude estimate is found by scaling the results from Farihi et al. (2018) by the inferred accretion rate and Gaia-DR2-measured parallax such that where is the parallax and the subscript G29 denotes the derived parameters for G29-38. CU19 demonstrated that the inclusion of convective overshoot in the modelling of debris accretion increases the inferred accretion rate by up to an order of magnitude. A similar hypothesis was also made by Bauer & Bildsten (2019) who claimed that the inclusion of thermohaline mixing could increase the accretion rates. For comparison with the X-ray fluxes shown in the figure, the EPIC detector on XMM-Newton may be sufficiently sensitive to detect a soft X-ray flux (0.5-2 keV) of 10 −15 erg s −1 cm −2 with exposures in excess of 300 ks (see Fig. 3; Watson et al. 2001). And according to the Chandra Proposers' Observatory Guide 1 , the point source sensitivity of the ACIS instrument on the Chandra X-ray Observatory is 4 × 10 −15 erg s −1 cm −2 , for an exposure time of 104 ks in the energy band 0.4-6.0 keV. Thus, current and future X-ray detecting facilities could allow the accretion rates at G29-38 and a handful of the other metal-polluted white dwarfs to be inferred via X-ray detections from sufficiently deep pointings. We intend to explore this line of research in a subsequent study. The X-ray emission from accreting white dwarfs would serve as a powerful constraint on the geometry and time dependence of accretion as well as convective overshoot effects at these objects.

CONCLUSIONS
We have presented the first diffusion coefficients of stellar surfaceplane transport derived from 3D tracer experiments for white dwarfs. We have utilised 3D RHD simulations for pure-hydrogen DA and pure-helium DB white dwarfs, which collectively span most of the parameter space of convective atmospheres, to examine the journey of trace metals once they are accreted to the white dwarf surface.  Figure 16. Spectroscopic observations of the Si II resonance doublet of SDSS 1043+085 taken with HST/COS shown in blue. The black curve shows the synthetic spectrum assuming the metals are distributed homogeneously across the surface. The purple-to-yellow lines show synthetic spectra if the same amount of metals are restricted to a central spot of a given radius. The radius refers to that of a circle in the 2D projection of the stellar surface, such that the radius of the spot lies in the range 0 < ≤ WD . The radii (0.1, 0.3, 0.5, 0.7, 0.8, 0.9 and 1.0 WD ) of the spots are given in the panels and correspond to surface areas of = 0.01, 0.09, 0.25, 0.49, 0.64, 0.81 and 1.00 WD , where WD refers to the visible surface area of the white dwarf in the 2D plane-projection. The metal abundance of a given spot, , is scaled inversely with the area using Eq. 30 and the spot metallicity fraction, , is given in brackets, next to the radius. Panels A-D vary only in the fixed background abundance of =0.0001, 0.01, 0.1 or 0.5. The top panel shows the spatial extent of each spot with enhanced metals, where the colours correspond to those of the synthetic spectra in the other panels. The observed spectrum (blue) has been smoothed using a Hann window function of size = 5, which made use of the Python package PyAstronomy (  Cooling Age [yrs] Figure 17. Bottom: Predicted total X-ray flux reaching the Earth from 48 polluted H-atmosphere white dwarfs with published accretion rates inferred from spectroscopic observations (Farihi et al. 2016). The fluxes were calculated using Eq. 31, by scaling the flux upper limit for G29-38 calculated by Farihi et al. (2018) with the inferred accretion rate and measured parallax from Gaia-DR2 (Gaia Collaboration et al. 2018). The fluxes computed using either the 1D or 3D model atmospheres (CU19) and are shown in open or filled circles, respectively. The green dashed lines indicate the impact of using either of the contending models for each object with a convective atmosphere ( eff < 18 000 K). Top: Change in accretion rate between the 1D and 3D overshoot models. The y-axis values are Δ = overshoot / 1D such that no overshoot correction corresponds to Δ = 1 (dotted line).
We find a strong temperature dependence in the ability of a white dwarf to homogenize its surface composition, a finding that somewhat challenges the often employed assumption that white dwarf surface layers can be assumed to be homogeneous in chemical composition. This assumption is verified only for cooler white dwarfs, a consequence of the increasing sinking time. In H-atmosphere stars we find that surface homogeneity can be explained by the horizontal diffusion due to convective motions below 10 500 ± 500 K where the stated uncertainty takes into account the accretion geometry and time dependence. For He-atmosphere white dwarfs the surface becomes homogeneous below about 30 000 K. We have demonstrated that hydrogen-rich DAZ white dwarfs above 13 000 K are predicted to show metal lines only in the regions of the stellar surface close to where accretion occurred. In contrast, UV observations of warm DA white dwarfs where strong and saturated metal lines are observed suggest rather homogeneous metal coverage at the surface. This discrepancy can only be explained by a near-spherical geometry of accretion or a yet unknown surface mixing process that would be several orders of magnitude faster than convective mixing. We have also provided a crude estimate of disc lifetimes around white dwarfs based on 3D RHD stellar atmosphere models, assuming a Jura-like disc and constant accretion rate, to be log( disc /[yr]) = 6.1 ± 1.4. This is only slightly increased compared to previous studies because the calculation depended on updated atmospheric models and a larger sample of polluted white dwarfs.

APPENDIX A: SDSS 1043+085 PARAMETER FITTING
This paper has been typeset from a T E X/L A T E X file prepared by the author.