# Convective Overshoot and Macroscopic Diffusion in Pure-Hydrogen Atmosphere
White Dwarfs

Tim Cunningham,1 Pier-Emmanuel Tremblay,1 Bernd Freytag,2 Hans-Günter Ludwig,3
and Detlev Koester4
1Department of Physics, University of Warwick, Coventry, CV4 7AL, UK
2Department of Physics and Astronomy, Uppsala University, Box 516, 751 20
Uppsala, Sweden
3Zentrum für Astronomie der Universität Heidelberg, Landessternwarte,
Königstuhl 12, 69117 Heidelberg, Germany
4Institut für Theoretische Physik und Astrophysik, University of Kiel, 24098
Kiel, Germany E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a theoretical description of macroscopic diffusion caused by
convective overshoot in pure-hydrogen DA white dwarfs using three-dimensional
(3D), closed-bottom, radiation hydrodynamics CO5BOLD simulations. We rely on a
new grid of deep 3D white dwarf models in the temperature range 11 400 K $\leq
T_{\mathrm{eff}}\leq 18\,000$ K where tracer particles and a tracer density
are used to derive macroscopic diffusion coefficients driven by convective
overshoot. These diffusion coefficients are compared to microscopic diffusion
coefficients from one-dimensional structures. We find that the mass of the
fully mixed region is likely to increase by up to 2.5 orders of magnitude
while inferred accretion rates increase by a more moderate order of magnitude.
We present evidence that an increase in settling time of up to 2 orders of
magnitude is to be expected which is of significance for time-variability
studies of polluted white dwarfs. Our grid also provides the most robust
constraint on the onset of convective instabilities in DA white dwarfs to be
in the effective temperature range from 18 000 to 18 250 K.

###### keywords:

white dwarfs – convection – hydrodynamics

††pubyear: 2019††pagerange: Convective Overshoot and Macroscopic Diffusion in
Pure-Hydrogen Atmosphere White Dwarfs–Convective Overshoot and Macroscopic
Diffusion in Pure-Hydrogen Atmosphere White Dwarfs

##  1 Introduction

For decades, understanding the diffusion of heavy elements through stellar
plasmas has been central to research in solar and stellar evolution (Michaud,
1970; Thoul et al., 1994), including the cooling of white dwarfs (Schatzman,
1945). In the latter case, large surface gravities ($7.0\lesssim\log g\lesssim
9.0$) imply the rapid settling of the initial composition into a structure
stratified according to atomic weight with an outer shell made of hydrogen and
helium.

The related scenario of metal accretion onto white dwarfs has now matured into
an extremely active field of research around evolved planetary systems (Veras,
2016). For a small fraction of these systems it is possible to study the
debris disk or gas around the white dwarf (Jura, 2003; Farihi et al., 2009;
Manser et al., 2016; Xu et al., 2018; Manser et al., 2019; Wilson et al.,
2019) but in most cases the primary evidence is metal pollution from tidally
disrupted asteroids in the stellar atmosphere (van Maanen, 1917; Zuckerman et
al., 2003; Gänsicke et al., 2012; Vanderburg et al., 2015). To transform
photospheric metal abundances into the parent body properties requires an
accurate model describing the volume of the stellar envelope over which this
abundance is prevalent. The accretion-diffusion model is usually employed,
where microscopic diffusion removes the heavy elements out of the observable
layers (Paquette et al., 1986a, b; Pelletier et al., 1986; Koester, 2009).

Constraining the mass and composition of the accreted material has wide-
reaching implications for the chemical compositions of other rocky worlds and
planetary evolutionary models (see, e.g., Zuckerman et al., 2007; Dufour et
al., 2012; Farihi et al., 2013; Farihi et al., 2016; Wilson et al., 2016; Xu
et al., 2017; Melis & Dufour, 2017). Since the emergence of space-based UV
spectroscopy the ability to detect metal pollution has significantly
increased, putting the fraction of degenerate stars with evolved planetary
systems at $\approx 50\%$ (Zuckerman et al., 2010; Koester et al., 2014).
These advances have allowed to describe the frequency and properties of
accretion events across white dwarf cooling age, mass, and spectral type, yet
the possibly observed correlations are not fully understood (Farihi et al.,
2012; Koester et al., 2014; Hollands et al., 2018), and systematic model
effects have been debated (Wachlin et al., 2017; Kupka et al., 2018; Bauer &
Bildsten, 2018, 2019). In this work we concentrate on theoretical aspects of
the accretion-diffusion scenario in the context of hydrogen-atmosphere DA
white dwarfs.

Young DA white dwarfs (cooling ages $t_{\mathrm{cool}}<10^{8}$ yr
corresponding to effective temperatures $T_{\rm eff}\gtrapprox 20\,000$ K)
have an atmosphere dominated by radiative energy transport. In this
temperature regime, radiative levitation can still maintain some heavier
chemical elements in the photosphere (Chayer et al., 1995a; Chayer et al.,
1995b), although ongoing accretion from planetesimals is invoked in a majority
of cases (Koester et al., 2014). Once DA white dwarfs have cooled sufficiently
to develop a convection zone, small characteristic convective turnover
timescales imply that metals become fully mixed within these turbulent layers,
and that microscopic diffusion only takes place at the base of the convection
zone. As this boundary layer is denser than the surface layers, this implies
longer settling times of the metals of up to 10 000 yrs (Koester, 2009; Bauer
& Bildsten, 2019).

Going back to the work of Prandtl & Tietjens (1925), convection has often been
modeled using 1D mixing-length theory (MLT). This model requires that the
Schwarzschild criterion for convective instability be satisfied and where it
is not, the models show no convective motions (Böhm-Vitense, 1958; Tassoul et
al., 1990). In the 1D MLT picture the base of the convection zone is at the
Schwarzschild boundary. Beneath this, the dominant process for mass transport
is that of microscopic diffusion. This microscopic diffusion comprises
contributions from gravitational settling, thermal diffusion, radiative
diffusion and diffusion driven by concentration gradients, unless the
diffusion concerns tracer particles, in which case the latter term is
neglected (Paquette et al., 1986b; Koester, 2009). The characteristic
velocities of microscopic diffusion near the Schwarzschild boundary are
$v_{\mathrm{diff}}\sim 10^{-7}$ km s-1. This can be up to seven orders of
magnitude less than the characteristic convective velocities of $\sim
1\,\mathrm{km\,s^{-1}}$ in the adjacent layer above.

The discontinuity between convective and non-convective layers in the 1D model
is unphysical (Spiegel, 1963; Roxburgh, 1978; Zahn, 1991) and multi-
dimensional numerical simulations have also confirmed this for the case of
white dwarfs (Freytag et al., 1996; Tremblay et al., 2015; Kupka et al.,
2018). Whilst the layers beneath the convectively unstable region are not able
to accelerate material to greater depths, convective cells accelerated at the
base of the unstable region can have sufficient momentum to penetrate the
deeper layers before dissipating their kinetic energy. This process is known
as convective overshoot and it has been shown from earlier 2D simulations to
be capable of mixing material (Freytag et al., 1996). Constraining precisely
the mass fraction of mixed material is crucial for our understanding of
accretion onto white dwarfs.

We present the first direct tests of mixing due to convective overshoot using
tracer particles and trace density in 3D radiation-hydrodynamics (RHD) CO5BOLD
simulations in the context of DA white dwarfs with shallow surface convection
zones. Building upon the work of Freytag et al. (1996, 2012), we estimate
macroscopic diffusion coefficients below the Schwarzschild unstable region
from direct 3D experiments (Sections 2-3). We employ the dependence of
macroscopic diffusion on mean convective velocities from these direct
experiments to infer the increase in convectively mixed mass for a wide range
of atmospheric parameters, using a grid of CO5BOLD simulations covering 11 400
K $<T_{\rm eff}<$ 18 000 K at $\log g=8.0$ (Section 4). We discuss the
implications on accretion rates and diffusion timescales (Section 5) and
conclude in Section 6.

##  2 Numerical Setup

The simulations presented in this work have been run using the 3D RHD code
CO5BOLD as described in Freytag et al. (2012). On a Cartesian grid, CO5BOLD
solves the coupled equations of compressible hydrodynamics and non-local
radiation transport, implicitly respecting conservation laws of energy, mass
and momentum. The vertical grid spacing is depth dependent, allowing for
better resolution of radiative transfer in upper layers, whilst the grid
spacing does not vary in the horizontal plane. For simulations presented
throughout this study radiation is handled with opacity and EOS tables from
Tremblay et al. (2013a, c), using both grey (Table 1) and non-grey schemes
(Table 2). The upper (surface) boundary is open to mass flows and radiation
(see Freytag 2017). The lower boundary is closed, requiring velocities to be
zero, and with a fixed radiative flux for a given effective temperature. The
four horizontal boundaries have periodic boundary conditions.

To increase stability around discontinuities a numerical reconstruction scheme
is required by the hydrodynamics solver. For all simulations presented in this
study the chosen scheme is designated as FRweno using a reconstruction by 2nd
order polynomials (Freytag, 2013). For the handling of trace density arrays
with CO5BOLD, which will be discussed in detail in the following section, a
piecewise-parabolic reconstruction scheme is used (Colella & Woodward, 1984).

The effective temperature is an input parameter for our simulations, as its
value determines the radiative flux put into the simulation lower boundary.
Following relaxation, the effective temperature is confirmed by the spatially
and temporally averaged emergent radiative flux. For brevity the effective
temperatures in Tables 1-2 will be referred to by rounding to the nearest 100
throughout this study, though any calculations involving the effective
temperature will use the more precise values shown in the tables.

###  2.1 Diffusion Coefficient Experiments

Table 1 shows the numerical setups for the simulations analysed in the direct
study of macroscopic diffusion in Section 3. All simulations have been built
from previously relaxed 3D simulations and given sufficient time to relax
after any parameter change.

Table 1: Numerical setup for grey opacity simulations analysed directly for
macroscopic diffusion with methods of CO5BOLD tracer density and path
integration. All simulations have a pure-hydrogen composition, closed bottom
boundary and $\log g$ = 8.0.

Sim. ID | $T_{\mathrm{eff}}$ | z | x,y | log $\tau_{\mathrm{R,top}}$ | log $\tau_{\mathrm{R,bot}}$ | $\Delta H_{p}$ | $z_{\mathrm{S,top}}$ | $H_{p,\mathrm{S,bot}}$ | $v_{\mathrm{z,S,bot}}$ | $[x,y,z]$ | $t_{\mathrm{diff~{}exp}}$ | $N_{\mathrm{tr~{}dens}}$
---|---|---|---|---|---|---|---|---|---|---|---|---
| [K] | [km] | [km] |  |  |  | [km] | [km] | [kms-1] | (grid points) | [s] |
A1 | 11 992 | 4.7 | 7.5 |  $-$4.28 | 3.63 | 5.19 | 1.06 | 0.53 | 4.48 | 250, 250, 250 | 6.02 | 10
B1 | 13 000 | 4.5 | 7.5 |  $-$3.35 | 3.01 | 5.06 | 0.66 | 0.40 |  | 150, 150, 150 | 9.05 | 20
B2 | 12 999 | 4.5 | 7.5 |  $-$3.33 | 3.01 | 5.02 | 0.43 | 0.38 | 2.78 | 250, 250, 250 | 3.16 | 10
C1 | 13 498 | 3.6 | 7.5 |  $-$3.19 | 2.69 | 4.37 | 0.52 | 0.37 |  | 150, 150, 150 | 6.20 | 20
C2 | 13 498 | 3.6 | 7.5 |  $-$3.18 | 2.69 | 4.37 | 0.51 | 0.38 |  | 250, 250, 150 | 7.91 | 20
C3 | 13 498 | 3.6 | 7.5 |  $-$3.20 | 2.69 | 4.39 | 0.50 | 0.35 | 1.52 | 250, 250, 250 | 8.41 | 10
C1-2 | 13 499 | 3.6 | 15.0 |  $-$3.19 | 2.69 | 4.39 | 0.50 | 0.37 |  | 150, 150, 150 | 8.14 | 20

Notes: The vertical (horizontal) extent of the simulation box is indicated by
$z$ ($x,y$). The number of pressure scale heights from the horizontally
averaged layer $\tau_{\mathrm{R}}=1$ and the simulation base is given by
$\Delta H_{p}$. All simulations have been fully relaxed prior to the addition
of tracer density at which point the diffusion experiment begins and its
duration is denoted by $t_{\mathrm{diff~{}exp}}$. The number of unique depths
at which tracer densities are added is given by $N_{\mathrm{tr~{}dens}}$. We
also include the geometric position of the upper Schwarzschild boundary,
$z_{\mathrm{S,top}}$, relative to the lower Schwarzschild boundary which is
set to be at the origin of the depth axis, i.e., where $z=0$ km. Denoted by
$v_{\mathrm{z,S,bot}}$ is the time-averaged RMS vertical velocity at the lower
Schwarzschild boundary. $T_{\rm eff}$ is derived from the time and spatially
averaged outgoing radiative flux. $\tau_{\mathrm{R,top}}$ and
$\tau_{\mathrm{R,bot}}$ are geometrical and time averages of the Rosseland
optical depth at the boundary layers of the simulation box. Characteristic
granule sizes for all simulations range from 1.0–1.2 km and data sampled from
simulations with period $\Delta t=0.01$ s.

###  2.2 Extended Grid of Closed 3D Simulations

Diffusion experiments using the direct methods laid out in Section 3 are
computationally expensive, in terms of both random-access memory during the
simulation run proper (Section 3.1) and data storage after the fact (Section
3.2). Furthermore, direct simulation of diffusion is hampered by numerical
waves as we shall discuss in Section 4, and thus in this work we also infer
indirectly about the macroscopic diffusion properties by using wave-filtered
velocities as a proxy. Towards this goal we employ a newly computed grid of
deep simulations across the effective temperature range 11 400 K $\leq$
$T_{\mathrm{eff}}$ $\leq$ 18 000 K. Numerical details are shown in Table 2.
The simulations detailed with 12 009 K $\leq$ $T_{\mathrm{eff}}$ $\leq$ 17 004
K had nominal input effective temperatures to the nearest 500 K. The input
effective temperatures for the two coolest (warmest) simulations were 11 400 K
and 11 600 K (17 525 K and 18 025 K). As an indicator of the quality of
relaxation of our simulations we find the maximum discrepancy between input
effective temperatures and those determined at the surface to be 0.4%, whilst
the majority of the grid ($T_{\mathrm{eff}}$ $\geq$ 12 000 K) has a maximum
difference of 0.1%.

These models are similar to the closed box simulations introduced in Tremblay
et al. (2013c, 2015) but were extended to deeper layers and are using a larger
number of grid points. The new grid has a vertical and horizontal resolution
of 30–50 m and 50–80 m, respectively. This is slightly less resolved than in
the previous grid, which had resolutions of $\sim$20 m and 30–50 m in the
vertical and horizontal directions, respectively; a compromise for probing
deeper layers. Furthermore, the grid was extended to significantly cooler and
hotter temperatures, with our coolest closed bottom simulation at $T_{\rm
eff}\approx 11\,400$ K overlapping with an open lower boundary simulation in
Tremblay et al. (2013c). All these simulations are sufficiently deep to
contain the entirety of the convectively unstable layers and deeper overshoot
layers, providing a unique opportunity to constrain convective velocities in
layers deeper than previously done with CO5BOLD, similarly to what has been
performed in Kupka et al. (2018).

Table 2: Extended grid of CO5BOLD simulations with pure-H composition and
$\log g$ = 8.0. All simulations have 2003 grid points, closed bottom
boundaries and non-grey opacities (Tremblay et al., 2013a).

$T_{\mathrm{eff}}$ | log $\tau_{\mathrm{R,top}}$ | log $\tau_{\mathrm{R,bot}}$ | $\Delta H_{p}$ |  $x$,$y$ | $z$
---|---|---|---|---|---
$[\mathrm{K}]$ |  |  |  | [km] | [km]
11445 |  $-$10.1 | 4.55 | 4.21 | $10.1$ | $10.0$
11651 |  $-$9.24 | 4.38 | 5.18 | $16.9$ | $9.8$
12009 |  $-$6.99 | 4.44 | 7.36 | $16.8$ | $10.6$
12514 |  $-$6.34 | 4.03 | 7.49 | $14.9$ | $8.7$
13005 |  $-$6.20 | 3.63 | 6.94 | $15.0$ | $7.1$
13503 |  $-$5.29 | 3.34 | 6.52 | $15.0$ | $6.1$
14000 |  $-$5.06 | 3.19 | 6.37 | $15.0$ | $5.8$
14498 |  $-$5.10 | 3.18 | 6.82 | $15.0$ | $6.2$
15000 |  $-$5.01 | 3.15 | 7.01 | $15.0$ | $6.4$
15501 |  $-$5.01 | 3.12 | 7.23 | $15.0$ | $6.7$
16002 |  $-$5.04 | 3.09 | 7.48 | $15.0$ | $6.9$
16503 |  $-$5.06 | 3.05 | 7.67 | $15.0$ | $7.0$
17004 |  $-$4.97 | 3.02 | 7.92 | $15.0$ | $7.1$
17524 |  $-$5.11 | 2.80 | 7.58 | $15.0$ | $6.4$
18022 |  $-$5.10 | 2.56 | 7.34 | $15.0$ | $5.7$

Notes: The number of pressure scale heights from $\tau_{\mathrm{R}}=1$ and the
simulation base is given by $\Delta H_{p}$. $\tau_{\mathrm{R,top}}$ and
$\tau_{\mathrm{R,bot}}$ are geometrical averages of the Rosseland optical
depth at the boundary layers of the simulation box. Further details on the
columns are given in the caption to Table 1.

##  3 Simulating Macroscopic Diffusion

The motivation for our research is to characterise the macroscopic diffusion
of trace metals in polluted white dwarfs for regions beneath the convection
zone, where 1D models currently predict no such mixing. In particular, we aim
at characterising the depth dependence of the efficiency of macroscopic mixing
in these overshoot regions. We do this by implementing a statistical, ensemble
averaged model for the mixing caused by overshoot in a finite region directly
beneath the convectively unstable layers, which we outline in the following.

The concentration of a fluid, $\phi$, can be described by the well known
diffusion equation, or Fick’s second law, in the absence of a source term and
with a spatially dependent diffusion coefficient as

$\frac{\partial\phi({\bf r},t)}{\partial t}-\nabla(D({\bf r})\nabla\phi({\bf
r},t))=0~{},$ (1)

where $D$ is the diffusion coefficient. We neglect any time dependence of the
diffusion coefficient as we are working with a system in a statistically
steady state - at least concerning transport properties in the vertical
direction. Given the spherical symmetry of a white dwarf surface we exclude
the possibility of horizontal dependence in the diffusion coefficient for the
purposes of this study. This allows the horizontal information from our 3D
simulations to be averaged, providing a robust statistical characterisation of
trace element distributions and ultimately giving a one dimensional diffusion
problem.

For a system undertaking a true diffusion process in one dimension the mean
displacement, $\langle z\rangle$, of an ensemble of particles, representing
the distance the concentration has travelled, is expected to evolve according
to

$\langle z^{2}\rangle=2Dt$ (2)

where $t$ is the time over which the system evolves. The objective of this
study is to calculate the depth dependent diffusion coefficients associated
with trace particles beneath the convectively unstable layers.

We use two methods to quantify the mean displacement of diffusing particles,
both of which will be described in the following sections. Method I utilises a
module built into CO5BOLD, designed for the study of dust formation in stellar
and solar environments, whilst Method II follows more closely the methodology
of Freytag et al. (1996) using tracer particles and serves as an independent
test of diffusion.

###  3.1 Method I: Tracer Density Arrays in CO5BOLD

Utilising the dust module included in CO5BOLD (Freytag et al., 2012) we can
add extra passive scalar density arrays to relaxed simulations. These density
fields have no mass, no opacity, and should thus only be advected by the local
velocity fields. We note, however, that CO5BOLD will adjust the time steps so
that the additional density fields can be reconstructed adequately, indirectly
influencing the overall numerical scheme. Microscopic diffusion velocity is
orders of magnitude smaller than macroscopic velocity fields, hence it is
ignored in these experiments. In any case we are only interested in the
stellar layers above which microscopic diffusion takes over. Via this method
we do not have actual tracer particles, but rather we store horizontally-
averaged number density distributions, providing significant advantages for
data storage and handling.

Additional and non-interacting density arrays are initially added to relaxed
simulations. To provide the most localised estimation of the diffusion
coefficient, a tracer density is inserted as a delta function (or horizontal
slab) in the z-direction such that

$\rho_{\mathrm{t}}(z)=\begin{cases}10^{5}~{}\mathrm{cm^{-3}},&\text{where
$z=z_{0}$}.\\\ 10^{-6}~{}\mathrm{cm^{-3}},&\text{where $z\neq
z_{0}$}.\end{cases}$ (3)

where $\rho_{\mathrm{t}}(z)$ is the tracer density. In our setup, Eq. (1) is
linear in $\phi$ so a horizontal slab is ideal for probing vertical mixing.
Since these density fields are massless the actual number density is only
relevant for the precision of the numerical schemes and rounding errors. The
large range between the peak and background tracer densities is chosen to
minimise any signature of net downward mixing due to the atmospheric density
gradient. Horizontally-averaged tracer densities from simulation C3 are shown
in the top panel of Fig. 1 for distributions used at the first and last
timesteps of the diffusion coefficient computation. The initial distributions
for the same simulation can be seen in the top panels of Fig. 2 where a
vertical cross section of a relaxed simulation at
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ to which the massless tracer density
(green) was added is shown.

Our initial experiments showed that standing waves appeared in the base layers
of our simulations, driven by convection and trapped between the convectively
unstable layer and bottom box boundary. This can be seen in Section 3.2 (see
Fig. 12) by the shallowing of the gradient at depths of
$z\lessapprox-1.5\mathrm{\,km}$ (grey shaded regions) in the vertical velocity
components of the simulation. We find that across all simulations the modes
are either above or below the acoustic cut-off frequency, suggesting both
p-modes, oscillations with pressure as the primary restoring force, and
g-modes, oscillations with gravity as the primary restoring force, can be
excited. In particular for the cooler models ($T_{\mathrm{eff}}\leq$12 000 K),
in the region beneath the unstable layers, we favour the interpretation of
g-modes (Freytag et al., 2010). Similar wave effects have also been observed
in deep white dwarf simulations by other groups (Kupka et al., 2018) where the
authors found evidence to suggest the presence of both p- and g-modes. A
comprehensive discussion on the treatment of these waves and extracting a
conservative estimate of the diffusion coefficient in this region can be found
in Section 4.2.

Preliminary results showed that the tracer densities were susceptible to
enhanced mixing in the wave region. Multi-dimensional simulations of AGB stars
have provided evidence that g-modes may be capable of mixing material (Freytag
et al., 2010). In order to ascertain whether the additional mixing is borne of
waves or numerical diffusion a comprehensive study of the wave properties
would be required, something beyond the scope of this work. Hence for the
purposes of direct diffusion experiments we prioritise probing the region
directly above the top of the wave region. The mixing timescale in the
convectively unstable region, with convective velocities of $\sim$ 1 km s-1,
is expected to be extremely short compared to timescales of accretion -
effectively providing instantaneous mixing. As such, we prioritise probing the
region beneath this, where $z<0$ km. The convention used throughout this work
is that our coordinate system is defined such that the lower Schwarzschild
boundary lies at $z=0$ km.

We focus our attention on the seven simulations detailed in Table 1, all deep
enough to fully enclose the convectively unstable layers and at least 4
pressure scale heights beneath the convectively unstable layers, defined by
the region where the entropy gradient with respect to depth is negative, i.e.,
$z_{\mathrm{S,top}}$ $>z>0$ km. The tracer density method is implemented on
all simulations in the table, from which diffusion coefficients are directly
derived. Simulations A1, B2 & C3, corresponding to $T_{\mathrm{eff}}$
$=12\,000,13\,000\,\,\&\,13\,500\,\mathrm{K}$, respectively, provide results
for macroscopic diffusion across a range of temperatures with a fixed grid
size ($250^{3}$). As an independent test of the results, these three
simulations are analysed using the method of path integration, the methodology
of which will be discussed in detail in Section 3.2.

Simulations C1 & C2, with $T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$, are
identical to C3, with the exception of grid resolution - and any numerical
scheme adjustments, such as time steps, arising from this change in grid size.
They are both analysed with the tracer density and path integration methods.
The box sizes of $150^{3}$ and $250^{2}\times 150$, respectively, serve as a
convergence test of the results with respect to spatial resolution.
Simulations B1 & B2, with $T_{\mathrm{eff}}=13\,000\,{\mathrm{K}}$ and grid
sizes $150^{3}$ & $250^{3}$, respectively, provide a similar test at a
different effective temperature. Finally, simulation C1-2 is identical to C1,
except for an extension in the x and y directions by a factor of two. This can
provide a further convergence test and also inform us on the nature of the
standing waves present in the bottom of the box.

####  3.1.1 Tracer Density Temporal Evolution

A given tracer density array, inserted as a horizontal slab, is rapidly
($\sim\mu s$) smoothed by the velocity field, where the less deep and more
vigorously convective layers promote this smearing to a greater extent. This
is demonstrated in Fig. 2 for simulation C3 (see Table 1) at
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$. Here we show snapshots of a 2D
vertical slice through 8 s of simulation time from the beginning of the
diffusion experiment. Tracer densities (green) inserted 1.2 km (left) and 1.0
km (right) beneath the unstable layers can be seen to be advected by the
convective flows (blue). The figure shows a range in tracer density of
$1\leq\rho_{\mathrm{trace}}/\mathrm{[cm^{-3}]}\leq 10^{5}$ with darker green
corresponding to a higher tracer density. Clearly some of the distribution
which began closer to the boundary of instability becomes incorporated into
the unstable layers ($z\geq 0$ km) within $t\approx 8\,\mathrm{s}$, though as
the tracer density is depicted logarithmically, the positional mean of the
distribution remains near the initial position. We show the same evolution in
Figs. 3 & 4 for simulations B1 and A1 (13 000 K and 12 000 K), respectively.

From the analytic solution of Eq. (1), in an ideal diffusion scenario an
initially narrow distribution is expected to evolve in the form of a Gaussian
for times $t>0$. The spatial extent of the density fields in our simulations
is quantified by the density-weighted standard deviation of the tracers,
$z^{*}_{\mathrm{rms}}$, which is given by

$z^{*}_{\mathrm{rms}}=\left(\frac{n_{z}}{n_{z}-1}\sum_{z}\left(\rho_{\mathrm{t}}(z)(z-\langle
z^{*}\rangle)^{2}\right)\middle/\sum_{z}\rho_{\mathrm{t}}(z)\right)^{1/2}$ (4)

where $\rho_{\mathrm{t}}(z)$ is the horizontally-averaged, depth dependent
tracer density, $n_{z}$ is the number of cells in the vertical dimension and
$\langle z^{*}\rangle$ is the density-weighted mean depth of the ensemble
given by

$\langle
z^{*}\rangle=\left(\sum_{z}z\rho_{\mathrm{t}}(z)\right)\cdot\left(\sum_{z}\rho_{\mathrm{t}}(z)\right)^{-1}$
(5)

To ascertain whether the ensemble of particles is evolving through a true
diffusion process we plot the density-weighted standard deviation,
$z^{*}_{\mathrm{rms}}$, as a function of time (see Fig. 1; middle panel). It
is expected that for true diffusion, the square of the ensemble spread will
evolve proportionally with time and we show the $\sqrt{t}$ fits made to the
$z^{*}_{\mathrm{rms}}$ evolutions with dashed lines. The depth dependent
diffusion coefficient, $D(z)$, is thus derived from Eq. (2) to be

$2\log_{10}(z^{*}_{\mathrm{rms}}(z,t))=\log_{10}(t)+\log_{10}(2D(z))~{}.$ (6)

(a)

(b)

(c)

Figure 1: Analysis of massless density arrays implemented with the CO5BOLD
tracer density for a simulation at $T_{\mathrm{eff}}$ = 13 500 K (Table 1; C3)
with $\log g$ = 8.0 and box size $250^{3}$. Top: Horizontally-averaged tracer
number density profile at first (solid) and last (dashed) time step used in
$\sqrt{t}$ fitting (see middle panel). Colours correspond to tracer densities
added at different depths. Middle: Evolution of the ensemble spread
characterised by the tracer density-weighted standard deviation,
$z^{*}_{\mathrm{rms}}$. Gaussian spread (solid) and a best fit of $\sqrt{t}$
(dashed). Bottom: Diffusion coefficients (circles) computed using Eq. 6.
Vertical velocity profiles $v_{\mathrm{z,rms}}$ and $v^{2}_{\mathrm{z,rms}}$
are shown in orange (dashed) and blue (solid), respectively. The depths on the
$x$-axis have been adjusted such that the lower Schwarzschild boundary lies at
$z=0$ km. We note that the full extent of the $x$-axis is considered
convectively stable under the Schwarzschild criterion, yet considerable mixing
is observed.

(a)  (b)

(c)  (d)

(e)  (f)

Figure 2: Demonstration of the tracer density implementation with CO5BOLD for
a simulation at $T_{\mathrm{eff}}$ = 13 500 K (Table 1, C3) with $\log g$ =
8.0 and box size $250^{3}$. Snapshots of a 2D vertical slice through the
simulation show the logarithmic tracer density (green) for two (left and
right) of the multiple density arrays added to the simulation. Only values of
$\rho_{\mathrm{trace}}\geq 1$ cm-3 are shown. Convective velocities are shown
in blue, where the magnitude is linear with line length. The depths on the
$y$-axis have been adjusted such that the lower Schwarzschild boundary lies at
$z=0$ km. Animated versions of this figure are currently available at
https://warwick.ac.uk/fac/sci/physics/research/astro/people/cunningham/movies/

(a)  (b)

(c)  (d)

(e)  (f)

Figure 3: Similar to Fig. 2 for a simulation at $T_{\mathrm{eff}}$ = 13 000 K
(Table 1, B1) with $\log g$ = 8.0 and box size $150^{3}$.

(a)  (b)

(c)  (d)

(e)  (f)

Figure 4: Similar to Fig. 2 for a simulation at $T_{\mathrm{eff}}$ = 12 000 K
(Table 1, A1) with $\log g$ = 8.0 and box size $250^{3}$.

The bottom panel of Fig. 1 shows the derived diffusion coefficients (circles)
for the region extending 1.5 km beneath the lowest convectively unstable layer
for simulation C3. Also shown are root mean square vertical velocity profiles
- $v_{\mathrm{z,rms}}$ (dotted, orange) and $v^{2}_{\mathrm{z,rms}}$ (solid,
blue) - time-averaged over the duration for which $\sqrt{t}$ fits are made
(middle panel - dashed). It can be seen that the decline of the diffusion
coefficient as a function of depth follows $D(z)=$ $v_{\mathrm{z,rms}}$
$\cdot$ 0.03 km in the near-overshoot region, extending 0.8 km beneath the
lower Schwarzschild boundary. For deeper layers, where $z<-0.8$ km, the
diffusion coefficient is well described by $D(z)=$ $v^{2}_{\mathrm{z,rms}}$
$\cdot$ 0.23 s.

Fig. 5 shows diffusion coefficients derived via this method for simulations B2
(left) at $T_{\mathrm{eff}}=13\,000\,{\mathrm{K}}$ and A1 (right) at
$T_{\mathrm{eff}}=12\,000\,{\mathrm{K}}$. We find that both of these
simulations exhibit similar behaviour to the
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ simulation (Fig. 1), with the depth
dependence of the diffusion coefficient being described by $D(z)=$
$v_{\mathrm{z,rms}}$ $\cdot d_{\mathrm{char}}$ in the near-overshoot region
and $D(z)=$ $v^{2}_{\mathrm{z,rms}}$ $\cdot t_{\mathrm{char}}$ in the deeper
layers.

We find the transition from $v_{\mathrm{z,rms}}$ to $v^{2}_{\mathrm{z,rms}}$
behaviour occurs at $z=$ $-$0.8, $-$1.2 and $-$1.6 km for simulations C3, B2,
and A1 (with $T_{\mathrm{eff}}=$ 13 500, 13 000 and 12 000 K), respectively.
The proportionality of transition depth with effective temperature is likely
due to the inverse proportionality of convective velocity, or kinetic energy,
and effective temperature in the range of these simulations (see Fig. 17). We
find that the characteristic distance required to fit the near-overshoot
region is relatively unchanging across the three temperatures, with
$d_{\mathrm{char}}\approx 0.03$ km, whilst the characteristic time required to
fit the far-overshoot region varies between $0.06\leq
t_{\mathrm{char}}/[\mathrm{s}]\leq 0.28$.

The vertical resolutions of simulations C3, B2 and A1 are 10–20 m, 10–30 m and
14–43 m, respectively, whilst the horizontal resolutions are all 30 m. As a
word of caution we point out that the characteristic distances invoked in the
near-overshoot region, $d_{\mathrm{char}}\approx 0.03$ km, are of the same
order as the grid resolution. At this stage we can not currently rule out the
possibility of some impact from numerical diffusion in the near-overshoot
region. This is independent of the diffusion coefficients derived with a
characteristic timescale in the far-overshoot region, which bares the most
significance for the ultimate determination of the mixed mass.

Previous results have shown that this characteristic time could be estimated
as $t_{\mathrm{char}}\sim H_{p}/v_{\mathrm{z,rms}}$ (i.e., Eq. (9) of Freytag
et al. 1996) where the pressure scale height, $H_{p}$, and RMS vertical
velocity are evaluated at the base of the convection zone. From an examination
of Table 1, which provides these two quantities evaluated at the lower
Schwarzschild boundary ($H_{p,\mathrm{S,bot}}$ and $v_{\mathrm{z,S,bot}}$), we
find estimates of a characteristic timescales of $t_{\mathrm{char}}=0.12$,
0.14 and 0.23 s for simulations A1, B2 and C3, respectively. These agree with
the characteristic times which rendered the best fit to the
$v^{2}_{\mathrm{z,rms}}$ lines to within a factor 0.5–2.0.

We find that all three simulations have sufficient motion from convective
overshoot to mix material for at least 1.5 km beneath the lowest formally
convective unstable layer, corresponding to at least 2.5 pressure scale
heights. We observe in all cases two distinct behaviours of the overshoot
diffusion coefficient. In the near-overshoot region, which extends 0.8–1.6 km
beneath the unstable layers, the diffusive efficacy decays with
$v_{\mathrm{z,rms}}$$\cdot d_{\mathrm{char}}$. In deeper layers the mixing
efficacy decays more rapidly, following instead a
$v^{2}_{\mathrm{z,rms}}$$\cdot t_{\mathrm{char}}$ profile.

The locations of the derived diffusion coefficients in Figs. 1 & 5 can be
compared with the y-axis positions in Figs. 2–4 where snapshots show the
evolution of tracer densities (green). In Fig. 2 the tracer densities were
placed in the far-overshoot region of simulation C3 (13 500 K) at
$z\approx-1.0$ km and $z\approx-1.2$ km. It is evident that overshoot plumes
are capable of penetrating this layer and mixing material in deeper layers
($z<-1.2$ km). These deeper layers correspond to the region which feels the
effect of waves trapped in the base of our simulations. The tracer density
method is susceptible to artificial diffusion driven by the increased
velocities here, typically leading to an overestimation of the diffusion
coefficients in this region. The increased velocities manifest as the upward
inflexion visible at $z\approx-1.9$ km in the right panel of Fig. 5. This
provides a lower limit on the depth at which a diffusion coefficient can be
directly derived via the method of tracer density. In Section 4.2 we will
discuss an approach which provides a conservative estimate of the decay of
mixing efficacy toward the stellar interior and beyond the region directly
simulated.

(a)  (b)

Figure 5: Similar to lower panel of Fig. 1 for simulations at
$T_{\mathrm{eff}}$ = 13 000 K (left, Table 1, B2) and $T_{\mathrm{eff}}$ = 12
000 K (right, Table 1, A1), both with $\log g$ = 8.0.

####  3.1.2 Effects of Numerical Resolution

The physical structures we wish to characterise include narrow overshoot
plumes which require grid points to be spaced sufficiently close. Furthermore,
numerical diffusion could adversely impact our results. We perform a
convergence test on the dependence of our results on spatial resolution by
using the tracer density to analyse simulations (C1, C2 & C3) at three
different resolutions ($150^{3}$, $250^{2}\times 150$ and $250^{3}$) and
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$. We note that the time steps are
changed accordingly to respect the Courant condition (Freytag et al., 2012).

Fig. 6 shows diffusion coefficients derived for all three simulations and,
qualitatively, we find that in the near-overshoot region the correlation of
the diffusion coefficient with the vertical velocity profile,
$v_{\mathrm{z,rms}}$, is insensitive to changes in spatial resolution for a
fixed box geometry. We do find a slight increase in the near-overshoot
diffusion coefficients for simulation C2, with grid size $250^{2}\times 150$.
It is possible that this is due to the increase in vertical grid spacing
compared to the $250^{3}$ simulation, meaning we cannot rule out that
numerical diffusion plays a role in this region and much higher resolution
simulations would be needed to comprehensively test this.

We observe that in the far-overshoot region (i.e., $z<-0.8$ km) the diffusion
coefficients tend to the $v^{2}_{\mathrm{z,rms}}$ profile as the horizontal
resolution increases. No significant change is seen for an increase in the
vertical resolution. We conclude that a resolution of $250^{2}\times 150$ is
adequate to resolve the physical processes involved with macroscopic diffusion
driven by deep convective overshoot for an atmosphere with
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ and $\log g$ = 8.

We now discuss a complimentary test to the above where the grid size is fixed
at $150^{3}$ and the box geometry is allowed to vary. The two simulations C1
and C1-2 from Table 1 have horizontal extents in the $x$ and $y$ direction of
7.5 and 15.0 km, respectively. This tests the spatial resolution in the
horizontal plane and as the vertical extent is kept constant across both
simulations this also serves as a test of varying aspect ratio. It can be seen
from Fig. 7 that for layers deeper than 1 km beneath the unstable layers the
diffusion coefficient exhibits similar behaviour, though decays less rapidly.
Within 1 km of the unstable layers the diffusion coefficient behaviour is
unchanged for an augmentation in the horizontal directions by a factor of two.
As a word of caution we point out that these convergence tests are not
exhaustive and, to be so, simulations with a significantly higher resolution
would be required. Were it for this tracer density method alone, these results
would be tentative, but a comparison with the path integration results -
presented in the following section - provides evidence that our estimates of
the size of the mixed region are physically robust.

(a)  Figure 6: Spatial resolution sensitivity test using tracer density
analysis for simulations C1, C2 and C3 and grid sizes $150^{3}$ (orange),
$250^{2}\times 150$ (red) and $250^{3}$ (blue), respectively, and
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ and $\log g$ = 8.0. The time-averaged
vertical velocity profiles $v_{\mathrm{z,rms}}$(dashed) and
$v^{2}_{\mathrm{z,rms}}$(solid) are shown for each simulation with the
diffusion coefficient extracted from each simulation (circles). The filled
error bars represent the standard deviation of the $\sqrt{t}$ fits.
Simulations C2 and C3 have been offset by $\log D$ = 2.5 and 5 for clarity.
(a)  Figure 7: The derived diffusion coefficients for simulations C1 and C1-2
from Table 1 with $T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$, $\log g$ = 8.0 and
grid size $150^{3}$. The simulations differ only in horizontal size where,
compared to simulation C1 (orange), the geometric extents in the $x$ and $y$
direction have been increased by a factor of two for simulation C1-2 (red).
The time-averaged vertical velocity profiles $v_{\mathrm{z,rms}}$(dashed) and
$v^{2}_{\mathrm{z,rms}}$(solid) are shown for each simulation with the
diffusion coefficient extracted from each simulation (circles). The filled
error bars represent the standard deviation of the $\sqrt{t}$ fits. Simulation
C1-2 has been offset by $\log D$ = 2.5 for clarity.

###  3.2 Method II: Path Integration

We discuss in the following the derivation of diffusion coefficients using a
method where the path of individual tracer particles is followed directly
through the simulation. In principle, path integration should give the same
results as the tracer density method presented in the previous section if we
were to horizontally average tracers, with any difference caused by numerical
schemes. Ostensibly, the way in which the diffusion coefficient is derived
from either method, path integration or tracer density, is also much the same.
We still look for a true diffusion process based on the spread evolving with
the square root of time.

Without an implemented tracer particle module, the path integration is
performed using the velocity field of pre-computed CO5BOLD snapshots. This
method is closer to the one presented in Freytag et al. (1996) and it allows a
more direct comparison to their results.

To ensure that the path of the seeds is fully interacting with the physical
processes within the simulation we are careful to sample velocity information
sufficiently frequently. To account for temporal changes in the velocity field
the sampling rate should be high enough that changes in the velocity field are
small. Spatial changes in the velocity field, which dictate the path taken by
the massless tracers, are handled by the integration time step being
sufficiently small.

Fig. 8 shows the depth dependent maximum velocity in units of grid points per
second for simulation C3, where the y-axis can be interpreted as the minimum
sampling frequency for any given depth. Ideally data would be sampled such
that massless particles are unlikely to travel much further than a few grid
points in one sampling period. The figure shows that to capture fully the
dynamical processes at every layer, velocity information would be required at
intervals of $\Delta t\sim 10^{-3}$ s, though this is only necessary if we
want to probe the most vigorous convective layers. It can be seen in Fig. 8
that an order of magnitude longer sampling period, $\Delta t\sim 10^{-2}$ s,
is sufficient and this also represents the chosen integration period.

Fig. 9 illustrates the vertical displacement of seeds for the simulation C3 at
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ with grid size $250^{3}$ over 8
seconds. At five depths beneath the unstable layers, $25\times 25$ seeds were
traced, at every tenth grid point. This lower number of particles helps to
clearly demonstrate the depth dependence of the spread of particles, with
those in deeper layers spreading less. For the path integration analysis
proper, which we will discuss in the next section, we use at least two orders
of magnitude more seeds at every depth.

For a given depth we place seeds at every grid point in a horizontal slice,
which for the highest resolution simulations presented here (A1, B2 and C3)
corresponds to 62500 seeds per depth. The three dimensional position,
$r_{x,y,z}(t)$, of a tracer beginning at grid point $r=[x_{0},y_{0},z_{0}]$ is
given by

$r_{x,y,z}(t)=r_{x_{0},y_{0},z_{0}}(t_{0})+\int v_{x,y,z}(t)~{}dt~{},$ (7)

where $v_{i}$ represents the $i^{\mathrm{th}}$ component of the velocity. This
integration is performed using a forward Euler method, where $dt=\Delta
t_{\mathrm{diff~{}exp}}$ is the time interval between consecutive snapshots of
extracted velocity information (see Table 1). Although the quantity of
interest is the depth dependence of diffusive efficacy it is necessary to
track the three-dimensional position of each tracer so that all the physical
processes present in our simulations have the opportunity to influence the
trajectory of the particles.

For the high resolution simulations, Eq. (7) is solved 62500 times at each
depth and time step such that we follow $\approx 10^{7}$ tracers per
simulation. As our primary focus is the region beneath the unstable layers we
present here tracers placed at every second vertical grid point for the region
defined by $z<0$ km, corresponding to $\approx 10^{6.5}$ tracer particles.

Figure 8: Depth dependence of maximum grid point displacement per second for
simulation C3 from Table 1 with $T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ and
$\log g$ = 8.0. Vertical (blue) and horizontal (orange) displacements are
shown independently. Displacements have been maximised over the final 25 ms of
the diffusion experiment. Figure 9: Tracer particle evolution for simulation
C3 with $T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ and $\log g$ = 8.0. Vertical
displacement of $25\times 25$ seeds (placed at every tenth horizontal grid
point for clarity), each placed at five different atmospheric depths in the
overshoot region. The colour intensity increases linearly with seed density.

As discussed in Section 3.1, to derive a diffusion coefficient we need to
quantify and fit the ensemble spread. With the method of path integration we
track the actual position of seeds and thus are able to access the mean square
displacement, $z_{\mathrm{rms}}$, given by

$z_{j,\mathrm{rms}}(t)=\left(\langle(z_{j}(t)-\langle
z_{j}(t)\rangle_{x,y})^{2}\rangle_{x,y}\right)^{1/2}~{},$ (8)

where $z_{j}(t=0)$ refers to the initial vertical position of a ’tray’ of
seeds placed at depth $j$ and $\langle z_{j}\rangle_{x,y}$ to the mean
vertical position of the seeds.

Eq. (8) is the standard deviation of the ensemble, making the retrieval of a
meaningful spread computationally simple and statistically robust. We can then
derive a diffusion coefficient from Eq. (2) as

$2\log_{10}(z_{\mathrm{rms}}(z,t))=\log_{10}(t)+\log_{10}(2D(z))~{}.$ (9)

The evolution of this quantity is shown in the top panel of Fig. 10 for
simulation C3. The final half of the total evolution time is fitted with lines
of $z_{\mathrm{rms}}\sim t^{1/2}$ (thick). All of the spread evolutions shown
correspond to tracers which began beneath the unstable layers, with the
distance beneath the lower Schwarzschild boundary indicated on the figure.

(a)

(b)

Figure 10: Path integration analysis of simulation C3 with
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ and grid size 2503, with plots akin
to the two lower panels of Fig. 1. Top: Standard deviation and $\sqrt{t}$ fit
(thick) for 15 of the 62 unique depths at which seeds were placed. The mean
initial depth of each set of $249\times 249$ seeds is indicated, for every
other distribution plotted. Bottom: Diffusion coefficients computed from the
method of path integration (solid green) are shown with vertical velocity
profiles ($v_{\mathrm{z,rms}}$, dotted-orange, and $v^{2}_{\mathrm{z,rms}}$,
solid-blue). Also shown are the results from the tracer density experiment for
the same simulation (black circles) from the bottom panel of Fig. 1. The
layers impacted by waves (see Section 4.2) are indicated by the grey shaded
region. No tracer density results are shown in this region due to the
artificial enhancement of the diffusion coefficients (see Section 3.1).

####  3.2.1 Characterising Diffusion

Our experiments with path integration allow us to derive diffusion
coefficients that can be compared to the method of density arrays (Section
3.1). We first discuss the results for simulation C3 with
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ and grid size $250^{3}$ in Fig. 10.
For the region within the convectively unstable layers
($0.5>z/[\mathrm{km}]\geq 0.0$) it is expected that seeds are mixed rapidly
(Freytag et al., 1996) with the spread increasing linearly with time. Just
below the convectively unstable region ($z<0.0$ km) seed ensembles are likely
to couple with the convectively unstable layers over a short timescale, $t<1$
s. This manifests in the top panel of Fig. 10 as a plateau in ensemble spread,
$z_{\mathrm{rms}}$. Once the upper edge of a seed ensemble reaches the
convectively unstable region, particles are carried in bulk by convective
velocities and the statistical, ensemble averaged model of diffusion is no
longer justified. Deriving a diffusion coefficient for the upper layers would
require fitting the region in time before the seeds couple with the
convectively unstable layer. There is significantly less data in this region
and for this reason no diffusion coefficients are computed here for $z>-0.2$
km.

The smoothness of the $z_{\mathrm{rms}}$ evolution demonstrates that the
resolution of this simulation provides a more than adequate sample of tracers
from which robust statistics are drawn. The fits by $\sqrt{t}$ functions are
excellent, implying that the average distribution of these tracers, whilst
moved by convective overshoot, can unequivocally be described as a diffusion
process. The diffusion coefficients computed from this evolution via Eq. (8)
are shown in the lower panel of Fig. 10 in green. This plot shows a strong
correlation of the diffusion coefficient with $v_{\mathrm{z,rms}}$ in the
region $-0.3>z/\mathrm{[km]}>-0.8$, which is in agreement with the results
from Method I (black circles) for the same simulation (Fig. 1; lower panel).
The diffusion coefficients computed by the method of path integration are also
in agreement in the far-overshoot region ($z<-0.8$ km) where a strong
correlation with $v^{2}_{\mathrm{z,rms}}$ is evident. This is again in
agreement with the prediction made by Freytag et al. (1996) where the authors
used the method of path integration for a simulation at
$T_{\mathrm{eff}}=13\,400\,{\mathrm{K}}$.

Figs. 11 & 12 show the results for simulations B2 and A1, respectively, which
have the same resolution as C3 and $T_{\mathrm{eff}}=13\,000\,{\mathrm{K}}$
and 12 000 K. It is again clear that the diffusion coefficients derived using
the method of path integration are in good agreement with the tracer density
method at all depths for which results are available. In the wave-dominated
region (grey shaded), where tracer density results proved inaccessible, the
method of path integration is able to separate diffusive action of the
overshoot motions from the essentially reversible motions due to waves. This
is evident from the flattening out of $v_{\mathrm{z,rms}}$ profile whilst the
diffusion coefficient continues to follow an exponential decay towards the
core. This effect becomes more obvious as the effective temperature decreases,
with simulation A1 at 12 000 K showing a sharp drop at the top of the wave
region. In simulations A1 and B2 the diffusion coefficient flattens out at
$z=-2.4$ and $-$2.5 km, respectively. This likely represents the limit of the
path integration method to disentangle the two kinds of motion and further
work would be required to ascertain the physical meaning of this behaviour.

(a)
Figure 11: Results of path integration analysis of simulation B2 with
$T_{\mathrm{eff}}=13\,000\,{\mathrm{K}}$ and grid size 2503. Akin to lower
panel of Fig. 10. (a)
Figure 12: Results of path integration analysis of simulation B2 with
$T_{\mathrm{eff}}=12\,000\,{\mathrm{K}}$ and grid size 2503. Akin to lower
panel of Fig. 10.

####  3.2.2 Effects of Numerical Resolution

We present here a test of the diffusion coefficient’s dependence on spatial
resolution derived via the method of path integration. Fig. 13 shows the
results for the two simulations C1 (top panel) and C2 (bottom panel) which
were prepared with grid size $150^{3}$ and $250^{2}\times 150$, respectively,
to compare to the $250^{3}$ simulation shown in Fig. 10 with the same
effective temperature. Qualitatively, both simulations retrieved similar
behaviour in the diffusion coefficient compared to simulation C3 - with a
diffusion coefficient scaling with $v_{\mathrm{z,rms}}$ for 0.8 km beneath the
unstable layers and with $v^{2}_{\mathrm{z,rms}}$ for deeper layers with
$z<-0.8$ km.

We observe a similar behaviour to the resolution test performed in Section
3.1.2 for the tracer density method with diffusion coefficient tending to
$v^{2}_{\mathrm{z,rms}}$ in the far-overshoot region as the horizontal
resolution increases. Qualitatively we also see no change as the vertical
resolution increases and thus conclude the results are derived from
sufficiently resolved simulations.

The three simulations are also in good agreement on the absolute value of the
diffusion coefficient which is most readily seen by examining the
characteristic distances and times required to fit the $v_{\mathrm{z,rms}}$
and $v^{2}_{\mathrm{z,rms}}$ profiles to the diffusion coefficients.

(a)

(b)

Figure 13: Similar to lower panel of Fig. 10 for simulations C1 (top) and C2
(bottom) from Table 1 with $T_{\mathrm{eff}}$ = 13 500 K, $\log g$ = 8.0 and
resolutions $150^{3}$ and $250^{2}\times 150$, respectively.

What is evident from the analysis of the five simulations discussed throughout
this section is that the vertical velocity profile is likely to be of great
significance to the characterisation of diffusive efficacy beneath the
convectively unstable layers, especially if one wishes to not perform a direct
diffusion experiment at every temperature of interest. The scaling of the
diffusion coefficient with $v_{\mathrm{z,rms}}$ immediately beneath the
unstable layers and with $v^{2}_{\mathrm{z,rms}}$ in deeper layers found in
both of tracer experiments is in good agreement with the dependence observed
by Freytag et al. (1996). Both approaches agree in terms of the observable
prediction, which is that trace elements are fully mixed in the overshoot
region $\approx$ 2.5 – 3.5 pressure scale heights immediately beneath the
convection zone, with mixing timescales still orders of magnitudes shorter
than those typical of microscopic diffusion. Furthermore, we emphasise that
both methods are impacted by standing waves and do not currently allow to
directly estimate the full extent of the mixed regions and the total mass of
trace metals in a white dwarf, which is therefore the goal of Section 4.

##  4 Results

In Section 3 we detailed the two methodologies used in this study to probe the
mixing efficiency of macroscopic diffusion due to convective overshoot, namely
the inbuilt CO5BOLD tracer density and path integration of seeds. While
standing waves near the bottom of the simulations prevented us to simulate
macroscopic diffusion over the full convective overshoot region, our results
highlighted the possibility of using $v_{\mathrm{z,rms}}$ and
$v^{2}_{\mathrm{z,rms}}$ as proxies for macroscopic diffusion. We therefore
use the extended grid of CO5BOLD simulations presented in Section 2.2 to study
$v_{\mathrm{z,rms}}$ in detail for a wide range of depths and $T_{\rm eff}$
values. Then, in Section 4.2 we develop a framework to lessen the contribution
of waves. We adopt this approach as fully quantifying the contribution of
waves to mixing is beyond the scope of this work. Though there exists some
evidence to suggest waves may contribute to the vertical mixing of material
(Freytag et al., 2010), our philosophy is that by removing their contribution
we can provide a lower limit on the mixed mass.

(a)

(b)

(c)

(d)

Figure 14: Dynamical and thermodynamic quantities extracted from the grid of
deep 3D simulations (Table 2) as a function of average Rosseland optical
depth. Top to bottom: Vertical component of convective velocities, horizontal
component of convective velocities, entropy and effective temperature.
Effective temperature is calculated as
$T_{\mathrm{eff}}=\sqrt[4]{F_{\mathrm{rad}}/\sigma}$ where $F_{\mathrm{rad}}$
is the radiative flux and $\sigma$ is the Stefan-Boltzmann constant.

(a)

(b)

(c)

Figure 15: Similar to Fig. 14 but for models exhibiting sizable convection
($11\,400\leq$ $T_{\mathrm{eff}}$ / [K] $\leq 13\,500$) and as a function of
geometric depth.

###  4.1 Dynamical and Thermal Properties

We rely on the 15 deep and closed simulations presented in Table 2 and
described in Section 2.2, covering the range 11 400 K $\leq$
$T_{\mathrm{eff}}$ $\leq$ 18 000 K. Selected mean quantities are shown in
Figs. 14-15 for all simulations computed in the extended grid.

Radiative energy transport becomes fully dominant for DA white dwarf
atmospheres above $T_{\mathrm{eff}}$ $\approx 14\,000$ K as illustrated by the
bottom panel of Fig. 14. This warmer end of the grid is convectively unstable
according to entropy profiles (Fig. 14, third panel from the top) but not
convectively efficient i.e., negligible energy is transported in convective
flows. In contrast, the three lowest effective temperature simulations
($T_{\mathrm{eff}}=$ 11 400, 11 600 and 12 000 K) can be seen in Fig. 15 to
have deep wells in their radiative flux profiles which implies a large
proportion of the energy carried in convective flows. Fig. 16 shows the time-
averaged Péclet number, evaluated at $\tau_{\mathrm{R}}$ = 1, as a function of
$T_{\rm eff}$. The Péclet number is given by Tremblay et al. (2013b) as

$\mathrm{Pe}=\frac{t_{\mathrm{rad}}}{t_{\mathrm{adv}}}=\frac{\langle\rho\rangle\langle
c_{p}\rangle v_{\mathrm{rms}}\tau_{\mathrm{e}}}{16\sigma\langle
T\rangle^{3}}\left(1+\frac{2}{\tau^{2}_{\mathrm{e}}}\right)~{},$ (10)

where $t_{\mathrm{rad}}$ is the radiative timescale and $t_{\mathrm{adv}}$ the
advective, or convective turnover, timescale, $c_{p}$ is the specific heat per
gram, $T$ the temperature, $\sigma$ the Stefan-Boltzmann constant, and
$\tau_{e}=\langle\kappa_{\mathrm{R}}\rangle\langle\rho\rangle H_{p}$ the
characteristic optical thickness of a disturbance with size $H_{p}$, for
$\kappa_{\mathrm{R}}$ the Rosseland mean opacity per gram. It represents the
ratio of energy being carried by bulk movement of fluid and that being carried
by radiation. For Pe $\ll$ 1, convection will not be responsible for carrying
much of the energy. It can be seen from Fig. 16 that Pe < 1 for
$T_{\mathrm{eff}}$ $\geq$ 12 000 K and thus it describes the transition to
inefficient convection at larger temperatures. We emphasise that it does not
imply a change in the mixing properties of the atmosphere at that temperature.
In fact, Tremblay et al. (2015, see figure 13) have shown that
$v_{\mathrm{z,rms}}$, and thus mixing capabilities, remain large for 3D
simulations up to $T_{\mathrm{eff}}$ $=14\,000$ K, the hottest model in their
grid at $\log g=8.0$. This behaviour is in agreement with the mixing-length
theory, which predicts that thermally inefficient but dynamically active
convective zones are found in DA white dwarfs up to 18 000 K, with the exact
$T_{\mathrm{eff}}$ value strongly dependent on the assumed mixing-length
parameter (ML2/$\alpha$ = 0.7 in Tremblay et al. 2015). Furthermore, the 1D
model can only suggest whether or not there is mixing and cannot provide
realistic mixed masses due to its neglect of overshoot.

Figure 16: Péclet number evaluated at $\langle\tau_{\mathrm{R}}\rangle=1$ as a
function of effective temperature and averaged over the final 50 ms of the
deep 3D grid presented in Table 2.

The entropy profiles of our extended grid of 3D simulations in Fig. 14
demonstrate that pure-hydrogen atmospheres become unstable to convection for
$T_{\mathrm{eff}}$ $\lessapprox 18\,000$ K. We attempted a $T_{\mathrm{eff}}$
$=18\,250$ K simulation, not shown on the panel, but it did not develop a
convection zone. While this is similar to the 1D prediction of the onset of
convection assuming a ML2/$\alpha$ = 0.7 parameterisation (Tremblay et al.,
2015), the 3D simulations provide a more robust threshold. Fig. 17 shows the
maximum convective velocity in the convectively unstable region, time-averaged
over the final 2.5 s, for all simulations in our deep grid (see Table 2) as
well as Tremblay et al. (2015). It can be seen that the grid presented in this
work fully connects with the previous, shallower one. It is also observed that
1D MLT under-predicts convective velocities for lower temperatures
($T_{\mathrm{eff}}$ $<14\,000$ K) and over-predicts for higher temperatures
($T_{\mathrm{eff}}$ $>14\,000$ K), although this could be corrected by using a
variable mixing-length parameter as a function of $T_{\rm eff}$. The first
convective model at $T_{\mathrm{eff}}$ = 18 000 K is already showing large 1D
velocities, illustrating the sharp and unphysical transition between
convective and radiative atmospheres. In 3D, the warmer convective models have
much lower velocities, in line with a smoother onset of convection in evolving
DA white dwarfs. This effect may be enhanced by numerical viscosity which, in
the regime where radiative energy transport is significantly more efficient
than convection, may reduce the 3D convective velocities.

Figure 17: Maximum convective velocities, time-averaged over the final 2.5 s
for simulations in Table 2 (green, half-filled) with 13 503 $\leq$
$T_{\mathrm{eff}}$ / [K] $\leq$ 18 002. Also shown are previous 3D results
(orange, filled) and 1D ML2/$\alpha$ = 0.7 mixing-length theory predictions
(blue, open) from Tremblay et al. (2015).

That the atmosphere is unstable is of critical importance for the capacity for
DA white dwarfs in this temperature range to mix material, specifically
accreted metals. Fig. 14 demonstrates that all simulations have large
velocities, although with a contribution both from convection and waves. It
therefore requires further investigation to fully characterise the mixing
capabilities.

###  4.2 Mixed Mass

The depth dependence of the diffusion coefficient just above the wave region
is well described by the square of the RMS vertical velocity,
$v^{2}_{\mathrm{z,rms}}$, (see Section 3) with the characteristic timescale
$t_{\mathrm{char}}\approx 0.18$ s taken from the
$T_{\mathrm{eff}}=13\,500\,{\mathrm{K}}$ simulation (see Fig. 10). Ideally our
simulations would extend to sufficiently deep layers that convection driven
velocities decay enough that microscopic diffusion takes over, bridging the
discontinuity between macroscopic and microscopic diffusion highlighted by
Tremblay et al. (2017). This would be computationally intractable given the
typical microscopic diffusion velocities of $v_{\mathrm{diff}}\sim 10^{-7}$ km
$\mathrm{s}^{-1}$ and a significant hinderance to this approach is the
presence of waves in the base of the simulation.

The upward inflexion in the velocity profile at the bottom of the overshoot
region highlighted in earlier figures (see Figs. 10 & 12, and upper panel of
Fig. 15) is due to these standing waves. Multi-dimensional simulations of
brown dwarf atmospheres have provided evidence that waves may be capable of
mixing material (Freytag et al., 2010). However, using
$v^{2}_{\mathrm{z,rms}}$ in the wave region as a proxy for macroscopic
diffusion could lead to an overestimation as standing waves may be amplified
by boundary conditions. Alternatively, extrapolating the diffusion coefficient
from the top of the wave region, e.g. with an exponential decay as a function
of the pressure scale height (Tremblay et al., 2015), would lead to an
extremely shallow decay of macroscopic mixing as a function of depth that may
not be realistic (Kupka et al., 2018). This suggests that more investigation
is needed to obtain sensible mixed masses.

####  4.2.1 Filtering

We employ a filtering technique to remove the waves from the velocity profile
which has been implemented in convective studies of M-type, main-sequence and
pre-main-sequence objects (Ludwig et al., 2006).

Fig. 18 shows the logarithmic vertical velocities at a layer within each of
the three regions with distinct physical characteristics; the convectively
unstable region (top), overshoot region (middle) and layer where plumes and
waves overlap (bottom), for simulation B1 from Table 1 with
$T_{\mathrm{eff}}=13\,000\,{\mathrm{K}}$, $\log g$ = 8.0 and grid size 1503.
In the convectively unstable region the convective cells, with characteristic
size $d_{\mathrm{gran}}\approx 1$ km, can be seen with upward (teal) and
downward (bronze) velocities. Beneath this layer in the overshoot region one
can identify around twenty distinct overshoot plumes, which in the horizontal
slice manifest as small circles with characteristic size $\sim 100$ m. These
circles are surrounded by large scale patches of size $\sim 5$ km and
relatively low velocity. These patches correspond to the standing waves. In
the region below the deepest overshoot plume all that remains is the
contribution from periodic oscillations.

The contrast in spatial extent of the overshoot plumes and standing waves is
of great benefit as the two occupy different regions of the power spectrum in
spatial frequency. This is depicted in Fig. 19 where the power as a function
of horizontal wavenumber is shown for vertical velocities. The large scale
wave features manifest as the blue peak at low wavenumber ($k<4\cdot
10^{-5}$), whilst the overshoot plumes correspond to the purple peak at higher
wavenumbers. The masking of the blue region (low wavenumbers) allows the
removal of the major contribution from waves in frequency space, while leaving
the overshoot plumes ostensibly untouched.

Figure 18: Logarithmically coloured vertical velocity snapshots with (right)
and without (left) spatial frequency filtering for regions demonstrating:
convective instabilities (top), overshoot plumes (middle) and transition from
plumes to waves (bottom), for simulation B1 with
$T_{\mathrm{eff}}=13\,000\,{\mathrm{K}}$, $\log g$ = 8.0 from Table 1. The
vertical distances from the lower Schwarzschild boundary are indicated in the
panels. (a)
Figure 19: Power spectrum for spatial frequencies present in lower panel of
Fig. 18. The unfiltered power spectrum (lower left panel of Fig. 18) can be
interpreted as a combination of the blue and purple regions. The filtered
power spectrum (lower right panel of Fig. 18) corresponds to the purple region
only. The vertical distance from the lower Schwarzschild boundary is indicated
in the panel and corresponds to the layer showing a transition from plumes to
waves.

Fig. 20 shows the diffusion coefficient computed from the
$v^{2}_{\mathrm{z,rms}}$ profiles, time-averaged over the final 2 s, of the
simulations from Table 2 with $T_{\mathrm{eff}}=12\,000\,{\mathrm{K}}$ and 13
500 K. The characteristic time required for the computation is taken from the
results of simulation C3 (see Fig. 1). We note that the scaling of the
diffusion coefficient with $v_{\mathrm{z,rms}}$ in the near-overshoot region
is of little relevance for deriving the total mixed mass. The effect of
removing increasingly small scale structures from the power spectrum is shown
by the lines becoming bluer. There is a clear asymptotic behaviour which
emerges below enclosed mass $\log M_{\mathrm{above}}/M_{\mathrm{star}}=-14$ in
both cases. The simulations also show significant influence from the lower
boundary condition which enforces zero mass flux. Clearly, filtering in
$k$-space against low spatial frequencies removes a large contribution from
waves yet results in little disturbance of the velocities within the
convective zone or the upper overshoot region. Any increase in filtering
beyond what is shown in Fig. 20 leads to a significant loss of velocities
within the convection zone which is an unphysical scenario we seek to avoid.

In the study by Ludwig et al. (2006) the filtering technique was also
implemented in the temporal domain, ultimately making a cut in $k$\- and
$\omega$-space, where $k$ is the horizontal wavenumber. In our white dwarf
atmospheres the standing waves have typical periods of 0.1 s (see also Kupka
et al., 2018), similar to the typical decay time of granulation within the
convection zone (Tremblay et al., 2013c). Therefore, it is difficult to
disentangle convection from waves in $\omega$-space. Given the effectiveness
of the $k$-space filtering technique discussed above, we simply neglect time
filtering in the following.

(a)  (b)

Figure 20: Diffusion coefficients as a function of enclosed stellar mass from
deep simulations of Table 2 with $T_{\mathrm{eff}}$ = 12 000 K (left) and 13
500 K (right). The macroscopic diffusion coefficients computed from the
directly simulated $v^{2}_{\mathrm{z,rms}}$ profile are shown as a solid
purple line. Those computed with an increasingly filtered
$v^{2}_{\mathrm{z,rms}}$ profile (i.e., increasingly small structures removed)
are shown in bluer colours, with adjacent lines increasing the filtering in
Fourier space by $\Delta k=1$ in box size units. The characteristic time scale
of 0.18 s found in Section 3 is employed. Our proposed extrapolation of
macroscopic diffusion coefficients is shown as a black dot-dashed line. The
convectively unstable region defined by the Schwarzschild criterion is
indicated by the orange shaded region. Also shown are the microscopic
diffusion coefficients for Mg, Ca and Fe from updated tables of Koester (2009,
dotted lines), taken as characteristic elements. The base of the mixed region
is where the macroscopic and microscopic diffusion coefficients are equal.
Diffusion coefficients at the base of the mixed regions defined in 1D (black
circle) and with our proposed parameterisation of convective overshoot for the
three characteristic elements (coloured circles) are also presented.

We find that the asymptotic filtered $v^{2}_{\mathrm{z,rms}}$ is linear in
$\log D$ \- $\log M_{\rm above}/M_{\rm star}$ space, allowing us to perform a
linear fit immediately beneath the start of the wave region. We can then
extrapolate this behaviour to deeper layers that are out of reach with 3D
simulations. Fig. 20 demonstrates that the decay of convective velocities in
the wave region is much more rapid with depth than in the overshoot region
above. Indeed the velocity scale height in the wave region (where overshoot
structures still exist) is much less than the pressure scale height,
reflecting similar behaviour which is observed in simulations of convection in
brown dwarf atmospheres (Freytag et al., 2010), but we cannot confirm whether
this change in convective velocity profile is impacted by waves, e.g. with
wave damping convective overshoot. For this reason and also our neglect of
waves as a possible contributor to mixing, our estimate of the mixed mass is a
lower limit.

We estimate the minimum amount of mixed mass in accreting white dwarfs by
looking for the intersection of the asymptotic filtered macroscopic diffusion
coefficient with microscopic diffusion coefficients. In Fig. 20 we estimate
these mixed masses using the magnesium, calcium and iron diffusion
coefficients (dotted) from Koester (2009); using updated
tables.111http://www1.astrophysik.uni-
kiel.de/$\sim$koester/astrophysics/astrophysics.html Here it is important to
emphasise that 3D corrections to mixed masses can be read directly from Fig.
20. The crude 1D estimate is simply the intersection of the base of the
convection zone with the diffusion coefficient (black circle), and the more
robust 3D picture is based on the intersection of the asymptotic filtered
macroscopic diffusion coefficient with the microscopic diffusion coefficients
(coloured circles). The fact that convective velocities decay rapidly with
depth in the wave region implies that while we extrapolate over a sizeable
range in velocity space, it corresponds to a fairly modest additional amount
in mixed mass. In other words, our direct simulations in Section 3 have
robustly shown that a sizeable amount of mass is unambiguously mixed in the
overshoot region, and the additional amount added by our filtering procedure
is fairly moderate.

The convection zone mass after one includes the full extent of the
convectively unstable and overshoot regions is shown in Fig. 21 with solid
lines for Mg (blue), Ca (purple) and Fe (red) and described in Table 3. The
figure also shows the mass at the base of the convection zone in 3D described
from the Schwarzschild criterion for our deep grid of simulations (green
line). Those are in fairly good agreement with the results of Tremblay et al.
(2015) (green triangles) and 1D structures (grey and orange points), and this
corresponds to the mixed mass in the 1D picture, i.e. neglecting the
contribution from overshoot. The lower limit on the 3D correction in mass is
at least 1 dex but can reach $\approx$ 2 dex for the intermediate $T_{\rm
eff}$ values in our grid. We will go on to consider the implications of this
result for the accretion-diffusion picture.

Recently Kupka et al. (2018) also presented a prediction for the range of
mixed mass for a simulation at $T_{\mathrm{eff}}=11\,800\,{\mathrm{K}}$ to be
$-12.7\gtrapprox\log q\gtrapprox-13.9$ based on the plume penetration depth
and horizontal to vertical velocity ratios, respectively. An interpolation of
Fig. 21 reveals this to be in good agreement with the predicted mixed mass
from our 3D grid which estimates the mixed mass at
$T_{\mathrm{eff}}=11\,800\,{\mathrm{K}}$ to be $\log q=-13.0$. We also find
agreement with the predicted size of the convectively unstable zone that the
authors found to be $\log q=-15.2$, which compares well with the green line in
Fig. 21.

Figure 21: Mixed mass in units of enclosed mass, $\log q=\log M_{\rm
above}/M_{\rm star}$, as a function of effective temperature. Mixed masses
determined by the intersection of macroscopic and microscopic diffusion
coefficients (see Fig. 20) for Mg (blue), Ca (magenta) and Fe (red) are shown
in solid for all simulations detailed in Table 2. Also shown is the 3D
definition of the unstable region from the grid presented in this study
(green, solid) and selected models from the earlier lower resolution grid
(green, triangles) from Tremblay et al. (2015). This corresponds to the
definition of mixed mass without overshoot. The two coolest models plotted
from the latter grid are open bottom simulations, and show good agreement with
the closed bottom simulations calculated in this work. In comparison, the
orange and grey lines show the 1D mixed mass of the unstable layers according
to 1D structures of Koester and Fontaine, respectively, using
ML2/$\alpha=0.8$.

##  5 Discussion

Chemical abundances in polluted white dwarf atmospheres as well as
characterisation of circumstellar dust and gas disks for a fraction of those
are the main observational properties with which we can constrain evolved
planetary systems. On the one hand, the chemical properties and time
variability of circumstellar material is independent of the white dwarf
convection model. These systems can be used to understand the variability of
predicted accretion rates onto the white dwarfs (Wyatt et al., 2014). On the
other hand, it is the accretion-diffusion scenario, currently under the 1D MLT
assumption of no convective overshoot, that is used to predict parent body
chemical compositions, the total accreted mass over the characteristic
diffusion timescale, and in some cases, estimate the time elapsed since the
last accretion event (Hollands et al., 2018). In the case of the hydrogen
abundance in helium-atmospheres, one can also infer about the full past
history of accretion since hydrogen does not diffuse out of the convective
layers (Gentile Fusillo et al., 2017). Time variable abundances in the
photosphere, although so far not observed (Farihi et al., 2018; Debes & López-
Morales, 2008), could also be used to understand the accretion process. All of
those inferred quantities need to be re-evaluated in light of 3D convection.
We have determined that the mixed mass, or the total accreted mass over a
characteristic diffusion timescale, is 1–2.5 orders of magnitude larger than
the crude 1D estimate and we dedicate this section to an overview of the other
quantities derived from the accretion-diffusion model.

One starting point is the impact of convective overshoot on the settling
times, $t_{\mathrm{diff}}$, of trace elements in polluted white dwarfs.
Accretion calculations are usually performed under the assumption that a
steady state has been reached between accretion at the white dwarf and
diffusion out of the base of the convection zone, for at least a few diffusion
timescales. This implies that the abundance in the mixed region is homogeneous
and in this scenario the accretion rate is given by

$\dot{M}_{i}=\frac{M_{\mathrm{cvz}}X_{i}}{t_{\mathrm{diff},i}}=4\pi r^{2}\rho
v_{\mathrm{diff}}$ (11)

where $M_{\mathrm{cvz}}$ is the mixed mass and, for an element $i$,
$\dot{M}_{i}$ is the accretion rate, $t_{\mathrm{diff},i}$ the diffusion
timescale, $X_{i}$ the elemental atmospheric abundance, $\rho$ the density at
the base of the mixed region and $r$ the stellar radius.

We use the 1D diffusion velocities presented in Koester (2009) to calculate
the change in settling timescale when one redefines the base of the mixed
region. Fig. 22 shows 3D corrections to diffusion timescales (top) and
accretion rates (bottom). Convective overshoot implies that microscopic
diffusion takes over as the dominant process of mass transport in deeper,
denser layers where diffusion timescales are longer. For an inference of
accretion rates we have the competing effects of larger masses and larger
diffusion timescales that partially cancel each other out as the base of the
mixed layer goes deeper. The net result from this competition is an increase
in accretion rate of up to an order of magnitude (see bottom panel of Fig.
22), which is not as large as the increase in mixed masses.

(a)

(b)

Figure 22: Diffusion timescale (top panel) and inferred accretion rates (lower panel) when using the mixed masses, shown in Fig. 21, borne from diffusion coefficients derived from our 3D grid using $v^{2}_{\mathrm{z,rms}}$ (solid). For simplicity the abundance, $X_{i}$, in Eq. (11) used to infer an accretion rate is assumed to be unity. Also shown are 1D structures from Koester (solid-circles). All quantities are provided for Mg (blue), Ca (magenta) and Fe (red). The strong upward inflexion in the 1D accretion rates for models with $T_{\mathrm{eff}}<12\,000$ K is a symptom of the strong downward inflexion seen in the mixed masses for the same models (see Fig. 21). Table 3: Input parameters for accretion-diffusion calculations derived for our full grid of deep 3D simulations (see Table 2). The mixed mass in units of enclosed mass, $M_{\mathrm{cvz}}$, diffusion velocity, $v_{\mathrm{diff}}$ and diffusion coefficient, $D$, are all evaluated at base of the 3D mixed region. Values are shown for mixed regions defined by diffusion coefficients with $D=v^{2}\cdot t_{\mathrm{char}}$ (see Fig. 21). For brevity, all values are expressed as a mean across those computed for magnesium, calcium and iron. $T_{\mathrm{eff}}$ | $\log\langle M_{\mathrm{cvz}}\rangle$ | $\log\langle v_{\mathrm{diff}}\rangle$ | $\log\langle D\rangle$
---|---|---|---
$[\mathrm{K}]$ |  |  |
11445 |  $-$12.82 |  $-$7.62 |  $-$10.28
11651 |  $-$12.89 |  $-$7.60 |  $-$9.97
12009 |  $-$13.19 |  $-$7.43 |  $-$9.90
12500 |  $-$13.41 |  $-$7.30 |  $-$9.64
13005 |  $-$13.64 |  $-$7.17 |  $-$8.88
13503 |  $-$13.89 |  $-$7.03 |  $-$8.76
14000 |  $-$14.09 |  $-$6.91 |  $-$8.68
14498 |  $-$14.58 |  $-$6.63 |  $-$8.51
15000 |  $-$14.49 |  $-$6.68 |  $-$8.52
15501 |  $-$14.60 |  $-$6.61 |  $-$8.47
16002 |  $-$15.01 |  $-$6.42 |  $-$8.39
16503 |  $-$15.24 |  $-$6.29 |  $-$8.31
17004 |  $-$15.73 |  $-$6.02 |  $-$8.19
17524 |  $-$16.12 |  $-$5.83 |  $-$8.11
18022 |  $-$16.33 |  $-$5.70 |  $-$8.04

To demonstrate the impact of 3D corrections derived from our numerical study
we recalculate accretion rates for a sample of 48 polluted DAZ white dwarfs
(Bergfors et al., 2014). As a function of effective temperature the change in
accretion rate is given by

$\Delta\dot{M}(T_{\mathrm{eff}})=\frac{1}{3}\sum_{\mathrm{Mg,Ca,Fe}}(\dot{M}^{\mathrm{3D}}_{el}-\dot{M}^{\mathrm{1D}}_{el})$
(12)

where, for simplicity, the individual change in accretion rate from each
element, $el$, considered - Mg, Ca and Fe - is weighted equally in the mean,
corresponding to an equal abundance of all three elements which is reasonably
close to solar abundance ratios.

Fig. 23 shows the DAZ and DBZ samples taken from Bergfors et al. (2014) (open
circles and blue filled circles, respectively) along with the 3D corrected DAZ
sample (black, filled circles). It can be seen that within the effective
temperature range considered in this paper (11 400 - 18 000 K), 3D corrections
bring the mean accretion rates between the DAZ and DBZ samples into closer
agreement. Furthermore, for the same region we find that the standard
deviation of the corrected sample is contained within the standard deviation
of the DBZ sample. This is shown in Fig. 24, where the mean (solid line) and
one standard deviation (filled error regions) are plotted for the original
sample from Bergfors et al. (2014) in the top panel and with the inclusion of
macroscopic mixing due to overshoot shown in the lower panel. The lack of a
correction below 11 000 K is due to the computational limitations of
simulating lower temperatures (due to their deep convection zones) and above
18 000 K we predict no correction due to overshoot, owing to the non-
convective, radiative atmospheres of these warmer objects. A small correction
due to overshoot is expected for the DBZ sample, but to a lesser extent than
for the DAZs at comparable temperatures.

Figure 23: Sample of polluted white dwarfs comprising DAZs (black, open) and
DBZs (blue) from Bergfors et al. (2014). The DAZ sample is also shown after
the inclusion of convective overshoot (black, filled). Figure 24: Accretion
rates from Bergfors et al. (2014) and Fig. 23 mean-averaged within eight
effective temperature bins. Shaded areas correspond to one standard deviation
for the overshoot-corrected DAZ (green), original DAZ (red) and DBZ (blue)
samples.

##  6 Conclusions

Using the RHD code CO5BOLD, we have presented the first macroscopic diffusion
coefficients for DA white dwarfs in the range 12 000 K $\leq
T_{\mathrm{eff}}\leq 13\,500$ K derived directly from 3D simulations. We have
shown that the phenomenon of convective overshoot is capable of mixing
material over a much greater extent than the MLT approximation, as suggested
before (Freytag et al., 1996; Tremblay et al., 2015; Kupka et al., 2018).
These direct simulations have demonstrated that macroscopic diffusion decays
with the vertical component of the velocities immediately beneath the unstable
layers and with the same quantity squared for at least $\approx$ 2.5 – 3.5
pressure scale heights beneath the unstable layers. For deeper layers, the 3D
diffusion experiments were inconclusive because of the presence of waves
amplified by boundary conditions, though the path integration method was less
susceptible to artificial mixing in the wave-dominated regions.

We have developed an approach to acquire a lower limit on the mixed mass in
the convection zone which consists of using wave-filtered convective
velocities as a proxy for macroscopic diffusion. We identify the points at
which the diffusion coefficients derived from the filtered velocity squared
profiles intersect the microscopic diffusion coefficients for magnesium,
calcium and iron, providing a conservative and elemental dependent estimate of
the increase in mixed mass.

We have found evidence to suggest that settling times for pure-hydrogen
atmosphere white dwarfs with $11\,400\leq T_{\mathrm{eff}}/[\mathrm{K}]\leq
18\,000$ are likely to increase by 1.5 – 3 orders of magnitude (see Fig. 22;
top). This has implications for predicting the period over which we are
unlikely to observe pollution variability in white dwarfs, providing the
possibility for constraining this model. The often discussed instantaneity of
accretion at warm DAZ white dwarfs is still valid compared to the Myr
timescales of accretion at DBZ white dwarfs (Wyatt et al., 2014), but we have
presented evidence that it may be slower than typically thought. Despite this
large increase in settling time, we predict that derived accretion rates may
increase by only up to an order of magnitude. In the immediate future we
intend to create a semi-analytic model for DA white dwarfs cooler than the 3D
grid presented here extends to, as well as for DB stars which should
experience similar effects from convective overshoot.

Our treatment of tracers as passive scalars contains the implicit
approximation that the particles have no mass or charge. For future work we
consider that this approximation could be lifted to provide the opportunity to
investigate the role of other instabilities which may arise from, for example,
the presence of a chemical gradient.

Ultimately the intention behind this research is to provide robust constraints
on accretion processes by improving the characterisation of the mass reservoir
in white dwarf convection zones, and should serve as an important step to
better understand the evolution of planetesimals in old planetary systems.

## Acknowledgements

The research leading to these results has received funding from the European
Research Council under the European Union’s Horizon 2020 research and
innovation programme n. 677706 (WD3D). HGL acknowledges financial support by
the Sonderforschungsbereich SFB 881 “The Milky Way System” (subprojects A4) of
the German Research Foundation (DFG). This research was supported in part by
the National Science Foundation under Grant No. NSF PHY-1748958. We are
grateful to Boris Gänsicke, Evan Bauer, Lars Bildsten and Christopher Manser
for helpful discussions in regard to the implications for planetary accretion
and stellar evolution.

## References

  * Bauer & Bildsten (2018) Bauer E. B., Bildsten L., 2018, ApJ, 859, L19
  * Bauer & Bildsten (2019) Bauer E. B., Bildsten L., 2019, ApJ, 872, 96
  * Bergfors et al. (2014) Bergfors C., Farihi J., Dufour P., Rocchetto M., 2014, MNRAS, 444, 2147
  * Böhm-Vitense (1958) Böhm-Vitense E., 1958, Z. Astrophys., 46, 108
  * Chayer et al. (1995a) Chayer P., Fontaine G., Wesemael F., 1995a, ApJS, 99, 189
  * Chayer et al. (1995b) Chayer P., Vennes S., Pradhan A. K., Thejll P., Beauchamp A., Fontaine G., Wesemael F., 1995b, ApJ, 454, 429
  * Colella & Woodward (1984) Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54, 174
  * Debes & López-Morales (2008) Debes J. H., López-Morales M., 2008, ApJ, 677, L43
  * Dufour et al. (2012) Dufour P., Kilic M., Fontaine G., Bergeron P., Melis C., Bochanski J., 2012, ApJ, 749, 6
  * Farihi et al. (2009) Farihi J., Jura M., Zuckerman B., 2009, ApJ, 694, 805
  * Farihi et al. (2012) Farihi J., Gänsicke B. T., Wyatt M. C., Girven J., Pringle J. E., King A. R., 2012, MNRAS, 424, 464
  * Farihi et al. (2013) Farihi J., Gänsicke B. T., Koester D., 2013, Science, 342, 218
  * Farihi et al. (2016) Farihi J., Koester D., Zuckerman B., Vican L., Gänsicke B. T., Smith N., Walth G., Breedt E., 2016, MNRAS, 463, 3186
  * Farihi et al. (2018) Farihi J., et al., 2018, Monthly Notices of the Royal Astronomical Society, 481, 2601
  * Freytag (2013) Freytag B., 2013, Memorie della Societa Astronomica Italiana Supplementi, 24, 26
  * Freytag (2017) Freytag B., 2017, Mem. Soc. Astron. Italiana, 88, 12
  * Freytag et al. (1996) Freytag B., Ludwig H.-G., Steffen M., 1996, A&A, 313, 497
  * Freytag et al. (2010) Freytag B., Allard F., Ludwig H.-G., Homeier D., Steffen M., 2010, A&A, 513, A19
  * Freytag et al. (2012) Freytag B., Steffen M., Ludwig H.-G., Wedemeyer-Böhm S., Schaffenberger W., Steiner O., 2012, Journal of Computational Physics, 231, 919
  * Gänsicke et al. (2012) Gänsicke B. T., Koester D., Farihi J., Girven J., Parsons S. G., Breedt E., 2012, MNRAS, 424, 333
  * Gentile Fusillo et al. (2017) Gentile Fusillo N. P., Gänsicke B. T., Farihi J., Koester D., Schreiber M. R., Pala A. F., 2017, MNRAS, 468, 971
  * Hollands et al. (2018) Hollands M. A., Gänsicke B. T., Koester D., 2018, MNRAS, 477, 93
  * Jura (2003) Jura M., 2003, ApJ, 584, L91
  * Koester (2009) Koester D., 2009, A&A, 498, 517
  * Koester et al. (2014) Koester D., Gänsicke B. T., Farihi J., 2014, A&A, 566, A34
  * Kupka et al. (2018) Kupka F., Zaussinger F., Montgomery M. H., 2018, MNRAS, 474, 4660
  * Ludwig et al. (2006) Ludwig H.-G., Allard F., Hauschildt P. H., 2006, A&A, 459, 599
  * Manser et al. (2016) Manser C. J., et al., 2016, MNRAS, 455, 4467
  * Manser et al. (2019) Manser C. J., et al., 2019, Science, 364, 66
  * Melis & Dufour (2017) Melis C., Dufour P., 2017, ApJ, 834, 1
  * Michaud (1970) Michaud G., 1970, ApJ, 160, 641
  * Paquette et al. (1986a) Paquette C., Pelletier C., Fontaine G., Michaud G., 1986a, ApJS, 61, 177
  * Paquette et al. (1986b) Paquette C., Pelletier C., Fontaine G., Michaud G., 1986b, ApJS, 61, 197
  * Pelletier et al. (1986) Pelletier C., Fontaine G., Wesemael F., Michaud G., Wegner G., 1986, ApJ, 307, 242
  * Prandtl & Tietjens (1925) Prandtl L., Tietjens O., 1925, Naturwissenschaften, 13, 1050
  * Roxburgh (1978) Roxburgh I. W., 1978, A&A, 65, 281
  * Schatzman (1945) Schatzman E., 1945, Annales d’Astrophysique, 8, 143
  * Spiegel (1963) Spiegel E. A., 1963, ApJ, 138, 216
  * Tassoul et al. (1990) Tassoul M., Fontaine G., Winget D. E., 1990, ApJS, 72, 335
  * Thoul et al. (1994) Thoul A. A., Bahcall J. N., Loeb A., 1994, ApJ, 421, 828
  * Tremblay et al. (2013a) Tremblay P.-E., Ludwig H.-G., Steffen M., Freytag B., 2013a, A&A, 552, A13
  * Tremblay et al. (2013b) Tremblay P.-E., Ludwig H.-G., Freytag B., Steffen M., Caffau E., 2013b, A&A, 557, A7
  * Tremblay et al. (2013c) Tremblay P.-E., Ludwig H.-G., Steffen M., Freytag B., 2013c, A&A, 559, A104
  * Tremblay et al. (2015) Tremblay P.-E., Ludwig H.-G., Freytag B., Fontaine G., Steffen M., Brassard P., 2015, ApJ, 799, 142
  * Tremblay et al. (2017) Tremblay P.-E., Ludwig H.-G., Freytag B., Koester D., Fontaine G., 2017, Mem. Soc. Astron. Italiana, 88, 104
  * Vanderburg et al. (2015) Vanderburg A., et al., 2015, Nature, 526, 546
  * Veras (2016) Veras D., 2016, Royal Society Open Science, 3, 150571
  * Wachlin et al. (2017) Wachlin F. C., Vauclair G., Vauclair S., Althaus L. G., 2017, A&A, 601, A13
  * Wilson et al. (2016) Wilson D. J., Gänsicke B. T., Farihi J., Koester D., 2016, MNRAS, 459, 3282
  * Wilson et al. (2019) Wilson T. G., Farihi J., Gänsicke B. T., Swan A., 2019, MNRAS,
  * Wyatt et al. (2014) Wyatt M. C., Farihi J., Pringle J. E., Bonsor A., 2014, MNRAS, 439, 3371
  * Xu et al. (2017) Xu S., Zuckerman B., Dufour P., Young E. D., Klein B., Jura M., 2017, ApJ, 836, L7
  * Xu et al. (2018) Xu S., et al., 2018, MNRAS, 474, 4795
  * Zahn (1991) Zahn J.-P., 1991, A&A, 252, 179
  * Zuckerman et al. (2003) Zuckerman B., Koester D., Reid I. N., Hünsch M., 2003, ApJ, 596, 477
  * Zuckerman et al. (2007) Zuckerman B., Koester D., Melis C., Hansen B. M., Jura M., 2007, ApJ, 671, 872
  * Zuckerman et al. (2010) Zuckerman B., Melis C., Klein B., Koester D., Jura M., 2010, ApJ, 722, 725
  * van Maanen (1917) van Maanen A., 1917, PASP, 29, 258