The Morphology of the Asteroidal Dust around White Dwarf Stars: Optical and Near-infrared Pulsations in G29-38

More than 36 yr have passed since the discovery of the infrared excess from circumstellar dust orbiting the white dwarf G29-38, which at 17.5 pc it is the nearest and brightest of its class. The precise morphology of the orbiting dust remains only marginally constrained by existing data, subject to model-dependent inferences, and thus fundamental questions of its dynamical origin and evolution persist. This study presents a means to constrain the geometric distribution of the emitting dust using stellar pulsations measured at optical wavelengths as a variable illumination source of the dust, which reradiates primarily in the infrared. By combining optical photometry from the Whole Earth Telescope with 0.7 – 2.5 μ m spectroscopy obtained with SpeX at NASA ’ s Infrared Telescope Facility, we detect luminosity variations at all observed wavelengths, with variations at most wavelengths corresponding to the behavior of the pulsating stellar photosphere, but toward the longest wavelengths the light curves probe the corresponding time variability of the circumstellar dust. In addition to developing methodology, we ﬁ nd the pulsation amplitudes decrease with increasing wavelength for principal pulsation modes, yet increase beyond ≈ 2 μ m for nonlinear combination frequencies. We interpret these results as combination modes derived from the principal modes of identical ℓ values and discuss the implications for the morphology of the warm dust. We also draw attention to some discrepancies between our ﬁ ndings and theoretical expectations for the results of the nonlinearity imposed by the surface convection zone on mode – mode interactions and on the behavior of the ﬁ rst harmonic of the highest-amplitude pulsation mode.


Introduction
White dwarf stars are known to host circumstellar dust that is widely interpreted as tidally disrupted, minor rocky bodies (Jura & Young 2014;Farihi 2016;Veras 2016;Guidry et al. 2021).The standard scenario requires that one or more minor bodies are gravitationally perturbed by major planets (Frewen & Hansen 2014;Petrovich & Muñoz 2017;Smallwood et al. 2018), exciting eccentricities and eventually resulting in nearly radial orbits.If a rocky body passes within ≈1 R e of the white dwarf, it will be tidally shredded, leaving orbiting dust that then has a large cross section to self-collisions and Poynting-Robertson drag (Bochkarev & Rafikov 2011;Veras et al. 2015;Malamud & Perets 2020).
In the limit where the circumstellar debris is the result of just one disrupted body, this dust is expected to be in a plane and its orbit is assumed to circularize eventually, although this process is not yet fully understood (Nixon et al. 2020;Malamud et al. 2021).If the dust is derived from multiple minor bodies originating on independent orbits, then there is no a priori reason why the circumstellar material should be confined to a single plane.Therefore, constraining the dust morphology Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
should constrain its origin, and thereby the dynamical evolution of the contributing bodies.While the typical dust emission detected at 3-5 μm implies orbital radii within or around the Roche limit of the star, currently there are only modeldependent inferences, and no compelling empirical constraints on the precise distribution and geometry of the circumstellar material (Bonsor et al. 2017;Swan et al. 2020).
Giclas 29-38 (G29-38 = ZZ Psc = WD 2326+049) is both a well-known luminosity variable (Shulov & Kopatskaya 1974) of the ZZ Ceti class and the prototype white dwarf hosting circumstellar dust.Zuckerman & Becklin (1987) discovered its infrared excess, and subsequently Graham et al. (1990) and Patterson et al. (1991) compared pulsations observed in the B, J, K (and L in the latter case) photometric bands, both finding consistency with a relatively planar disk geometry.
The idea behind using stellar pulsations in this manner, which forms the basis of our updated effort, is that geometric oscillations are seen differently by Earth-based observers than by the circumstellar dust, and that these two views may, in principle, be disentangled using multiwavelength photometry.Luminosity variations arising from the photosphere are prominent primarily in the optical bands where the stellar emission peaks (Brassard et al. 1995), whereas any dust illuminated by the stellar pulsations will reradiate effectively only at infrared wavelengths.
The infrared instruments of these earlier observations recorded only weak pulsation signals, and conclusions regarding the pulsation modes were not possible, thus limiting the interpretation of the circumstellar dust morphology.Spitzer infrared observations utilizing all three instruments, including photometry and spectroscopy, were used to characterize the overall infrared dust continuum, and identify broad 10 μm emission from small silicate dust grains (dominated by olivines), but ultimately the data are consistent with a range of possible circumstellar dust configurations (Reach et al. 2005(Reach et al. , 2009;;Ballering et al. 2022).
Current constraints on the geometrical distribution of dust come purely from modeling the dust component of the spectral energy distribution (e.g., Jura 2003;Reach et al. 2005;Ballering et al. 2022).There are model-dependent inferences from luminosity and calcium equivalent width amplitudes for two large pulsations excited in 2008 that suggest polar accretion (Thompson et al. 2010).And while no stellar magnetic field has been detected toward G29-38 (1σ = 0.5 kG; Bagnulo & Landstreet 2021), a field as modest as 0.01 kG is more than sufficient to truncate a high-rate accretion flow of 10 9 g s −1 , and result in accretion along field lines (Farihi et al. 2018;Cunningham et al. 2021).The stellar optical polarization has been sensitively measured at 275.3 ± 31.9 ppm (Cotton et al. 2020), and in the context of a dust disk model, is consistent either with a nearly face-on, optically thin disk or dust with a low Bond albedo.While highly valuable, these studies do not offer further insight into the circumstellar dust geometry, beyond the dust emission modeling from infrared observations.
The approach taken in this work is an updated version of the multiwavelength pulsation studies (Graham et al. 1990;Patterson et al. 1991), and leverages the fact that infrared detectors and instrumentation have improved dramatically in the last four decades.The primary instrument of this study, SpeX, saw first light in 2000 (Rayner et al. 2003), and its infrared array was upgraded in 2014.In addition to technological advances in sensitivity and readout compared to prior decades, a key advantage of SpeX is the synchronous multiwavelength coverage possible with its low-resolution spectroscopic mode.Simultaneous observations of photosphere-and dust-emitted light obviates the need for phasing detector clocks and readouts.Concurrent, multiwavelength data are particularly valuable for G29-38, as its power spectrum is not stable and can change dramatically from one observing season to another (Kleinman et al. 1998).These significant improvements in observing capability allow what is essentially a new technique for studying the circumstellar dust morphology of G29-38.In this paper, we describe that technique and although we are unable to place firm constraints on the circumstellar dust morphology, we are able to connect the photospheric pulsations and dust responses to viable dust geometries.In addition, we find that some pulsation properties appear to be discrepant with theoretical expectations as a function of wavelength, with relevance to all ZZ Ceti white dwarf stars.

SpeX Infrared Spectroscopy
We observed G29-38 at the NASA Infrared Telescope Facility (IRTF), using the near-infrared imager and spectrograph SpeX (Rayner et al. 2003), during two runs consisting of three contiguous partial nights each, on 2018 September 4-6 and again on 2020 November 8-10.SpeX was operated in spectroscopic prism mode, simultaneously covering approximately 0.7-2.5 μm.In 2018, the standard slit was opened to 3″ width (R ≈ 20), and in 2020 this was changed to 1 6 (R ≈ 40).The sky conditions for the 2018 observing run were variable with light cirrus on two nights, and the seeing ranged from 0 7 to 1 1.In 2020, the conditions were mostly excellent, with seeing between 0 3 and 0 8.In both cases, the wide slits captured the majority of the stellar light, minimizing atmospheric dispersion losses that could alter the true spectral slope.
Integrations on the array were taken continuously during two or three observing blocks each night, preceded by, interspersed with, and immediately followed by calibration frames and standard star observations.The science target was nodded along the slit in the standard ABBA pattern for point sources, for total durations of 45 minutes to 3 hr, typically 1-2 hr.Integration times were 14.8 s in 2018 and 5.10 s in 2020, where the total time on-source was approximately 2.2 to 5.6 hr per night.Calibration frames consisted of arc lamps for wavelength calibration, flat-field images, and telluric standard stars.A total of 1671 and 2200 spectra of G29-38 were obtained in 2018 and 2020, respectively.
The data were reduced with SpeXtool version v4.1 (Vacca et al. 2003;Cushing et al. 2004).All reduction steps were carried out in the standard manner, including flat-fielding, AB pair subtraction, spectral extraction, and wavelength calibration.The science target frames were all extracted individually, while those of the telluric standards were combined into a single spectrum of higher signal-to-noise ratio (S/N).The individual spectra of G29-38 were telluric corrected on a nightly basis, in a uniform manner using a single set of standard observations corresponding most closely to the median airmass of the science target.This choice was motivated by the fact that time-varying signals would be introduced by using several telluric standards (at different airmasses) within a single night.Instead, a single telluric standard was used to correct each nightly data set, and while this does not produce the best correction of telluric features for individual frames, it introduces a minimally time-variable (and well-characterized) signal corresponding to the gradual change in airmass.
The advantage of deriving frequency measurements from SpeX spectroscopy, obtained with a wide slit, is that such measurements are differential both in wavelength and time, and therefore insensitive to variable observing conditions such as seeing or transparency changes.After the individual spectra were reduced, which corrects for nonlinearities in the detector, removes night sky lines, etc., we created a simple procedure to interpolate across thorium radiation events and cosmic rays in the detector.Only 51 spectra, or slightly over 1%, were affected by such artifacts.During this stage, spectra prior to MJD = 58365.48(the beginning of 2018 September 4) were removed owing to clouds, together with some spectra from 2020 November 9 that were significantly fainter, likely due to guiding issues.This procedure resulted in 1576 and 2184 spectra of G29-38 from 2018 and 2020, respectively.
We defined 10 bandpasses for synthetic photometry, along with two further bandpasses dominated by strong telluric lines, which were not analyzed further.Table 1 lists the 10 defined bandpasses, their central wavelength designation, and the actual S/N range in each of these resulting from synthetic photometry (the bandpasses were designed to achieve S/ N ≈ 100 per bandpass per observation).The synthetic photometry S/N values are typically higher per observation for the 2018 spectra because these integrations were three times longer than for the 2020 spectra.We discuss the subsequent normalization steps in the time-series analyses of these data below.

Moris Optical Photometry
In order to extend the wavelength coverage of the IRTF observations, the SpeX observations used the Moris CCD imager for guiding purposes (Gulbis et al. 2011).The Moris images were taken through the Sloan Digital Sky Survey gband filter, over a 60 × 60 arcsec 2 field of view on an 512 × 512 pixel 2 CCD with 0 12 pixels.These data were obtained simultaneously with the SpeX data, but no effort was made to phase the readout times of these two instruments, as the guider was necessarily obtaining shorter exposures of 3.0 s than employed with SpeX.Due to their limited temporal coverage, these data offered little additional help in identifying the stellar pulsation modes, and we did not continue this approach during the 2020 IRTF run.

WET Optical Photometry
In both 2018 and 2020, we organized contemporaneous international, ground-based optical photometry campaigns to support the IRTF observations.The telescopes, detectors, filters, and time-series durations are listed in Tables 2 and 3, where coverage totaled 112.2 hr in 2018, and 173.3 hr in 2020, and data reduction followed the prescription outlined in Provencal et al. (2012).In brief, raw images were calibrated and aperture photometry performed using the MAESTRO photometry pipeline (Dalessio 2010), where each image was corrected for bias and thermal background, and normalized by its flat field.MAESTRO performed photometry for a range of aperture sizes for the target and comparison stars.The combination of aperture and comparison star(s) that resulted in the highest quality raw light curve was chosen for each individual run from each observing site.
The next data reduction step employed the WQED pipeline (Thompson & Mullally 2009).WQED examined individual light curves for photometric quality, removed outlying points, divided by appropriate comparison stars, and corrected for   (2009).The final product is a compiled optical light curve for G29-38 spanning the 2018 and 2020 SpeX observations.

Frequency Identification
PERIOD04 (Lenz & Breger 2005) and PYRIOD (Bell 2021) were used to derive discrete Fourier transforms (FTs) of the time-series photometry, with amplitudes in units of millimodulation amplitude (1 mma = 0.1%).Figure 1 presents the FTs of the 2018 WET photometry (upper left and lower right panels) as well as the FTs of the 2018 SpeX synthetic phototometry with wavelength shown in microns, as indicated in each panel.For the synthetic photometry at 1.51 μm, data from the first night of this 2018 run were not included because of excess noise at this wavelength, presumably due to excess sky noise in the adjacent sky band.The pulsation frequencies at 1113 and 1211 μHz are clearly visible in the optical and nearinfrared bandpasses though weaken significantly at 1.6 μm.
The pulsation frequencies at 788 and 2225 μHz display the opposite behavior and appear strongest at the longest wavelengths.Figure 2 presents similar FTs for the 2020 data.At this epoch, all pulsation frequencies were weaker than in 2018, though the peak near 2000 μHz is clearly visible.The amplitude of this frequency decreases as a function of wavelength out to 1.6 μm, then it appears to increase.The lower pulsation amplitude of G29-38 is within its previously observed range (Kleinman et al. 1998).
The 2018 and 2020 optical light curves form the basis for the identification of the pulsation frequencies during the SpeX observations.A statistically significant frequency is defined as a frequency with an amplitude at least four times the average noise level, representing a 99.9% probability that the frequency represents a true signal in the data and is not the result of random noise.The noise level is defined as the the average amplitude remaining in the FT after prewhitening (discussed further below) by the dominant peaks (Provencal et al. 2012), and is a conservative estimate, especially for the 2018 data.In  2018, G29-38 was in a fairly high-amplitude pulsation state and exhibited amplitude modulation, and thus it was impossible to remove all of the significant power completely before calculating the noise level.The mean noise over a broad range of pulsation frequencies in 11 different frequency ranges (WET optical plus the 10 SpeX bands) for 2018 and 2020 are presented in Figure 3.During the 2018 campaign, the WET observations had a typical mean noise level of 1 mma.The IRTF observations, covering fewer pulsation cycles, had somewhat higher noise levels, rising up to ≈4 mma in a few cases.In particular, for the bluest SpeX band, at 0.78 μm, there was additional noise near 1000 μHz, and for the two reddest SpeX bands, at 2.23 and 2.41 μm, there was additional noise near 2300 μHz.The likely cause of this additional noise at these isolated frequencies is our inability to remove the signal produced by amplitude modulation completely.We also calculated Monte Carlo simulations for the WET data using the Monte Carlo routine in PERIOD04.This routine generates a series of simulated light curves using the original times, the fitted frequencies and amplitudes, and also containing Gaussian noise.Each simulated light curve is subjected to a least-squares fit.The uncertainties are produced by the distribution of fit parameters.The Monte Carlo results are consistent with a noise level of 0.5 mma for the 2018 optical data and 0.2 mma for the 2020 optical data.
Armed with an understanding of the noise properties of these data, frequency identification is now possible.In both data sets, the noise in the FTs derived from the WET runs is significantly lower than the FT noise from the SpeX synthetic photometry, which covered fewer stellar oscillations.Given the additional higher frequency resolution of the optical data, we use the WET optical frequencies as the basis to identify the frequencies in the SpeX FTs.
PERIOD04 and PYRIOD were used to prewhiten the FTs.The process involves identifying the largest amplitude resolved peak in the FT, or the Lomb-Scargle periodogram in the case of PYRIOD, fitting the data set with a sinusoid of that frequency, subtracting the fit from the light curve, recomputing the FT, examining the residuals, and repeating the process until no significant power remains.This process works well for stable pulsators, but is more complicated in the presence of amplitude modulation.We illustrate the process in Figure 4 with an example of prewhitening of the dominant power in the 2018 WET data.In the top panel of Figure 4, the red arrow identifies the dominant 1211 μHz peak.The middle panel shows the result of removing a sinusoid with that frequency, and identifies the next highest-amplitude frequency at 1112 μHz.The bottom panel shows the result of removing both the 1211 μHz and 1112 μHz peaks from the data.Significant power remains, but a resolved frequency cannot be clearly identified.This is the signature of amplitude modulation, and the remaining power cannot be prewhitened.
Table 4 contains the resulting frequency identifications, within the constraints imposed by amplitude modulation, for the WET 2018 and 2020 data sets.The table presents frequency and its uncertainty in units of μHz, amplitude and its uncertainty in units of mma, and the S/N of the measured amplitude, which is the ratio of the amplitude of the frequency to the amplitude of the noise in that frequency range.We use the symbols f and g to identify frequencies in the 2018 and 2020 data sets, respectively.The next step was to fit the pulsation amplitudes and phases for the 10 bands of synthetic SpeX photometry in order to determine how these varied as a function of wavelength.This was performed independently for the 2018 and 2020 data sets because G29-38 is, as stated above, not a stable pulsator on a yearly timescale.We fixed the frequencies to those derived from the WET time-series photometry because of the much greater time base of these data, and then used PERIOD04 to fit the amplitude and phases in the SpeX data.

Patterns Among Pulsation Properties
Figure 1 displays a surprising result.The frequencies at 1211, 1112, and 1999 μHz ( f 0 , f 1 , and f 2 in Table 4, respectively) are clearly visible in the optical and near-infrared bandpasses and weaken significantly at 1.6 μm.However, the frequencies at 788 and 2225 μHz display the opposite behavior.These two frequencies have been identified as combination modes, defined as exact numerical combinations (by addition or subtraction) of larger-amplitude, principal modes.In pulsating white dwarf stars, these combination modes are thought to arise from mode-mode interactions in the convection zone (Brickhill 1992b;Wu 2001).The two combination frequencies at 788 and 2225 μHz (corresponding to f 2f 0 and f 0 + f 3 in Table 4) increase in amplitude as a function of wavelength.
Figure 2 presents similar FTs for the 2020 data.At this epoch, G29-38 was in a lower-amplitude pulsation state.These lower pulsation amplitudes provided S/N 10 for the amplitude measurements of the g 0 (2001 μHz) mode in the Spex data, with weaker modes even harder to measure.For this reason, we drop further analysis of the 2020 data and focus on the amplitude and phase behavior of the pulsations in the 2018 data.
Figure 5 presents the amplitudes (top panels) and times of amplitude maxima (bottom panels) as a function of wavelength for pulsation frequencies with sufficient S/N to discern patterns among their properties.The left-hand panels display properties for the three highest-amplitude g-modes (principal modes f 0 , f 1 , and f 2 ) and the right-hand panels display the four other harmonic or nonlinear combination frequencies (2f 0 , f 0 + f 1 , f 0 + f 3 , and f 2 − f 0 ).For the principal modes, the amplitudes decrease with increasing wavelength, as expected for photospheric g-modes (see Brassard et al. 1995, Figure 3) with no significant luminosity contribution via dust absorption and reradiation.Additionally, the time of maxima for f 0 , f 1 , and f 2 are approximately constant versus wavelength except for the optical WET measurement for f 0 which lags the time of maxima of f 0 (period of 825 s) at most other wavelengths by about 1 minute.This one phase offset could be caused by a change in the orientation of the f 0 mode captured during the 2018 WET run, which had a longer temporal coverage than the 2018 IRTF observations.Except for this one WET time lag, the three principal modes again behave as one would expect for atmospheric pulsations with little to no contribution by reradiated light from the dust.Note.The identified frequencies are labeled sequentially in amplitude from f 0 for 2018 and from g 0 for 2020.The amplitudes of the combination frequencies ( f 0 + f 1 , f 0 + f 3 , and f 2 − f 0 ), follow a different pattern.The observed amplitudes of these combination modes generally decrease from the visible through near-infrared, then increase, particularly beyond 2 μm.While the increase in mode amplitude at the longest wavelengths is what one would expect for modes strongly reradiated by dust, the spectral energy distribution modeling of Reach et al. (2009) estimated that 34% of the flux at 2.2 μm is due to dust, with the remaining flux emitted by the star's photosphere.Assuming these pulsation amplitudes have the same fraction of dust contribution as the time-averaged star, even after correcting for the dust contribution the photospheric pulsation amplitudes increase through 2.5 μm.Although tentative, we note that this is in contrast to the prediction of Wu (2001), who finds that the ratios of the amplitudes of parent modes to the amplitudes of their resultant combination modes are independent of wavelength.The 2f 0 harmonic (discussed further below) may also increase in the infrared or may be consistent with a near-constant amplitude at all observed wavelengths.The time of maxima for these frequencies varies considerably, from essentially constant ( f 0 + f 3 ) to decreasing as a function of wavelength (2f 0 and f 0 + f 1 ).One combination mode in the fourth panel, f 2 − f 0 , has had 1100 s subtracted from its times of maxima to aid visualization.

Discussion
We now examine whether the optical and near-infrared mode characteristics (amplitude and phase) observed in 2018 are broadly consistent with a few simple models.We do not make this comparison to the 2020 data because the star's pulsations were too weak at the observed epoch.We seek to find a simple model that can explain how the three highest-amplitude pulsation frequencies ( f 0 , f 1 , and f 2 ), which are all principal pulsation modes, decrease in amplitude as a function of increasing wavelength, while the combination modes display an increasing amplitude with wavelength, particularly beyond 2 μm (Figure 5).Making the interpretation more difficult is that the mode identifications (the specific spherical harmonics) for any of the observed modes are not well known.The general dearth of consistent multiplets in G29-38, along with its empirically variable frequencies and amplitudes with time, have made mode identification challenging (Kleinman et al. 1998).Recent work using long-baseline Transiting Exoplanet Survey Satellite data has revealed some possible mode identification (Uzundag et al. 2023), but since those observations are separated from our IRTF data by many months (and G29-38 has shown changes on shorter timescales), we have not adopted those identifications here.

Information Encoded in the Combination Modes
The combination modes (all combinations with f 0 ) show the strongest infrared signatures.In general, combination modes have different geometries (in terms of the angular dependence of their surface brightness variations) than their underlying parent modes.Currently, with insufficient knowledge of the mode ℓ identifications, it is not possible to model the stellar surface geometry uniquely and potentially constrain the dust geometry.We search instead for a simple picture that unites these infrared signals from the combination models.
Principal pulsation modes have surface variations that can be described via spherical harmonic functions, ( Y , where Θ and Φ are, respectively, the angles of latitude and longitude with respect to some chosen axes in the star, ℓ is the spherical harmonic order, or total number of surface nodal lines, and m identifies the number of nodes along a line of latitude.Note that for stellar g-modes ℓ = 1, 2, 3,K and |m| ℓ for each mode.For a given mode with order ℓ, the relative amplitudes of the components with values m depend on the axes chosen.If the axis of rotation is known, it is commonly used to define the (Θ, Φ) coordinate system.Alternatively, if circumstellar dust were primarily in a plane, that and a perpendicular axis could define the longitude and latitude for mathematical convenience.There is, as yet, no such known natural coordinate system for G29-38.
The angular dependence of a combination mode with individual dependencies Y ℓ m 1 1 and Y ℓ m 2 2 is the product of these two functions (Brickhill 1992b;Brassard et al. 1995;Wu 2001), which equals the sum over a set of Y ℓ m with ℓ running from |ℓ 1 − ℓ 2 | to ℓ 1 + ℓ 2 .While the resulting geometry can be complicated, if ℓ 1 = ℓ 2 , then there is always a component in the sum which has ℓ = 0. (See also Brassard et al. 1995, Appendix C, where they demonstrate that either the product of two ℓ = 1 modes, one m = + 1, the other m = −1, or the product of two ℓ = 2 and m = 0 modes have a spherically symmetric Y 0 0 component, but the product of an ℓ = 1 mode and an ℓ = 2 mode has no such component.)Therefore, if the combination mode results from two modes with the same ℓ value, there will be a component of the combination mode that has no angular dependence (see also Kurtz 2005).This component is an isotropic pulsation, much like a radial-mode pulsator such as an RR Lyrae or Cepheid.In such a case, whatever the dust distribution, such a component will always provide a signal in the reradiated flux.We propose the possibility that the observed 2018 combination modes were visible in the reradiated flux precisely because each contributing mode had the same value of ℓ.This suggestion offers a prediction and possibly a constraint for future mode identification.

Historic Data
Historic time-series optical and infrared data were presented by Graham et al. (1990) and Patterson et al. (1991).The data presented by Patterson et al. (1991) seems to have slightly higher quality, judging by the number of pulsation frequencies detected, but are in agreement with those of Graham et al. (1990).Patterson et al. (1991) detected pulsations in B with periods of 615 s (1620 μHz), 186, 243, and 268 s.They also detected pulsations in K with periods of 186, 242, and 272 s, which are consistent with the pulsations periods seen in B. They interpreted the 615 s period as a principal stellar mode, with no corresponding dust response.They envisaged the dust as being in the form of an optically thin circular disk (see their Section 4.4.1), and suggested that this mode (with respect to the axes aligned with the disk) had m ≠ 0, so that it elicited no response from the dust disk.They argued that the amplitudes of the other pulsation frequencies were too high in K to be just stellar g-modes (indeed Graham et al. 1990 found these pulsations only in K and not in B), and thus that they must represent a dust response.However, since the theory of combination modes, and in particular their surface structures, was only just being developed (Brickhill 1992b;and later Brassard et al. 1995;Goldreich & Wu 1999;Wu 2001), they assumed the stellar surface brightness variations to be of the form relevant to the principal modes, viz, Thus in order to explain why the pulsations seen in K were not seen in B they needed to appeal to a special geometry for the assumed stellar g-modes relative to the disk and to the observer.In particular, the infrared-strong modes all needed to have ℓ = 2 and m = 0 in axes relative to the disk, and the disk needed to be oriented at an angle of ≈55°relative to the observer.In our interpretation, while the 615 s period arises from a principal stellar mode (most likely with ℓ = 1; Kleinman et al. 1998), the other detected pulsation frequencies that are all strong in the infrared would be combination modes, whose fundamental contributors all have the same value of ℓ.In this interpretation a special geometric configuration is not required.

The First Harmonic
The oscillation seen at frequency 2f 0 is the first harmonic of the parent mode f 0 , which gives rise to the combination modes that we find to have strong visibility in the near-infrared.
G29-38 is well known as a high-amplitude oscillator, and as such it would be no surprise if the dominant mode that we find, the f 0 mode, were of sufficiently high amplitude that its underlying variability is slightly nonlinear.Such nonlinearity would primarily manifest itself as a first harmonic in the light curve with period 2f 0 .If the f 0 mode has surface angular dependence Y ℓ m , then this contribution to the variability at frequency 2f 0 would also have the same Y ℓ m surface distribution.As an additional effect, for a single principal mode (Wu 2001), the sinusoidal oscillation at the base of the convection zone is modified as it passes through the zone.Thus the emergent oscillation is no longer purely sinusoidal, but can be represented as a Fourier series with harmonic frequencies 2f 0 , 3f 0 ,....These are seen because the original oscillation has been distorted.The magnitude of the coefficient for the first harmonic 2f 0 contribution caused by this effect is proportional to the square of the amplitude of the original mode, and is therefore nonnegative.This implies that, if the original mode f 0 has a surface distribution of the form Y ℓ m , the first harmonic contribution 2f 0 caused by this effect has surface distribution ´. Thus this contribution to the first harmonic always has a surface component that is isotropic.For this reason the first harmonic frequency 2f 0 can, in principle, be more prominent than the underlying parent mode.As Wu (2001) comments: "This arises because the apparent amplitudes of the higher ℓ modes suffer stronger cancellation when integrated over the stellar disk, while the harmonics of these modes do not."In this respect the first harmonic induced by passage through the convection zone is analogous to the behavior of combination modes with the same value of ℓ (Section 4.1).We note that in our data the near-infrared response elicited by the 2f 0 harmonic is significantly less than for the combination modes that we detect (Section 3.2).It is possible that this comes about because of the two possible sources of the oscillation at frequency 2f 0 , which, as we have noted, have two different surface brightness distributions.
The wavelength behavior of 2f 0 may aid future modeling of g-mode stellar pulsations in constraining the properties of the convection zone (Montgomery et al. 2020), and in that spirit we present additional analysis of this harmonic.For example, in contrast to expectations (Wu 2001), we find that the ratio of the amplitudes of the harmonic (2f 0 ) to the fundamental ( f 0 ) depends strongly on wavelength.Figure 6 presents the f 0 mode (blue dotted sinusoid) over two periods along with superpositions (additive combinations) of f 0 and 2f 0 at different wavelengths.For the pulsations measured at 0.5 μm, the frequencies, amplitudes, and phases are all fit from the 2018 WET time-series photometry.For pulsation frequencies measured at longer wavelengths, the frequencies are fixed at the values determined from the 2018 WET data, but the amplitudes and phases are fit from the 2018 SpeX synthetic photometry.At 0.5 μm, the superposition of 2f 0 and f 0 creates a curve that departs slightly from a sinusoid with a slower rise and steeper fall.This contrasts with the prediction of Wu (2001; see also Brickhill 1992a, Figure 7) that such distortion should lead to peaked light curves with sharp ascents and shallow descents.At 1.51 μm the superposition departs more clearly from a sinusoid and by 2.41 μm the amplitudes of f 0 and 2f 0 are similar (see also Figure 5).Measuring the mode properties of the principal and harmonic frequencies across such a range of wavelengths may offer new constraints on the properties of surface convection zones or their time-dependent depth changes for these stars.

Simple Dust Morphology Models
We note that our observational data extend only to 2.5 μm, where emission from the hottest dust, nearest to sublimation (at the inner edge of a disk or the dust orbital periastron) should peak.Thus any morphological information that can be derived from the pulsations constrains only the inner edge of the circumstellar environment.Additionally, at these wavelengths, the dust is only marginally the dominant contributor to the flux (e.g., Xu et al. 2018); careful modeling will be required to disentangle the atmospheric and dust contributions properly.With these caveats in mind, we proceed to an initial exploration of how the infrared amplitudes of the combination modes constrain the morphology of the circumstellar dust.We start by assuming that the dust is more or less in a plane, of some indeterminate thickness.The dust might be a circular disk, an eccentric disk, or even in a parabolic stream (e.g., Nixon et al. 2020).Given the theories for the origins for the dust, there is no reason to suppose that the rotation axes of the dust and of the white dwarf are in any way aligned.In considering the illumination of the dust by the white dwarf, the natural coordinates have the Θ = 0 axis (direction toward the north pole) perpendicular to the dust plane.In this case, for example, a (2, 0) mode on the star (relative to the starʼs rotation axis, which is unlikely to be aligned perpendicular to the dust plane) would look to the dust as a collection of the multiplets (2, −2), (2, −1), (2, 0), (2, +1), and (2, +2).Similarly, any ℓ = 1 mode on the star would look to the dust as a collection of the multiplets (1, −1), (1, 0), and (1, +1).

Optically Thin Dust
We consider first the simplest case in which the dust is optically thin at all wavelengths.In what follows, we continue to work in a coordinate system based on the dust plane.We consider some simple dust configurations (see Cotton et al. 2020): 1.If the dust is distributed uniformly and spherically symmetrically then the only dust response will be to the (0, 0) component.This surface brightness distribution of the modulation only occurs as a combination mode.2. If the dust is distributed as a uniform circular disk of finite vertical thickness, then by symmetry none of the m ≠ 0 components produce a net dust response.In addition the (1, 0) component produces no net response.Yet the (2, 0) component can produce a nonzero net response (see Figure 9 in Graham et al. 1990).In this case, if the underlying basic modes are all ℓ = 1, we expect zero dust response from these, and the only nonzero dust response to be from the combination modes.But if the underlying modes are all ℓ = 2, then we expect a dust response from both the underlying modes and from the combination modes.3.If instead the dust is in an eccentric disk (of which an extreme example would be a parabolic stream), then as seen from the star, the dust density will have an m = 1 ( i.e., cos Φ) component.In this case there will also be a response from the underlying m = +1/−1 components but not from the underlying m = +2/−2 components.Note that in this case there can be a phase difference between the stellar pulsations and dust response depending on the geometry.
Thus for these simple dust configurations, one should always expect a response from combination modes (of equal ℓ), as they always have a (0, 0) component.We might expect a response from a (2, 0) component.Note that this could be a (2, 0) component of a principal pulsation mode or part of the combination response from (1, 0) × (1, 0) or (2, 0) × (2, 0) (see the appendix of Brassard et al. 1995).If the dust distribution is circular, then we expect no response to ℓ = 1 components.But if it is eccentric then there can be a response to such components, with the amplitude dependent on the degree of eccentricity.
Overall, if we see only a response to combination modes and none to the principal modes, the best guess would be that the underlying modes have ℓ = 1 and the dust distribution is circular.But, importantly, there are uncertainties here depending on the expected amplitudes of the response and the observational uncertainties.For example, the f 0 , f 1 , and f 2 modes do show amplitudes at 2.5 μm; one would need to be able to disentangle how much of these oscillations at this wavelength might be due to the dust response.Such detailed modeling is beyond the scope of the present paper.

Optically Thick Dust
If the dust is optically thick, and indeed that seems probable given the large infrared excess of G29-38, then there are too many possibilities to be considered here.We restrict the discussion to two obvious simple generalizations, based on previous theoretical ideas and spectral modeling: 1. Consider the dust to be arranged in a circular, infinitesimally thin disk, which is optically thick in the vertical direction.In this case none of the m ≠ 0 components contribute to the infrared response.But both the (1, 0) and the (2, 0) components would give rise to a response.2. Consider the dust to be in a circular disk of finite thickness so that the main dust response comes from the inner edge of the disk (see Ballering et al. 2022).In this case, assuming an intermediate inclination angle, one expects no contribution from (1, 0), but contributions from (2, 0) as well as from the m ≠ 0 components, but not from m = +2/−2 components.
While our observations of higher infrared contributions in the combination modes may favor some of these models, in order to rule out any of them definitively requires modeling the variations of amplitudes and phases as functions of wavelength.

Conclusions
We obtained near-infrared spectrophotometry of G29-38 using SpeX at the IRTF in 2018 and 2020, along with contemporaneous optical WET time-series photometry.We detected six principal pulsation modes, one harmonic, and three combination modes in the 2018 data set.Pulsations during the 2020 observations were too weak to see clear trends in amplitude as a function of wavelength.Among the 2018 pulsations with sufficient S/N to measure across the full range of observed wavelengths, the three principal modes all showed declining pulsation amplitudes with increasing wavelength.These pulsation amplitudes decreased from 20 to 35 mma in V, g, and the other WET optical bands to ≈5 mma at 2.5 μm.The three nonlinear combination modes from 2018 exhibited strikingly different behavior, rising from approximately 6-9 mma in the optical to 8-19 mma in the near-infrared.The harmonic of the principal f 0 mode may show similar behavior, though at a reduced level.
The basic theory of the origin and the behavior of harmonics and combination modes for stars with a surface convection zone, which includes all ZZ Ceti stars, has been set out by Brickhill (1992aBrickhill ( , 1992b) ) and by Wu (2001).We have found some discrepancies with the expectations laid out in those papers.These papers predict for a single mode that the ratio of the amplitude of the first harmonic to the amplitude of the underlying mode should be independent of wavelength and that the effect of the first harmonic is to give rise to a peaked light curve with a sharp ascent and shallow descent.It is evident from Figures 5 and 6 that our results disagree.Also, these papers predict that the ratios of the amplitudes of the principal modes to the amplitudes of their resultant combination modes are independent of wavelength.While we have not performed modeling to separate the dust versus photospheric contributions, from Figure 6 and the time-averaged spectral energy distribution modeling of Reach et al. (2009), it appears that our results disagree.These findings warrant further investigation.
We conclude that the dominant near-infrared response is most likely due to the isotropic component of the combination modes.Unfortunately, isotropic pulsation components are of no use in determining the geometric structure of the dust.To make future gains in determining the dust morphology, we need to measure the dust response to modes with ℓ 1, i.e., principal modes and perhaps nonisotropic components of the combination modes.This will require careful modeling.Additionally, in order to clarify what fraction of the infrared pulsations is from the dust response versus from the stellar photosphere, we need detailed modeling of atmospheric pulsation modes, which depend on identifying the spherical harmonics of these modes.An alternative approach would be to observe at wavelengths longer than 5 μm, where G29-38ʼs spectrum is almost entirely due to the heated dust (Reach et al. 2009).Because of the timevarying nature of the stellar pulsation modes, some epochs have only weak pulsations that are insufficient to constrain the dust environment.However, whenever G29-38 exhibits different pulsation frequencies these may provide new diagnostics and the dust distribution itself may vary with time.The phase behavior of the pulsation frequencies also warrants further investigation and may help resolve the dust and surface pulsation geometries.
While not a goal of this study, our analysis of the combination modes likely indicates that the observed principal modes all have the same ℓ value.This constraint is both a prediction and a potentially useful tool for future mode identification on this important star.
.4 Note.The ERAU 1.0 m is at Embry Riddle Aeronautical University in the USA.The Krakow 0.5 m is at the Astronomical Observatory of the Jagiellonian University in Poland.The Molėtai 0.35 and 1.65 m are at the Molėtai Astronomical Observatory in Lithuania.The SARA-CT 0.6 m is at the Cerro Tololo Inter-American Observatory in Chile.The SARA-RM 1.0 m is at the Roque de los Muchachos Observatory in Spain.The TSAO 1.0 m is at the Tien Shan Astronomical Observatory in Kazakhstan.differential extinction.The result was a series of light curves with times in seconds and amplitude variations represented as fractional intensity.The final step in the reduction process combined the individual light curves from the Whole Earth Telescope (WET) photometry and applied barycentric corrections.For this step, it was assumed that G29-38 oscillates around a mean light level.This important assumption allowed the correction of instrumental intensity offsets for any overlapping light curves.The question of the treatment of overlapping data is discussed in detail inProvencal et al.

Figure 1 .
Figure 1.FTs for the 2018 WET photometry, labeled "optical," and the 10 2018 SpeX bands (labeled with their central wavelengths).The WET optical FT is repeated in the lower right panel to aid in visualizing the amplitudes of the pulsation frequencies vs. wavelength.The pulsation frequencies at 1113 and 1211 μHz are clearly visible in the optical and near-infrared bandpasses, though weaken at wavelengths 1.6 μm.The pulsation frequencies at 788 and 2225 μHz display the opposite behavior and appear strongest at the longest wavelengths.

Figure 2 .
Figure 2. Same as Figure 1, except based on the 2020 November data.The strongest pulsation frequencies are approximately half the amplitudes of those in 2018 September.

Figure 3 .
Figure 3. Mean noise as a function of frequency for the optical photometry and the SpeX synthetic photometry.The increased noise in several of the SpeX channels, particularly in 2018 at 0.78 μm near 1000 μHz and at 2.23 and 2.41 μm near 2300 μHz, are likely due to our inability to remove the dominant frequencies completely.

Figure 4 .
Figure 4. Demonstration of prewhitening using the dominant 1211 μHz frequency in the 2018 WET data.The first panel is the original FT.The red arrow identifies the dominant 1211 μHz peak.The second panel has been prewhitened by 1211 μHz, and the third panel is prewhitened by both 1211 μHz and 1112 μHz.The green line indicates the formal noise level.

Figure 5 .
Figure5.The top two panels display the pulsation amplitudes measured during 2018 as a function of wavelength.The three highest-amplitude g-modes (lefthand panels) clearly decrease in amplitude with wavelength.The amplitude of the first harmonic (2f 0 in the right-hand panels) of the highest-amplitude gmode ( f 0 ) displays a different behavior.The amplitudes of the combination modes ( f 0 + f 1 , f 0 + f 3 , and f 2 − f 0 in the right-hand panels) increase beyond 2 μm.The bottom panels display the times of maxima for these pulsations as a function of wavelength.

Figure 6 .
Figure6.Comparison of the f 0 mode (blue dotted sinusoid) with the superposition of f 0 and 2f 0 at three different wavelengths (black curves, see legend).The gray shaded zone covers the ±1σ uncertainties in pulsation amplitudes and phases.

Table 1
Bandpasses Used in the SpeX Synthetic Photometry

Table 3
Journal of 2020 Time-series Photometry