# The Geometry of the G29-38 White Dwarf Dust Disk from Radiative Transfer
Modeling

Nicholas P. Ballering Virginia Initiative on Cosmic Origins Fellow.
Department of Astronomy, University of Virginia, Charlottesville, VA 22904,
USA  Colette I. Levens  Department of Atmospheric, Oceanic and Planetary
Physics, University of Oxford, Clarendon Laboratory, Oxford, OX1 3PU, UK  Kate
Y. L. Su Steward Observatory, University of Arizona, 933 North Cherry Avenue,
Tucson, AZ 85721, USA  L. Ilsedore Cleeves Department of Astronomy, University
of Virginia, Charlottesville, VA 22904, USA  Department of Chemistry,
University of Virginia, Charlottesville, VA 22904, USA Nicholas P. Ballering
<EMAIL_ADDRESS>
###### Abstract

Many white dwarfs host disks of dust produced by disintegrating planetesimals
and revealed by infrared excesses. The disk around G29-38 was the first to be
discovered and is now well-observed, yet we lack a cohesive picture of its
geometry and dust properties. Here we model the G29-38 disk for the first time
using radiative transfer calculations that account for radial and vertical
temperature and optical depth gradients. We arrive at a set of models that can
match the available infrared measurements well, although they overpredict the
width of the 10 $\micron$ silicate feature. The resulting set of models has a
disk inner edge located at 92–100 $R_{\text{WD}}$ (where $R_{\text{WD}}$ is
the white dwarf radius). This is farther from the star than inferred by
previous modeling efforts due to the presence of a directly illuminated front
edge to the disk. The radial width of the disk is narrow ($\leq$10
$R_{\text{WD}}$); such a feature could be explained by inefficient spreading
or the proximity of the tidal disruption radius to the sublimation radius. The
models have a half-opening angle of $\geq$1.4∘. Such structure would be in
strong contradiction with the commonly employed flat-disk model analogous to
the rings of Saturn, and in line with the vertical structure of main-sequence
debris disks. Our results are consistent with the idea that disks are
collisionally active and continuously fed with new material, rather than
evolving passively after the disintegration of a single planetesimal.

††software: RADMC-3D (Dullemond et al., 2012)

##  1 INTRODUCTION

Dust disks around white dwarfs are evident from their infrared (IR) excess,
and dozens of such disks have been discovered (e.g. Rebassa-Mansergas et al.,
2019; Dennihy et al., 2020a; Xu et al., 2020; Lai et al., 2021). The dust
likely originates from planetesimals scattered into highly eccentric orbits
then tidally disrupted upon close approach to the white dwarf star (Debes &
Sigurdsson, 2002; Jura, 2003; Veras et al., 2014). Actively disintegrating
planetesimals have been seen via their transit signatures in a growing number
of systems (Vanderburg et al., 2015; Vanderbosch et al., 2020, 2021; Guidry et
al., 2021; Farihi et al., 2022) although the process by which the resulting
debris forms a dust disk remains an active area of investigation (e.g., Veras
et al., 2015; Brouwers et al., 2022).

Disk material is fed onto the star itself and observed as metal-rich pollution
in the white dwarf photosphere. Hundreds of polluted white dwarfs are now
known. Analysis of polluted white dwarfs allows the elemental abundance of the
original planetesimals to be inferred (e.g., Xu et al., 2019). Gas disks have
also been found around 21 white dwarfs, seen via emission or absorption lines
of atomic Ca, Fe, O, Mg, Si, and Na (e.g., Gänsicke et al., 2006; Dennihy et
al., 2020b; Melis et al., 2020; Gentile Fusillo et al., 2021). How the dust
and gas disks interact, why gas disks are less common than dust disks, and the
mechanism by which material is accreted onto the star are also areas of
ongoing investigation (e.g., Rafikov, 2011; Metzger et al., 2012; Swan et al.,
2020, 2021; Li et al., 2021). See Farihi (2016) and Veras (2021) for recent
reviews of white dwarf disks and pollution.

G29-38—the subject of this study—was the first white dwarf with a detected IR
excess (Zuckerman & Becklin, 1987; Graham et al., 1990; Tokunaga et al.,
1990), and it remains one of the brightest and most well-observed dust disks.
As a result, the system has a wealth of data from broadband photometry to
spectroscopy, covering wavelengths from the near-IR to the submillimeter.
Various models have been constructed to fit the available data. An early model
was geometrically flat, analogous to the rings of Saturn, as expected from
collisional damping of the dust’s vertical motion (Jura, 2003). This model is
optically thick, resulting in featureless IR emission that depends on the
temperature profile between inner and outer disk radii and inclination. The
flat-disk model fit the available photometry at the time, and it is still
employed for many systems with limited data (e.g. Gentile Fusillo et al.,
2021). However, Spitzer/Infrared Spectrograph (IRS) spectroscopy revealed
mineralogical features, including a prominent silicate emission feature at
$\sim$10 $\micron$, which are inconsistent with entirely optically thick
emission (Reach et al., 2005). Several other bright disks with IRS spectra
also exhibit mineralogical features (Jura et al., 2007, 2009a), suggesting
that, aside from being brighter, the G29-38 disk may not be atypical.

Other models of the G29-38 disk have accounted for an optically thin component
in order to fit the mineralogical features. Purely optically thin models were
used by Reach et al. (2005, 2009), with the latter study conducting a detailed
compositional analysis of the dust. Their model did not explicitly prescribe
the dust location or disk geometry, but for a reasonable dust distribution the
inferred mass would not be entirely optically thin, so the model was not fully
self-consistent. Reach et al. (2009) also explored moderately optically thick
models and models combining optically thick and thin components. The optically
thick component was posited as a flat disk and the thin component either a
diffuse halo or warped outer disk, following the model setup Jura et al.
(2007) used to interpret the GD 362 dust disk. Xu et al. (2018a) also used a
two-component (optically thin plus optically thick) model to fit the G29-38 IR
spectral energy distribution (SED) and interpret the observed variability.
None of these two-component models represented a complete and self-consistent
picture of the dust disk.

The goal of this paper is to forward model the dust emission around G29-38
with a cohesive disk structure similar to a protoplanetary disk or main-
sequence debris disk. IR spectra of the model are generated using a radiative
transfer code that naturally accounts for radial and vertical temperature and
optical depth gradients in the disk. We describe the G29-38 data and how we
fit them with radiative transfer models in Section 2. In Section 3 we present
the model fits, illustrate the effect of each model parameter on the simulated
spectra, and show our derived model dust temperature. In Section 4 we discuss
the implications of the disk’s grain sizes, mass, vertical structure, and
radial structure, and we compare the results with the general properties of
main-sequence debris disks. We also discuss the limitations of our models. We
summarize our results in Section 5.

##  2 Methods

In short, our methodology entails compiling photometry and spectra of the
G29-38 system from the literature and fitting these data with spectra of a
model disk. The model geometry is parameterized to resemble a protoplanetary
disk or main-sequence debris disk. We use the RADMC-3D radiative transfer tool
(Dullemond et al., 2012) to generate the model spectra. We constrain the disk
parameters by fitting a grid of models to the data.

###  2.1 Archival Data Collection

We use two sets of IR spectra from the Spitzer Space Telescope (Werner et al.,
2004) IRS (Houck et al., 2004): AOR 13828096 from 2006 (SL2, LL2, and LL1
modules covering 5.2–8.7 and 14.0–38 $\micron$) and AOR 22957568 from 2007
(SL2 and SL1 modules covering 5.2–14.5 $\micron$). These spectra have been
previously presented in the literature (Reach et al., 2009; Farihi et al.,
2014; Xu et al., 2018a). We obtain the spectra from the Combined Atlas of
Sources with Spitzer IRS Spectra (Lebouteiller et al., 2011) online database.
The IRS data have uncertainties of $\sim$5% at $\lambda$ $<$ 12 $\micron$ and
larger uncertainties at longer wavelengths. We exclude the measurements at
wavelengths longer than 30 $\micron$ because they have a very low signal-to-
noise ratio. Another IRS spectrum covering 5.2–14.5 $\micron$ was acquired in
2004 (AOR 10184192), which we do not use in our analysis. It showed a 10%
weaker excess than the 2007 spectrum, and this difference was analyzed by Xu
et al. (2018a).

We also use spectra from the IRTF/SpeX instrument (Rayner et al., 2003)
originally published by Kilic et al. (2006). The spectra were provided by M.
Kilic (private communication). We use only the data from 1.3–2.5 $\micron$, as
shorter wavelengths trace the white dwarf photosphere, not the disk. We re-
scale this spectrum by a factor of 0.96 so the short end aligns with our white
dwarf photosphere model, as was also found to be necessary by Reach et al.
(2009). We assign uncertainty values of 5% to the points in this spectrum.

Finally, we compile IR photometry from the literature (Table 1). G29-38 is
known to be variable in the optical and near-IR (Graham et al., 1990; Kleinman
et al., 1998) due to stellar pulsation. This can lead to systematic offsets
between measurements that are not taken simultaneously and likely explains why
the J- and H-band photometry are somewhat fainter than our model photosphere
(see Section 3.1). The IRAC photometry was measured over multiple epochs (Xu
et al., 2018a); we adopt the median value and use the standard deviation of
the measurements for the uncertainty. No uncertainty estimates were provided
for the IRS peakup photometry (Reach et al., 2005), so we adopt a conservative
estimate of 10%. G29-38 was not detected at 70 $\micron$ with Spitzer/MIPS
(Jura et al., 2009b), at 100 or 160 $\micron$ with Herschel/PACS, nor in the
(sub)mm with ALMA (Farihi et al., 2014). We do not formally include these
measurements in our fitting process, but we note that our models that fit the
IR detections are also consistent with these upper limits.

Table 1: Photometric Data  $\lambda$ ($\micron$) |  $F_{\nu}$ (mJy) | Instrument | Reference
---|---|---|---
1.25 | 9.00 $\pm$ 0.24 | 2MASS | Cutri et al. (2003)
1.63 | 6.15 $\pm$ 0.24 | 2MASS | Cutri et al. (2003)
2.19 | 5.55 $\pm$ 0.16 | 2MASS | Cutri et al. (2003)
3.35 | 7.68 $\pm$ 0.17 | WISE | Cutri & et al. (2012)
3.55 | 8.08 $\pm$ 0.15 | IRAC | Xu et al. (2018a)
4.49 | 8.99 $\pm$ 0.19 | IRAC | Xu et al. (2018a)
4.60 | 8.73 $\pm$ 0.17 | WISE | Cutri & et al. (2012)
5.73 | 8.27 $\pm$ 0.11 | IRAC | Xu et al. (2018a)
7.87 | 8.48 $\pm$ 0.24 | IRAC | Xu et al. (2018a)
10.50 | 11.1 $\pm$ 2.2 | IRTF | Tokunaga et al. (1990)
11.60 | 7.39 $\pm$ 0.20 | WISE | Cutri & et al. (2012)
16.0 | 3.7 $\pm$ 0.37 | IRS | Reach et al. (2005)
22.10 | 2.4 $\pm$ 1.02 | WISE | Cutri & et al. (2012)
23.68 | 2.75 $\pm$ 0.05 | MIPS | Xu et al. (2018a)

###  2.2 Disk Model

We model the dust distribution around G29-38 as an azimuthally symmetric disk
with a geometry like that of a protoplanetary disk or main-sequence debris
disk. While the structures of protoplanetary and debris disks have been
revealed by resolved imaging, this is not possible for white dwarf disks owing
to their much smaller sizes, so white dwarf disk structure is much more
uncertain. Thus, we take a simplified approach and acknowledge reality may be
more complex.

The radial distribution of our dust disk model follows a power law between
$r_{\text{in}}$ and $r_{\text{out}}=r_{\text{in}}+\Delta r$. That is, the
surface density follows

$\Sigma(r)=\Sigma_{\text{in}}\left(\frac{r}{r_{\text{in}}}\right)^{-\gamma}.$
(1)

We fix $\gamma$ = 1, a typical value found for protoplanetary disks from
resolved submillimeter images (Andrews et al., 2010). We parameterize the
amount of dust in terms of total mass, $M$, so the surface density at the
inner edge of the disk is

$\Sigma_{\text{in}}=\frac{M(2-\gamma)}{2\pi
r_{\text{out}}^{2}\left(\frac{r_{\text{out}}}{r_{\text{in}}}\right)^{-\gamma}-2\pi
r_{\text{in}}^{2}}.$ (2)

The dust distribution decreases exponentially above and below the midplane,
with volume density

$\rho(r,z)=\frac{\Sigma(r)}{\sqrt{2\pi}H(r)}\exp\left[-\frac{1}{2}\left(\frac{z}{H(r)}\right)^{2}\right].$
(3)

The disk scale height, $H$, increases linearly with $r$, so we parameterize
the vertical structure using the half-opening angle of the disk, $\Theta$,
i.e., $H(r)=r\tan(\Theta)$. In protoplanetary disks, the scale height
increases with radius faster than linearly; that is, the disks are flared.
This is due to vertical support from the large gas reservoir in hydrostatic
equilibrium. No gas is detected in the G29-38 disk, so our disk model is not
flared.

While $r$ and $z$ are cylindrical coordinates, we use a spherical coordinate
grid to set up the radiative transfer simulation, as RADMC-3D does not support
cylindrical coordinates. We use a regular grid with 100 cells in the radial
direction, spaced logarithmically between $r_{\text{in}}$ and
$\sqrt{r_{\text{out}}^{2}+(4H_{\text{out}})^{2}}$, ensuring that the full
radial extent of the disk is sampled to four scale heights. We use 116 cells
in the polar coordinate, $\theta$, with higher density within four scale
heights of the midplane. Specifically, we use nine cells from
$0<\theta<\frac{\pi}{2}-4\Theta$, 98 cells from
$\frac{\pi}{2}-4\Theta<\theta<\frac{\pi}{2}+4\Theta$, and another nine cells
from $\frac{\pi}{2}+4\Theta<\theta<\pi$. The cells are linearly spaced within
each region. We use 20 cells, equally spaced, in the azimuthal direction.

We use amorphous astronomical silicates for the dust composition, with optical
constants from Draine (2003)111Accessed from
https://www.astro.princeton.edu/~draine/dust/dust.diel.html. Xu et al. (2018a)
used either a combination of astronomical silicates and crystalline forsterite
or olivine grains, along with an opaque disk component, to produce a decent
fit to the G29-38 spectra. The crystalline component contributes more subtle
and narrow spectral features compared to the astronomical silicates. We ignore
the crystalline component, as our goal is to fit the most prominent structure
in the IR spectrum. We generate a library of opacity spectra for the dust at a
range of grain sizes using their optical constants and Mie theory. We do this
using a code included in the RADMC-3D software package based on the algorithm
by Bohren & Huffman (1983).

We model a power-law distribution of grain sizes, $n(a)\propto a^{-3.5}$
between $a_{\text{min}}$ and $a_{\text{max}}$ (Mathis et al., 1977), as is
typically found for protoplanetary and debris disks. We fix $a_{\text{max}}$ =
5 $\micron$ but leave $a_{\text{min}}$ as a free parameter. To sample the
distribution, we populate the RADMC-3D grid with 30 grain sizes, each treated
as a separate “species.” All grain sizes follow the same spatial distribution.

The disk is heated by radiation from the central white dwarf star, which we
model as a 11,240 K spherical blackbody with a radius of $R_{\text{WD}}$ =
9.04 $\times$ $10^{8}$ cm (Xu et al., 2018a). The radiative transfer
simulation proceeds in two steps. First, the dust temperature is computed
throughout the disk using $10^{7}$ photon packets sampling wavelengths from
0.08 to 1000 $\micron$. The radiative transfer is accelerated in high-optical-
depth parts of the disk with a modified random walk procedure. Scattering is
not included in the radiative transfer. Second, to compare with observations,
the thermal emission spectrum of the star plus disk is calculated using
$10^{4}$ photon packets. The spectrum is computed on a custom wavelength grid
from 0.1–1500 $\micron$ sampled finely ($\Delta\lambda=0.05\,\micron$) from
5.1–38 $\micron$ and more coarsely elsewhere. Following Xu et al. (2018a), we
fix the disk inclination to $i$ = 30∘, and we discuss the effect of
inclination in Section 4.5. We normalize the resulting spectrum to a distance
of 17.5 pc.

###  2.3 Fitting Procedure

Table 2: Model Parameters  Parameter | Description | Best Fit | Good FitsaaGood fits are models with $\chi^{2}$ $<$ $1.2$ times the minimum $\chi^{2}$
---|---|---|---
$r_{\text{in}}$ ($R_{\text{WD}}$) | Inner radius | 96 | 92–100
$\Delta r$ ($R_{\text{WD}}$) | Radial width | 1 |  $\leq$10
$\Theta$ (deg.) | Half-opening angle | 5 |  $\geq$1.4
$M$ ($10^{18}$ g) | Dust mass | 4.5 | 4–5
$a_{\text{min}}$ ($\micron$) | Minimum grain size | 1.4 | 1–2

To search for the best-fit model, we vary five parameters: $r_{\text{in}}$,
$\Delta r$, $\Theta$, $M$, and $a_{\text{min}}$. The following parameters
remain fixed: $\gamma$ = 1, $a_{\text{max}}$ = 5 $\micron$, and $i$ = 30∘. We
parameterize $r_{\text{in}}$ and $\Delta r$ in units of $R_{\text{WD}}$, which
provides a convenient scale for the size of the disk. For reference,
$100R_{\text{WD}}$ = 0.006 au.

We evaluate the fit according to the $\chi^{2}$ metric. We first search for a
good fit by exploring parameter space in an ad hoc manner. We then build a
grid of models centered around our best-fit model. The grid entails $9408$
models. To quantify the range of parameters that yield good fits, we examine
the set of models with $\chi^{2}$ values within some factor of the lowest
$\chi^{2}$. We conclude that models with $\chi^{2}$ $<$ $1.2$ times the
minimum $\chi^{2}$ achieve a good fit to the data.

When computing $\chi^{2}$, we did not distinguish between photometric and
spectroscopic data points. Thus, the spectra, due to having more data points,
naturally influence the fit more so than does the photometry. Some studies
that fit IR SEDs counteract this by increasing the weight given to photometry
when computing $\chi^{2}$ (e.g., Ballering et al., 2013). This may be
necessary when the spectra and photometry sample non-overlapping wavelengths
(e.g., mid-IR and far-IR). In this case, the spectra and photometry sample
similar wavelengths, so up-weighting the photometry was not necessary.

##  3 Results

Figure 1: Measured infrared spectral energy distribution of G29-38 along with
the 41 good-fitting models (green). For clarity, we do not show error bars on
the spectroscopic data points.

###  3.1 Best-Fit Model

The parameters of the best-fit model and good-fitting models are summarized in
Table 2. Our best-fit model has an inner edge at 96 $R_{\text{WD}}$ and
relatively narrow radial width of 1 $R_{\text{WD}}$. The half-opening angle is
5∘. The dust mass is 4.5 $\times$ $10^{18}$ g and the minimum grain size is
1.4 $\micron$.

Figure 2: Lowest $\chi^{2}$ value found in the grid for each value of each
parameter, normalized to the overall minimum $\chi^{2}$. The green line
indicates the threshold of $1.2$ times the minimum $\chi^{2}$, below which
models are considered to be good fits.

The spectra of all $41$ good-fitting models from the grid are plotted in
Figure 1. While the models do not match the data perfectly in all wavelength
ranges, they do capture the most salient features of the IR spectrum: excess
starting at $\sim$1.7 $\micron$, a featureless excess peaking at $\sim$5
$\micron$, a prominent silicate feature at $\sim$10 $\micron$, a shoulder at
$\sim$18 $\micron$, and a steep decrease in flux toward longer wavelengths.
This is the first time that all these aspects of the IR spectra have been
reproduced from a single coherent physical disk model.

The most prominent discrepancy is around the silicate feature, with the model
feature being broader than the data. This could be due to our use of Mie
theory, which assumes grains are perfect spheres. Model opacity spectra of
non-spherical grains show differences in their IR solid state features and
generally better agreement with laboratory measurements (e.g., Mutschke et
al., 2009). Another discrepancy is that the J and K photometry are fainter
than the photosphere model, likely due to variability from stellar pulsations.
This has no impact on the fit, as the disk contributes little flux at these
wavelengths.

In Figure 2 we show the 1D $\chi^{2}$ contours for each parameter—that is, the
best $\chi^{2}$ value achieved in the grid with each parameter fixed at each
value. We show the 2D version of these constraints in Figure 3, with parameter
pairs that achieve a good fit marked with an $\times$.

Figure 3: Lowest $\chi^{2}$ value found for each pair of parameter values in
the model grid. Color corresponds to the base-10 log of the $\chi^{2}$
normalized to the minimum $\chi^{2}$ value. Parameter combinations marked with
an $\times$ have good-fitting models with $\chi^{2}$ $<$ $1.2$ times the
minimum $\chi^{2}$. Note that the parameter values are not uniformly spaced.

###  3.2 Impact of Each Parameter

We illustrate the effect of each parameter on the model spectra in Figure 4.
In the left column, one parameter at a time is varied while holding all other
parameters fixed at the overall best-fit values. The right column shows models
where one parameter is fixed while the other four parameters are free to
converge to the best-fitting model. This illustrates how degeneracies among
the parameters can lead to good fits for a range of parameter values, as well
as how models outside of this range are excluded by their poor fits.

The inner edge location strongly affects the shorter wavelengths of the model.
Smaller values place the dust closer to the star, leading to a brighter excess
arising at shorter wavelengths.

Wider belts (with the same inner edge location and total mass) shift more mass
to larger radii, resulting in the model being fainter at short wavelengths.
Compensating for this in the fitting results in models that are too bright at
longer wavelengths. The smallest $\Delta r$ value in our model grid (0.1
$R_{\text{WD}}$) yields a good fit, so our constraints on $\Delta r$ are best
described as an upper limit.

The half-opening angle determines the size of the disk as seen by the star and
thus the amount of radiation intercepted. Higher/lower values
increase/decrease the flux at all wavelengths. Our fitting results in a lower
limit on $\Theta$.

The total dust mass also increases or decreases the model flux at all
wavelengths. Lower masses result in a more prominent silicate feature relative
to the continuum, as more of the disk is optically thin. In contrast, high
masses have higher optical depths, resulting in a muted silicate feature.

Smaller grain sizes lead to a brighter continuum flux, especially at short
wavelengths, because there is more grain surface area for the same total mass,
and because small grains achieve a higher temperature than larger grains in
the same radiation field. Smaller grains also result in a more prominent
silicate feature. Larger values of $a_{\text{min}}$ can be compensated for in
the fitting to some extent with smaller values of $r_{\text{in}}$, and vice
versa. However, the largest value of $a_{\text{min}}$ considered (4 $\micron$)
results in a very broad silicate feature and flux at 20–30 $\micron$ that does
not decrease with wavelength as steeply as the data.

Figure 4: Left: effect on the model spectra by varying each parameter
individually while keeping the other parameters fixed at the values of the
overall best-fit model. Right: models when a single parameter is fixed while
allowing the other four parameters to converge to their best-fitting values.
The legend shown in the right panel also applies to the models shown in the
left panel for each row.

###  3.3 Dust temperature

The radiative transfer calculation computes the dust temperature throughout
the disk for each grain size in the model. In Figure 5 we show the $z$ versus
$r$ temperature distribution of our best-fitting model for three grain sizes.
We find that the temperature decreases radially across the disk, as expected,
and that small grains are warmer than large grains, also as expected. The
maximum grain temperature of $\sim$1000 K is somewhat less than the
sublimation temperature of silicates ($\sim$1200 K), although the sublimation
temperature depends strongly on the grain chemical composition and internal
structure (Mann et al., 2004), which are not well constrained here. This
confirms that our model is physically consistent (i.e. does not require dust
hotter than the expected sublimation temperature). However, it also leaves
open the question of what, if not sublimation, sets the inner boundary of the
dust distribution.

##  4 Discussion

###  4.1 Minimum Grain Size

Our model’s minimum grain size of 1–2 $\micron$ is larger than the value of
0.1 $\micron$ used in the models of Reach et al. (2009) and Xu et al. (2018a).
Radiation pressure is very weak around white dwarf stars, so the minimum grain
size is not set by the “blowout” size as in debris disks around many main-
sequence stars.222Radiation pressure is insignificant for grains on circular
orbits, as we assume for our model. Grains on very eccentric orbits around a
white dwarf are susceptible to radiation pressure (Brouwers et al., 2022).
Kenyon & Bromley (2017) note that small grains are likely to sublimate rapidly
and predict a minimum grain size of 0.1–1 $\micron$.

The dust clouds that transit WD 1145+017 exhibit wavelength-independent
extinction out to 4.5 $\micron$, indicating the grains must be larger than 1.5
$\micron$ in size (Xu et al., 2018b). Xu et al. (2018b) argue that smaller
grains are rapidly sublimated. This is similar to the range of values found by
our set of models for the minimum grain size for the G29-38 disk. Thus, our
results are consistent with a minimum grain size set by sublimation although,
as shown in Section 3.3, our model dust temperature does not quite reach the
expected sublimation temperature for silicates.

###  4.2 Disk Mass

Prior models that assumed a fully optically thick disk are insensitive to the
dust mass, so we cannot compare our results with such models. When comparing
with optically thin models, we must be aware of the assumed maximum grain
size, as most of the mass will be in the largest grains. The optically thin
model of Reach et al. (2009) had 2 $\times$ $10^{19}$ g in grains up to
$a_{\text{max}}$ = 10 $\micron$ (1.6 $\times$ $10^{19}$ g in grains up to 5
$\micron$ for their power-law distribution). Xu et al. (2018a) found 4.8
$\times$ $10^{18}$ g in grains up to 3 $\micron$ (6.2 $\times$ $10^{18}$ g in
grains up to 5 $\micron$) in their optically thin astronomical silicate model.
This does not include the forsterite or optically thick components of their
model. Overall, our derived dust mass of 4–5 $\times$ $10^{18}$ g is somewhat
lower than found by previous models but consistent within uncertainties.

A recent measure of the X-ray flux from G29-38 derived an instantaneous
accretion rate of 1.63 $\times$ 109 g s-1 of material onto the white dwarf
(Cunningham et al., 2022). At this rate, a disk with our best-fit mass of 4.5
$\times$ $10^{18}$ g would have a lifetime of 2.8 $\times$ $10^{9}$ s (87 yr).
This lifetime is only a few times longer than the time since this disk was
discovered in 1987, and no substantial mass depletion has been observed. (Note
that this lifetime is a lower limit, as our mass measurement only includes
grains up to 5 $\micron$.) Instead, the IR excess at 10 $\micron$ was observed
to increase by 10% from 2004 to 2007, likely due to an increase in dust mass
(Xu et al., 2018a).

White dwarf disks can be replenished before they are fully depleted. This
could be from a single large planetesimal that continuously sheds mass, such
as the actively disintegrating body around WD 1145+017 (Vanderburg et al.,
2015) or the remnant planetesimal inferred around SDSS J1228+1040 (Manser et
al., 2019). Alternatively, the disk could be resupplied by a continuous infall
of small planetesimals from larger orbital distances. This picture is
supported by the widespread IR variability of white dwarf disks (Swan et al.,
2020) as well as by a statistical analysis of white dwarf accretion rates
inferred from their atmospheric pollution (Wyatt et al., 2014).

The evolution of white dwarf dust disks was modeled numerically by Kenyon &
Bromley (2017). They find that collisions play a dominant role in disk
evolution, breaking larger solids down into small dust grains, which are then
destroyed by sublimation. This differs from the traditional flat-disk model
where destructive collisions have a negligible impact on disk evolution, and
Poynting-Robertson (P-R) drag drives mass loss as grains drift inward.

In collisional models where no new material is added to the disk, the disk
mass drops by several orders of magnitude after $10^{2}$–$10^{4}$ yr, with the
IR excess falling below detectable levels after 10–100 yr (Kenyon & Bromley,
2017). This timescale does not strictly rule out this scenario for the G29-38
disk, but it would require the disk to have formed not long before it was
discovered.

The scenario of ongoing disk replenishment was also modeled by Kenyon &
Bromley (2017), who found that the disk mass reaches a stable equilibrium.
This equilibrium mass (their Equation 11) depends on the white dwarf mass,
mass input rate, size of the infalling planetesimals, material density, and
the disk radial location, width, and eccentricity. Using the X-ray-measured
accretion rate for the mass input rate, 1 km for the planetesimal size, 3 g
cm-3 for the material density, 0.01 for the disk eccentricity, and our best-
fit model for the disk location and width, we find an equilibrium mass of 1.0
$\times$ 1018 g. This mass is within a factor of $\sim$5 compared to the disk
mass of our best-fit model, which is a fairly good agreement considering the
inherent uncertainties in many inputs to this calculation.

If disks are resupplied by multiple infalling planetesimals, which are
unlikely to have identical compositions, we might expect the composition of
the white dwarf pollution to vary over time. Yet, the equivalent widths of the
metallic lines in the G29-38 atmosphere appear stable (Debes & López-Morales,
2008). Similarly, the pollution of the white dwarf GD 56 was observed to be
stable over a 15 yr timescale despite its IR excess varying by $\sim$20%
(Farihi et al., 2018). This suggests disk and accretion processes are
decoupled to some degree. A detailed investigation of the implications of our
disk models on the accretion process is beyond the scope of this paper.

Figure 5: Dust temperature distribution of the best-fitting model at three
different grain sizes. The pixels in these images, which are linearly spaced
in cylindrical coordinate, do not correspond to radiative transfer model
cells, which are non-linearly spaced and in spherical coordinates. We plot the
temperature up to $z$ = $\pm 3H_{\text{out}}$. Note that the axis dimensions
are not to scale.

###  4.3 Vertical Structure

Prior models assumed the disk was vertically flat (and optically thick) or
optically thin with the vertical structure unconstrained, leaving little
ground for quantitative comparison. The scale height at the inner edge of our
model is $H_{\text{in}}$ $\geq$ 2.2 $R_{\text{WD}}$. This is larger than the
stellar radius, so the disk is primarily illuminated at its front edge, rather
than at its upper and lower surfaces. Our model thus corresponds to the
“physically thick” regime discussed by Reach et al. (2009). In this regime,
the half-opening angle must be greater than 0.8∘ in order to intercept enough
stellar radiation to account for the total disk luminosity (Reach et al.,
2009). Our models, with $\Theta$ $\geq$1.4∘, are consistent with this
constraint.

Flat disk models were motivated by the idea that grains should dynamically
settle (Jura, 2003). The disk evolution model of Kenyon & Bromley (2017) found
that the collisional timescale is shorter than the settling timescale. That
is, planetesimals are broken down into dust and lost via sublimation faster
than those particles can settle. Our set of models are consistent with the
idea that significant vertical settling has not taken place, which would
support the conclusion that collisions play a major role in white dwarf disk
evolution.

Ongoing infall of small planetesimals onto the disk, as discussed in Section
4.2, could maintain its vertically extended structure. In the WD 1145+017
system, clouds of dust released from actively disintegrating planetesimals
transit in front of the star, but the existing disk of dust does not (Xu et
al., 2018b). This demonstrates that new material can enter a disk at a
different inclination and thus contribute to the disk’s scale height.

With collisions playing a dominant role in the G29-38 disk, we can compare its
vertical structure with that of other systems dominated by collisional
evolution: main-sequence debris disks. It is most salient to compare with
observations of debris disks in the (sub)millimeter, which trace the large
grains that are less affected by radiation pressure and gas in these disks
(Olofsson et al., 2022).

The edge-on Kuiper belt-like debris disk around $\beta$ Pictoris has dust in
two overlapping vertical distributions: a dynamically cold component with a
half-opening angle of $\Theta$ = 0.8∘ and a dynamically hot component with
$\Theta$ = 6.3∘, as measured by ALMA (Matrà et al., 2019). The AU Mic debris
disk was imaged with ALMA to have a half-opening angle of $\Theta$ = 1.1–1.8∘,
depending on the wavelength and specific analysis (Daley et al., 2019; Vizgan
et al., 2022). Similar analyses found $\Theta$ = 2.3∘ for the disk around HR
4796 (Kennedy et al., 2018) and $\Theta$ = 2.9∘ for the disk around q1 Eri
(Lovell et al., 2021). These values are consistent with (or slightly lower
than) the half-opening angle values our models obtain for the G29-38 disk,
which would suggest this system is more similar to main-sequence debris disks
than to the rings of Saturn.

###  4.4 Radial Location and Width

The location of the disk inner edge in our model (92–100 $R_{\text{WD}}$) is
significantly farther from the star than in previous models, and the radial
width ($\leq$10 $R_{\text{WD}}$) is smaller. For comparison, the flat opaque
disk model of Jura (2003) had $r_{\text{in}}$ = 11 $R_{\text{WD}}$ and
$r_{\text{out}}$ = 31–69 $R_{\text{WD}}$, and von Hippel et al. (2007) found
similar values. The optically thick component from the model by Xu et al.
(2018a) spans 8–25 $R_{\text{WD}}$, while their optically thin component spans
8–834 $R_{\text{WD}}$. The best fitting “moderately optically thin” model of
Reach et al. (2009) had its inner edge at 50 $R_{\text{WD}}$. The more distant
inner edge inferred from our model is likely due to its vertical height. That
is, our model has a front edge that can be directly illuminated and heated by
the star, whereas flat disk models are only illuminated on their top and
bottom surfaces at a shallower angle. von Hippel et al. (2007) noted that the
absence of a directly illuminated front edge in the flat disk model would lead
to an underestimate of the inner edge location.

While no gas is detected in G29-38, the larger $r_{\text{in}}$ compared to
flat disk models has implications for the gas versus dust location in systems
that do possess both. In such systems, the dust disk location has been
inferred from optically thick models, and the gas and dust have been found to
be radially co-located (Brinkworth et al., 2009; Melis et al., 2010). This
finding is inconsistent with the hypothesis that the gas disk originates at
the inner edge of the dust disk via sublimation and extends inwards from there
to the stellar surface. Rather, it supports the idea that gas is produced
throughout the extent of the dust disk by sublimation or grain sputtering
(Jura, 2008). Revisiting such systems by fitting their dust disks with
radiative transfer models, as we have done for G29-38, may change the picture
of the relative locations of gas and dust and perhaps alter our understanding
of the origin of white dwarf gas disks. Unfortunately, at the moment, none of
the white dwarfs hosting gas disks have mid-IR spectra of their dust disks,
which are essential for fitting with radiative transfer models.

What might cause the G29-38 dust disk to be radially narrow? One possibility
is that the tidal disruption radius (the orbital location within which
planetesimals are destroyed by tidal forces), which may set the disk outer
edge, and the sublimation radius, which may set the inner edge, happen to be
close together for this system. The tidal disruption radius for G29-38 is
$\sim$115 $R_{\text{WD}}$ (Jura, 2003), which is near to the dust location of
our model. As discussed previously, the temperature of the small grains at the
inner edge is close to, but somewhat below, the expected sublimation
temperature. Another possibility is that the dust is produced only in a narrow
ring, perhaps just inside of the tidal disruption radius, and fails to spread
radially. P-R drag is the primary mechanism by which dust grains could drift
inwards, but in white dwarf disks the P-R drag timescale is longer than the
collision timescale (Farihi et al., 2008; Kenyon & Bromley, 2017; Farihi et
al., 2018). Thus, grains are destroyed by collisions (and sublimation) before
they can spread radially. The observed variability in the G29-38 excess can be
explained by a change in dust mass without varying its location (Xu et al.,
2018a). Similarly, the variability in the GD 56 excess is color-independent
(Farihi et al., 2018), implying the dust location is constant. This suggests
that while some aspects of a white dwarf dust disk are variable, there is some
process that stabilizes its location.

###  4.5 Model Limitations

Our model assumes an axisymmetric disk with a grain size distribution
independent of location. Theoretical work suggests that disks around white
dwarfs can be eccentric, especially during their initial formation from a
tidally disrupted asteroid (Nixon et al., 2020; Brouwers et al., 2022). While
our results show that the IR excess can be reproduced with a circular disk,
they do not rule out disks with inherently different geometries. We suggest
that radiative transfer be used to calculate the IR spectra of models derived
from detailed simulations of disk formation and evolution.

In our fitting, we fixed $\gamma$ = 1. Because a radially narrow disk
configuration is favored, $\gamma$ does not have a large effect on the model.
We also fixed $i$ = 30∘ in the fitting. Exploring the effect of inclination on
our best-fit model, we find the SED to be significantly fainter for $i$ =
0∘–20∘ because the illuminated front edge of the disk is less visible when
viewed top-down. For $i$ = 30∘–80∘ (the most likely range of inclinations),
the brightness increases somewhat with increasing inclination, but the effect
is not large. At $i$ = 90∘, the disk partially occults itself and the central
star, so the disk is fainter, especially at shorter wavelengths.

Lastly, we fixed $a_{\text{max}}$ = 5 $\micron$. To study this assumption, we
make a small number of additional models with $a_{\text{max}}$ = 10 $\micron$.
This larger maximum grain size leads to a somewhat fainter excess, as more of
the mass now resides in larger grains. We are able to recover a decent fit to
the data by also increasing $M$ to 5$\times$1018 g, decreasing $r_{\text{in}}$
to 92 $R_{\text{WD}}$, and decreasing $a_{\text{min}}$ to 1 $\micron$. These
values fall within the range of good fits from our original grid, so we
conclude that our choice of $a_{\text{max}}$ does not significantly bias the
fit.

##  5 Summary

We fit the IR measurements of the G29-38 white dwarf dust disk for the first
time using radiative transfer models based on a coherent disk structure. Prior
modeling efforts had employed an optically thick (vertically flat) model, a
purely optically thin model, or a linear combination of the two. We explore a
grid of models with five free parameters: the inner edge location, the radial
width, the vertical half-opening angle, the total dust mass, and the minimum
grain size.

We find that the disk mass (4–5 $\times$ 1018 g) is generally consistent with
(or somewhat lower than) prior models. The minimum grain size (1–2 $\micron$)
is somewhat larger than in prior models, but consistent with theoretical
expectations for dust sublimation. The half-opening angle ($\geq$1.4∘) is
inconsistent with the flat disk model and in line with the vertical structure
seen in debris disks around main-sequence stars. The inner edge location
(92–100 $R_{\text{WD}}$) is significantly farther from the star than found by
prior models because those models did not include a directly illuminated inner
edge. The radial width ($\leq$10 $R_{\text{WD}}$) is narrower than found by
prior models, perhaps due to inefficient spreading or the proximity of the
tidal disruption radius to the sublimation radius.

Overall, our models would be consistent with the growing picture of white
dwarf disks as collisionally active and continuously fed by fresh material.
The smallest grains are ultimately destroyed by sublimation. Short collisional
timescales can explain why the dust does not settle vertically or spread
radially. Thus, within our model framework, white dwarf disks would be similar
to main-sequence debris disks.

We thank the referee for providing helpful feedback on this paper. We thank
Dr. Mukremin Kilic for sharing with us the SpeX data. We thank Siyi Xu, Amy
Bonsor, Erik Dennihy, and Johan Olofsson for useful feedback on this project.
We also gratefully acknowledge the use of the University of Virginia Rivanna
High-Performance Computing system, which was utilized to run the radiative
transfer models presented in this work. N.P.B and K.Y.L.S acknowledge support
from NASA XRP grant 80NSSC22K0234. N.P.B. and C.I.L. acknowledge support from
the Virginia Initiative on Cosmic Origins (VICO).

## References

  * Andrews et al. (2010) Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241, doi: 10.1088/0004-637X/723/2/1241
  * Ballering et al. (2013) Ballering, N. P., Rieke, G. H., Su, K. Y. L., & Montiel, E. 2013, ApJ, 775, 55, doi: 10.1088/0004-637X/775/1/55
  * Bohren & Huffman (1983) Bohren, C. F., & Huffman, D. R. 1983, Absorption and scattering of light by small particles (New York: Wiley)
  * Brinkworth et al. (2009) Brinkworth, C. S., Gänsicke, B. T., Marsh, T. R., Hoard, D. W., & Tappert, C. 2009, ApJ, 696, 1402, doi: 10.1088/0004-637X/696/2/1402
  * Brouwers et al. (2022) Brouwers, M. G., Bonsor, A., & Malamud, U. 2022, MNRAS, 509, 2404, doi: 10.1093/mnras/stab3009
  * Cunningham et al. (2022) Cunningham, T., Wheatley, P. J., Tremblay, P.-E., et al. 2022, Nature, 602, 219, doi: 10.1038/s41586-021-04300-w
  * Cutri & et al. (2012) Cutri, R. M., & et al. 2012, VizieR Online Data Catalog, II/311
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, 2246, 0
  * Daley et al. (2019) Daley, C., Hughes, A. M., Carter, E. S., et al. 2019, ApJ, 875, 87, doi: 10.3847/1538-4357/ab1074
  * Debes & López-Morales (2008) Debes, J. H., & López-Morales, M. 2008, ApJ, 677, L43, doi: 10.1086/587550
  * Debes & Sigurdsson (2002) Debes, J. H., & Sigurdsson, S. 2002, ApJ, 572, 556, doi: 10.1086/340291
  * Dennihy et al. (2020a) Dennihy, E., Farihi, J., Gentile Fusillo, N. P., & Debes, J. H. 2020a, ApJ, 891, 97, doi: 10.3847/1538-4357/ab7249
  * Dennihy et al. (2020b) Dennihy, E., Xu, S., Lai, S., et al. 2020b, ApJ, 905, 5, doi: 10.3847/1538-4357/abc339
  * Draine (2003) Draine, B. T. 2003, ApJ, 598, 1026, doi: 10.1086/379123
  * Dullemond et al. (2012) Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library.  http://ascl.net/1202.015
  * Farihi (2016) Farihi, J. 2016, New A Rev., 71, 9, doi: 10.1016/j.newar.2016.03.001
  * Farihi et al. (2014) Farihi, J., Wyatt, M. C., Greaves, J. S., et al. 2014, MNRAS, 444, 1821, doi: 10.1093/mnras/stu1545
  * Farihi et al. (2008) Farihi, J., Zuckerman, B., & Becklin, E. E. 2008, ApJ, 674, 431, doi: 10.1086/521715
  * Farihi et al. (2018) Farihi, J., van Lieshout, R., Cauley, P. W., et al. 2018, MNRAS, 481, 2601, doi: 10.1093/mnras/sty2331
  * Farihi et al. (2022) Farihi, J., Hermes, J. J., Marsh, T. R., et al. 2022, MNRAS, 511, 1647, doi: 10.1093/mnras/stab3475
  * Gänsicke et al. (2006) Gänsicke, B. T., Marsh, T. R., Southworth, J., & Rebassa-Mansergas, A. 2006, Science, 314, 1908, doi: 10.1126/science.1135033
  * Gentile Fusillo et al. (2021) Gentile Fusillo, N. P., Manser, C. J., Gänsicke, B. T., et al. 2021, MNRAS, 504, 2707, doi: 10.1093/mnras/stab992
  * Graham et al. (1990) Graham, J. R., Matthews, K., Neugebauer, G., & Soifer, B. T. 1990, ApJ, 357, 216, doi: 10.1086/168907
  * Guidry et al. (2021) Guidry, J. A., Vanderbosch, Z. P., Hermes, J. J., et al. 2021, ApJ, 912, 125, doi: 10.3847/1538-4357/abee68
  * Houck et al. (2004) Houck, J. R., Roellig, T. L., van Cleve, J., et al. 2004, ApJS, 154, 18, doi: 10.1086/423134
  * Jura (2003) Jura, M. 2003, ApJ, 584, L91, doi: 10.1086/374036
  * Jura (2008) —. 2008, AJ, 135, 1785, doi: 10.1088/0004-6256/135/5/1785
  * Jura et al. (2009a) Jura, M., Farihi, J., & Zuckerman, B. 2009a, AJ, 137, 3191, doi: 10.1088/0004-6256/137/2/3191
  * Jura et al. (2007) Jura, M., Farihi, J., Zuckerman, B., & Becklin, E. E. 2007, AJ, 133, 1927, doi: 10.1086/512734
  * Jura et al. (2009b) Jura, M., Muno, M. P., Farihi, J., & Zuckerman, B. 2009b, ApJ, 699, 1473, doi: 10.1088/0004-637X/699/2/1473
  * Kennedy et al. (2018) Kennedy, G. M., Marino, S., Matrà, L., et al. 2018, MNRAS, 475, 4924, doi: 10.1093/mnras/sty135
  * Kenyon & Bromley (2017) Kenyon, S. J., & Bromley, B. C. 2017, ApJ, 844, 116, doi: 10.3847/1538-4357/aa7b85
  * Kilic et al. (2006) Kilic, M., von Hippel, T., Leggett, S. K., & Winget, D. E. 2006, ApJ, 646, 474, doi: 10.1086/504682
  * Kleinman et al. (1998) Kleinman, S. J., Nather, R. E., Winget, D. E., et al. 1998, ApJ, 495, 424, doi: 10.1086/305259
  * Lai et al. (2021) Lai, S., Dennihy, E., Xu, S., et al. 2021, ApJ, 920, 156, doi: 10.3847/1538-4357/ac1354
  * Lebouteiller et al. (2011) Lebouteiller, V., Barry, D. J., Spoon, H. W. W., et al. 2011, ApJS, 196, 8, doi: 10.1088/0067-0049/196/1/8
  * Li et al. (2021) Li, D., Mustill, A. J., & Davies, M. B. 2021, MNRAS, 508, 5671, doi: 10.1093/mnras/stab2949
  * Lovell et al. (2021) Lovell, J. B., Marino, S., Wyatt, M. C., et al. 2021, MNRAS, 506, 1978, doi: 10.1093/mnras/stab1678
  * Mann et al. (2004) Mann, I., Kimura, H., Biesecker, D. A., et al. 2004, Space Sci. Rev., 110, 269, doi: 10.1023/B:SPAC.0000023440.82735.ba
  * Manser et al. (2019) Manser, C. J., Gänsicke, B. T., Eggl, S., et al. 2019, Science, 364, 66, doi: 10.1126/science.aat5330
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425, doi: 10.1086/155591
  * Matrà et al. (2019) Matrà, L., Wyatt, M. C., Wilner, D. J., et al. 2019, AJ, 157, 135, doi: 10.3847/1538-3881/ab06c0
  * Melis et al. (2010) Melis, C., Jura, M., Albert, L., Klein, B., & Zuckerman, B. 2010, ApJ, 722, 1078, doi: 10.1088/0004-637X/722/2/1078
  * Melis et al. (2020) Melis, C., Klein, B., Doyle, A. E., et al. 2020, ApJ, 905, 56, doi: 10.3847/1538-4357/abbdfa
  * Metzger et al. (2012) Metzger, B. D., Rafikov, R. R., & Bochkarev, K. V. 2012, MNRAS, 423, 505, doi: 10.1111/j.1365-2966.2012.20895.x
  * Mutschke et al. (2009) Mutschke, H., Min, M., & Tamanai, A. 2009, A&A, 504, 875, doi: 10.1051/0004-6361/200912267
  * Nixon et al. (2020) Nixon, C. J., Pringle, J. E., Coughlin, E. R., Swan, A., & Farihi, J. 2020, arXiv e-prints, arXiv:2006.07639.  https://arxiv.org/abs/2006.07639
  * Olofsson et al. (2022) Olofsson, J., Thébault, P., Kral, Q., et al. 2022, MNRAS, 513, 713, doi: 10.1093/mnras/stac455
  * Rafikov (2011) Rafikov, R. R. 2011, MNRAS, 416, L55, doi: 10.1111/j.1745-3933.2011.01096.x
  * Rayner et al. (2003) Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362, doi: 10.1086/367745
  * Reach et al. (2005) Reach, W. T., Kuchner, M. J., von Hippel, T., et al. 2005, ApJ, 635, L161, doi: 10.1086/499561
  * Reach et al. (2009) Reach, W. T., Lisse, C., von Hippel, T., & Mullally, F. 2009, ApJ, 693, 697, doi: 10.1088/0004-637X/693/1/697
  * Rebassa-Mansergas et al. (2019) Rebassa-Mansergas, A., Solano, E., Xu, S., et al. 2019, MNRAS, 489, 3990, doi: 10.1093/mnras/stz2423
  * Swan et al. (2020) Swan, A., Farihi, J., Wilson, T. G., & Parsons, S. G. 2020, MNRAS, 496, 5233, doi: 10.1093/mnras/staa1688
  * Swan et al. (2021) Swan, A., Kenyon, S. J., Farihi, J., et al. 2021, MNRAS, 506, 432, doi: 10.1093/mnras/stab1738
  * Tokunaga et al. (1990) Tokunaga, A. T., Becklin, E. E., & Zuckerman, B. 1990, ApJ, 358, L21, doi: 10.1086/185770
  * Vanderbosch et al. (2020) Vanderbosch, Z., Hermes, J. J., Dennihy, E., et al. 2020, ApJ, 897, 171, doi: 10.3847/1538-4357/ab9649
  * Vanderbosch et al. (2021) Vanderbosch, Z. P., Rappaport, S., Guidry, J. A., et al. 2021, ApJ, 917, 41, doi: 10.3847/1538-4357/ac0822
  * Vanderburg et al. (2015) Vanderburg, A., Johnson, J. A., Rappaport, S., et al. 2015, Nature, 526, 546, doi: 10.1038/nature15527
  * Veras (2021) Veras, D. 2021, in Oxford Research Encyclopedia of Planetary Science, 1, doi: 10.1093/acrefore/9780190647926.013.238
  * Veras et al. (2014) Veras, D., Leinhardt, Z. M., Bonsor, A., & Gänsicke, B. T. 2014, MNRAS, 445, 2244, doi: 10.1093/mnras/stu1871
  * Veras et al. (2015) Veras, D., Leinhardt, Z. M., Eggl, S., & Gänsicke, B. T. 2015, MNRAS, 451, 3453, doi: 10.1093/mnras/stv1195
  * Vizgan et al. (2022) Vizgan, D., Hughes, A. M., Carter, E. S., et al. 2022, ApJ, 935, 131, doi: 10.3847/1538-4357/ac80b8
  * von Hippel et al. (2007) von Hippel, T., Kuchner, M. J., Kilic, M., Mullally, F., & Reach, W. T. 2007, ApJ, 662, 544, doi: 10.1086/518108
  * Werner et al. (2004) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1, doi: 10.1086/422992
  * Wyatt et al. (2014) Wyatt, M. C., Farihi, J., Pringle, J. E., & Bonsor, A. 2014, MNRAS, 439, 3371, doi: 10.1093/mnras/stu183
  * Xu et al. (2019) Xu, S., Dufour, P., Klein, B., et al. 2019, AJ, 158, 242, doi: 10.3847/1538-3881/ab4cee
  * Xu et al. (2020) Xu, S., Lai, S., & Dennihy, E. 2020, ApJ, 902, 127, doi: 10.3847/1538-4357/abb3fc
  * Xu et al. (2018a) Xu, S., Su, K. Y. L., Rogers, L. K., et al. 2018a, ApJ, 866, 108, doi: 10.3847/1538-4357/aadcfe
  * Xu et al. (2018b) Xu, S., Rappaport, S., van Lieshout, R., et al. 2018b, MNRAS, 474, 4795, doi: 10.1093/mnras/stx3023
  * Zuckerman & Becklin (1987) Zuckerman, B., & Becklin, E. E. 1987, Nature, 330, 138, doi: 10.1038/330138a0