Simulation of Stark-broadened Hydrogen Balmer-line Shapes for DA White Dwarf Synthetic Spectra

White dwarfs (WDs) are useful across a wide range of astrophysical contexts. The appropriate interpretation of their spectra relies on the accuracy of WD atmosphere models. One essential ingredient of atmosphere models is the theory used for the broadening of spectral lines. To date, the models have relied on Vidal et al., known as the unified theory of line broadening (VCS). There have since been advancements in the theory; however, the calculations used in model atmosphere codes have only received minor updates. Meanwhile, advances in instrumentation and data have uncovered indications of inaccuracies: spectroscopic temperatures are roughly 10% higher and spectroscopic masses are roughly 0.1 M ⊙ higher than their photometric counterparts. The evidence suggests that VCS-based treatments of line profiles may be at least partly responsible. Gomez et al. developed a simulation-based line-profile code Xenomorph using an improved theoretical treatment that can be used to inform questions around the discrepancy. However, the code required revisions to sufficiently decrease noise for use in model spectra and to make it computationally tractable and physically realistic. In particular, we investigate three additional physical effects that are not captured in the VCS calculations: ion dynamics, higher-order multipole expansion, and an expanded basis set. We also implement a simulation-based approach to occupation probability. The present study limits the scope to the first three hydrogen Balmer transitions (Hα, Hβ, and Hγ). We find that screening effects and occupation probability have the largest effects on the line shapes and will likely have important consequences in stellar synthetic spectra.


Introduction
White Dwarfs (WDs) represent the end stage of the vast majority of stars (>97%) in our universe (Van Horn 2015). Their intrinsic properties and ubiquity make them useful in a surprisingly wide range of astrophysical contexts. One of the most well known and historically well established of these is cosmochronology (Winget et al. 1987;Fontaine et al. 2001). Perhaps the most widely relevant application for the general astronomical community as a whole is the fact that hot WDs are used to calibrate many major observatories and instruments including Hubble, the Sloan Digital Sky Survey, ESO, and Kitt Peak (Narayan et al. 2019). In all of these contexts, the fidelity of the information we can extract is predicated on the accuracy of our model atmosphere calculations.
The inferences we can make using WDs rely on determinations of the fundamental stellar properties: effective temperature (T eff ) and surface gravity (log g). A grid of model atmospheres can be used to derive fundamental parameters of DA WDs via spectroscopic and photometric methods (Bergeron et al. 1992;Koester 2008;Hubeny et al. 2021). In the spectroscopic method, the individual absorption lines are fit to the synthetic spectra whereas, in the photometric method, fits are done using broadband integrated flux values. The best-fit solution gives the fundamental parameters of the star.
Since the first application of spectroscopic fitting to a large sample of DA WD stars by Bergeron et al. (1992), it has become the most widely used technique to determine the effective temperatures and surface gravities of these stars. The popularity of this method can be attributed to its precision. Characteristics of hydrogen line profiles are highly sensitive to changes in atmospheric conditions (Schulz & Wegner 1981). Therefore, relative changes in log g and T eff manifest as variations in the widths and strengths of the lines and are quantifiable to a high degree of precision.
The success of the spectroscopic technique is well established by the myriad high-impact results borne from its application to large samples of DA stars. Determinations of the shape of the DA population mass distribution, the boundaries of the ZZ Ceti instability strip, characterizations of the luminosity function, and the evolution and formation rate of DA WDs were all performed using the spectroscopic method (Bergeron et al. 1992;Liebert et al. 2005;Eisenstein et al. 2006;Gianninas et al. 2006;Kepler et al. 2007). Despite the litany of successful applications of this technique, its accuracy is fundamentally constrained by the fidelity of the physics used to construct the model atmospheres. It has become clear in recent years, as data quality improves, that the existing physical models are no longer sufficient to faithfully reproduce the observations. A particularly egregious example of this is the highly discrepant set of masses and temperatures obtained from fitting different lines in the same star (Fuchs 2017). Perhaps even more concerning is the emergence of a clear systematic offset between temperature and gravity determinations in largescale comparisons between the spectroscopic versus photometric methods. Spectroscopically inferred temperatures are roughly 10% higher and spectroscopic masses are roughly 0.1 M e higher than their photometric counterparts Genest-Beaulieu & Bergeron 2019). Given this, a reevaluation of the fundamental physical assumptions and models used in model atmospheres is well justified.
One likely source for the discrepancies described above is the way that the line shapes used in model atmospheres are calculated. The key to determining the stellar fundamental quantities is our ability to calculate the way that light interacts with the matter in the atmosphere before emerging from the surface of the star. WD atmospheres are highly ionized and dense enough that plasma effects have a significant effect on the bound-bound transitions. In this type of environment, pressure broadening dominates the line shapes and, by extension, impacts the shape of the absorption features in the emergent stellar spectrum. The line shapes are used to calculate cross sections of bound-bound transitions, which are then used to compute opacities and radiative rates. These values come into play when solving the radiative transfer equation governing the propagation of light through the layers of the atmosphere. Ultimately, this yields the final emergent stellar spectrum. Therefore, the set of line shapes provided to the code represents one of the most fundamental pieces of physics used to construct WD model atmospheres. They are one of the key pieces of constitutive physics that fundamentally dictate our ability to infer masses and temperatures and use them effectively in our astrophysical models.

Approaches to Line-shape Calculations and Recent Issues
In its most fundamental form, a line shape is a curve that provides a probability distribution for a bound-bound transition of a given photon energy. Model atmosphere codes use a grid of broadened hydrogen line profiles that spans the parameter space of possible transitions, temperatures, and electron densities found in WD atmospheres to synthesize emergent spectra. These codes have historically relied on line profiles calculated using the model of Vidal et al. (1970), known as the unified theory of line broadening (VCS). VCS unified the impact approximation and one-electron theory and calculated profiles from the cores out to the far wings for the first time in a self-consistent way. Since then, the VCS line shapes have been used as the standard set of hydrogen lineshape calculations in model atmosphere codes.
However, there are several complications that precluded model atmosphere codes ability from properly reproducing observed spectra using the VCS line shapes. In particular, Bergeron (1993) observed systematic discrepancies in the stellar effective temperature inferred using different individual lines. Fitting to higher Balmer lines returned cooler temperatures than lower ones. The discrepancies were attributed to the observation that the model continuum opacity was too large in the regions between the lines where the wings of adjacent lines overlap. Early attempts to fix the problem were done by incorporating "nonideal effects" in the line-shape calculations.
Physically, these nonideal effects have to do with the treatment of level dissolution, which is also often referred to in the literature as "continuum lowering" or "ionization potential depression" (Inglis & Teller 1939). In short, whereas isolated atoms have fixed energy levels that the electrons can occupy, atoms embedded in high-density plasmas experience large energy-level perturbations in the form of Coulomb interactions from the surrounding charged particles. When the perturbations become large enough, higher atomic levels are assumed to be "dissolved." Because an electron can no longer occupy the level, transitions to and from this level no longer exist. In practice, this level dissolution can be enforced by limiting the electric field magnitudes that atoms in the plasma can experience. Beyond some upper bound, the transition is not allowed to occur because the electron will not be able to occupy the level. The term "nonideal effects" is a reference to the fact that the formalism used to set the dissolution point originally comes from a nonideal equation of state that was intended to parameterize the partition function for a nonideal gas.
Incorporating this effect changes not only the intensity of the line but also its shape, with the largest change being a dramatic drop in the intensity of the wings. The diminished line opacity is then accounted for by extra continuum opacity ("pseudocontinuum") that can extend to lower energy than the continuum transition energy. One well-known formalism used to parameterize this effect is the occupation probability formalism of Hummer & Mihalas (1988, hereafter HM88).
In using the HM88 formalism to fit the experimental emission data of Wiese et al. (1972), Bergeron (1993) noticed that their model predicted too much emission between the lines, presumably because their calculated lines were too broad. In order to mimic the effects of occupation probability on line shapes while using VCS line profiles, they suppressed the flux in the region between the lines by minimizing the contribution of the pseudo-continuum. They achieved this through an ad hoc doubling of the maximum microfield required for state dissolution ("2 × β crit "). Originally intended as a temporary measure, this implementation became standard until the work of Tremblay & Bergeron (2009, hereafter TB09). TB09 combined the VCS approach with the prescription of HM88 in the line-profile calculation itself. This removed the arbitrary factor of 2 from β crit and made the new implementation selfconsistent.
The TB09 line shapes have since replaced VCS and are now well established as the standard set of line-shape calculations used for WD model atmosphere calculations. However, we note that the main difference between the two is a slight modification to the electric field distributions and an application of electron damping applied after the fact to the calculated profile. The fundamental assumptions of VCS are otherwise unchanged in TB09. A reevaluation of the line shapes used in model atmospheres must therefore start with an interrogation of VCS.

Theoretical Advancements and the Simulation Code Xenomorph
While VCS was pioneering work, it has limitations that cause it to fail at high density. Since VCS was first developed, there have been many developments that improve on the fundamental assumptions made by VCS and benchmark both theory and experiment. Many of these have been incorporated into various simulation codes (e.g., Stamm & Voslamber 1979;Gigosos & Cardenoso 1987;Stambulchik & Maron 2006;Gomez et al. 2016). While different codes incorporate different pieces of additional physics, there are three main changes that generate more accurate line shapes than VCS. The first is a higher-order multipole expansion of the Coulomb interaction between radiators and plasma perturbers compared with the dipole approximation of VCS. This allows for gradients in the electric field across the emitting atom. The second change allows for the expansion of the basis set of atomic states (within computational limitations), which leads to higher-order Stark terms. Finally, many simulation codes allow for dynamic ions as opposed to the static ion treatment in VCS. Perhaps surprisingly, none of the improved line-shape calculations mentioned above has been used to construct model atmospheres of WD stars. Xenomorph, the code initially developed by Gomez et al. (2016) and used for the work presented in this paper, implements all three theoretical advances.
Xenomorph is a semiclassical line-shape calculation code that performs the line-shape calculation in two main steps. First, it constructs a spherical ensemble of classical point particles that represent the plasma perturbers and evolves them along straight-line trajectories. The center of the simulation sphere is taken to be the position of the radiator, which in our case is a neutral hydrogen atom. The end product of the first portion of the code is the time history of the total Debyescreened electric potential at the position of the neutral hydrogen atom due to the ensemble of plasma perturbers. In the second portion of the calculation, the time sequence of electric potentials is used to integrate the Schrödinger equation and produce the final line profile. We demonstrate that when properly configured to match the approximations made by VCS, Xenomorph can reproduce the VCS line shapes. The additional pieces of physics therefore represent a credible improvement over VCS theory.
The statistical noise inherent in these calculations is affected by the choice of numerical method used to calculate the line profile from the time-evolved dipole moment operator. Another important consideration is the amount of randomness available to the simulation. We implement a novel approach to particle reinjection, which avoids strategies that artificially impose adherence to the Maxwellian velocity and impact parameter distributions.
This study will mark the first instance in which simulationbased hydrogen Balmer-line-shape calculations are developed and intended to be used in WD model atmospheres. We examine independently the effect of each of the new pieces of physics on the line shapes themselves and compare them to those of previous theories that precede the theoretical advancements. To the extent possible, we demonstrate that the inclusion of the new physics reproduces benchmark laboratory results better than previous theories. With the exception of the effect of ion dynamics on Hα, the theoretical advancements in the treatment of Stark broadening (the new pieces of physics included in Xenomorph) have not been propagated through a model atmosphere code and the effects on fundamental parameters have not been systematically investigated. We will explore this in greater detail in future work.
Given the significance of the effect of occupation probability on g log and T eff determinations, we incorporated two different prescriptions of occupation probability in Xenomorph based on the Hummer & Mihalas (1988) and Fisher & Maron (2002) approaches. This study marks the first instance, to our knowledge, in which occupation probability has been incorporated into simulation-based line-profile calculations. We also demonstrate that under the same occupation probability prescription used by TB09, we can reproduce their profiles, which are calculated using the VCS theory modified to incorporate the HM88 prescription of occupation probability.
The flexibility afforded by Xenomorph presents a valuable opportunity to explore the theoretical advancements made in line-shape theory since VCS. The code has been vetted against its predecessors and against laboratory data. It now allows us to parameterize the extent to which the new physics affects emergent stellar spectra as well as make direct comparisons against current standards used in DA WD model atmospheres and their spectroscopic and photometric fits. In this preliminary study, we limit the scope to low principal quantum number lines: Hα, Hβ, and Hγ (n = 3 → 2, 4 → 2, and 5 → 2, respectively). We will investigate higher-n lines in the future. We will also make the full grids of line profiles freely available to the WD community. Incorporating these effects into model atmosphere codes can potentially help make progress in the effort toward resolving the newly discovered systematic discrepancies in mass and temperature.
In short, the goal of this paper is to examine the myriad factors which affect hydrogen Balmer-line shapes at the conditions of interest in DA WD atmosphere calculations. In Section 2, we provide a brief overview of the theory of hydrogen line broadening in the context of WD atmospheres, discuss specifics of the simulation line-shape code Xenomorph, and describe its departures from other simulation codes and our choice of numerical method. In Section 3, we detail the new pieces of physics incorporated in the Xenomorph line-shape calculations that distinguish it from the VCS theory. We describe a preliminary set of fits to laboratory data in Section 4 and present our conclusions in Section 5.

Line Broadening in WD Atmospheres
DA WD atmospheres are made up of a dense plasma of essentially pure hydrogen. The radial extent of the atmosphere spans a range of number densities and temperatures along a continuous gradient. To construct synthetic spectra, the atmosphere is modeled as a set of thin, discretely stratified layers. The set of layers in aggregate will span the same number densities and temperatures that the actual continuously varying atmosphere does (see Figure 1).
Investigating the various depths of formation along a given absorption feature in the emergent stellar spectrum allows us to parameterize the plasma conditions at the layer of the atmosphere in which that portion of the line was formed. Figures 2 and 3 show wavelength-dependent electron density and temperature at the layer in the model atmosphere where the wavelength-dependent optical depth τ ν = 2/3, the depth point that most closely characterizes the emergent flux (Hubeny & Mihalas 2014). The depth of line formation can vary significantly for different stellar effective temperatures and surface gravities. In particular, the depth of formation of the core of the line is more sensitive to changes in g log , and the depth of formation in the wings of the line is more sensitive to changes in effective temperature.
Each set of conditions for each layer of the atmosphere requires a corresponding hydrogen line shape calculated (or interpolated) at the same conditions to compute the appropriate opacities and perform the radiative transfer. Because the fundamental parameters inferred by the spectroscopic method rely sensitively on the shape of the features in the emergent stellar spectrum, they also rely sensitively on the input Starkbroadened line shapes. Any inaccuracies in the line-shape calculations will propagate through to the mass and temperature determinations. Figures 1, 2, and 3 also show, in gray shaded regions, the typical plasma conditions achieved by the White Dwarf Photosphere Experiment on the Z-machine at Sandia National Laboratories (Rochau et al. 2014;Montgomery et al. 2015). While valuable theoretical improvements have been achieved, laboratory data available to validate the codes are sparse. Recently, measurements of hydrogen Balmer lines have been collected in both emission and absorption in a region of parameter space relevant to DA WD atmospheres Schaeuble et al. 2019). Efforts to validate Xenomorph against these experimental data are currently underway.
The primary pressure-broadening mechanism at conditions where there is significant ionization is collisions from electrons and ions. These collisions can be approximated by a timedependent Stark effect. Charged particles induce perturbations in the energy levels of the hydrogen atoms that correspond to the basis set of stationary state wave functions. The ensemble effect of the close-range electric field interactions between the charged plasma perturbers and neutral radiators produces the characteristic Stark-broadened line shapes.
The Hamiltonian for a neutral hydrogen atom (radiator or absorber) can be written as where we have placed the radiator at r = 0. V ext (r, t) is the timedependent Coulomb interaction between the atom and plasma; this is commonly treated as the net electrostatic potential produced by the plasma perturbers at the position of the radiator. In our case, for a pure hydrogen plasma, the perturbers are simply electrons and protons. H 0 is the unperturbed Hamiltonian of the isolated hydrogen atom in cgs units: where m e is the mass of the electron. The line shapes for spontaneous emission and absorption are assumed to be identical aside from the fact that they are inverse The blue shaded region corresponds to the parameter space covered by the VCS line profiles calculated by Lemke and used to calculate opacities for bound-bound hydrogen transitions (Lemke 1997). The vertical dotted lines show conditions at three different optical depths, and the horizontal dashed lines point to the conditions at an optical depth of τ ross = 2/3. The gray shaded region demonstrates the typical range of experimental conditions achievable on the Z-machine. Here, we plot the atmospheric structure for two different atmosphere models: one with the VCS line shapes and another with the new Xenomorph profiles containing all of the new features. This ensures the maximum amount of difference between the profiles. The atmospheric structure remains unchanged. The largest difference observed between the two models is less than a percent in the temperature at an optical depth of ≈10. Therefore, the differences in the profiles are not significant enough to induce a change in the temperature and density profiles of the atmospheres.
processes. We note that investigations of differences between emission and absorption profiles have been scrutinized under the framework of partial frequency redistribution (Hubeny & Lites 1995;Hubeny & Mihalas 2014, ch. 10 and 15). The studies conclude that the high densities in WD atmospheres make the effects of partial frequency redistribution completely negligible. Though recent laboratory experiments may suggest that this assumption may not always hold, we maintain for now the existing convention of calculating emission profiles and taking their inverse to be the absorption profiles.

Simulation
Xenomorph (Gomez et al. 2016), the simulation code we use for the line-shape calculations, performs a semiclassical calculation of the line shape where the atom is treated quantum mechanically and the plasma is treated as an ensemble of classical point particles with definite position and momentum. The line shape is produced by spontaneous emission from a neutral H atom (radiator) using a time sequence of perturbing electric fields from the evolution of a classical plasma. While some properties, such as electron capture (Gomez et al. 2020), can only be modeled by a full quantum mechanical treatment of the plasma perturbers, Smith et al. (1969aSmith et al. ( , 1969b show that a classical treatment of plasma perturbers is valid for weakly coupled plasmas. Treating the perturbers as classical particles is valid when the average distance between particles far exceeds the thermal de Broglie wavelength (Baranger 1958a). Therefore, the classical picture is a good approximation for the relevant plasma conditions in WD atmospheres. The code makes use of the Heisenberg picture, which evolves the operators in time, in contrast with the Schrödinger picture, which evolves the wave functions (or state kets). The mechanics involved in the computation of the line shape (which only requires knowing the time history of the dipole moment operator) using either the autocorrelation or power spectrum method (discussed in Section 2.4) makes the Heisenberg picture the natural choice.
The dynamics of the system are constructed first. The plasma perturbers are initially placed randomly throughout the spatial extent of a spherical cavity and are assigned velocities taken from a Maxwellian distribution at a given temperature. The perturbers are evolved along straight-path trajectories. This is not exactly physically accurate because interparticle Coulomb interactions should alter the trajectories. The effects of interparticle interactions are mimicked using a Debye-screened Coulomb potential in lieu of solving the full N-body problem (Stambulchik et al. 2007). The nuances of the choice of screening prescription are discussed in more detail in Section 2.4.
The time history of the evolution of the plasma is then used to calculate the net electrostatic potential, V(r, t), produced at the center of the simulation sphere (the position of the radiator). This net electrostatic potential is added to the unperturbed Hamiltonian, H 0 , for the isolated hydrogen atom to construct the full time-dependent Hamiltonian, H(t). The time-dependent Hamiltonian is then used to calculate the time-evolution Figure 2. Electron density at the depth of formation. Top panel: electron density at an optical depth of τ ν = 2/3 for three DA WD atmosphere models with T eff = 15,000 K at three different g log . Bottom panel: same as the top panel but for T eff = 20,000 K. The higher opacity in the line centers means light from these wavelengths is emitted from shallower depths, higher in the atmosphere where the density is lower. Conversely, the wings are emitted from denser regions deeper in the atmosphere.  Figure 2 but for temperature. At these stellar effective temperatures, the Balmer region of the spectrum is formed higher and at cooler temperatures than the stellar effective temperature. operator, U(t), which is a solution of the time-dependent Schrödinger equation: While the exact choice of time step is arbitrary, it should be small enough to properly resolve the change in the electric microfield and avoid discontinuities in its time sequence. We choose a time step such that an electron with the thermal velocity traverses a fraction (1/20) of the average particle separation: The N used to calculate the average particle separation is the number density of the plasma. Decreasing the time-step factor to as low as 1/100 results in no significant improvements. We therefore choose 1/20 to gain a commensurate speedup in run times while still sufficiently resolving the time sequence of the microfield. The timeevolution operator, U, is then used to compute the timeevolution of the dipole moment operator D in time steps of size Δt: where α and β refer to the upper (initial) and lower (final) states of the transition, respectively. Once we have the time sequence of the dipole moment operator in hand, we use the preferred numerical method (autocorrelation or power spectrum) to compute the Starkbroadened line profile (Anderson 1949;Rosato et al. 2020).

Randomness: Impact Parameter and Maxwellian Velocity Distributions
The artificial conditions of the simulation must be designed to match those of the real physical conditions to the extent possible. The spatial extent of the simulation sphere is dictated by the range at which Coulomb forces are nonnegligible and must be set such that the particle count is high enough to faithfully reproduce the statistical properties of the actual plasma . Finite bounds to the simulation are also imposed by computational limitations. Throughout the duration of the simulation, some fraction of the perturbing particles will inevitably "leave" the bounds of the simulation and must be reinjected. The reinjection scheme is one of the most delicate pieces of the simulation and must be carefully handled to avoid introducing subtle, unphysical behavior. Some of the unanticipated effects can substantially affect the accuracy of the final line shapes.
Regardless of the specifics of the reinjection scheme, the simulation particles must preserve both the impact parameter (b) and Maxwellian velocity (v) probability distribution functions. Historically, simulations have reinjected particles with velocities identical to their exit velocities and with impact parameters resampled from the same narrow bin (or cell) of values containing the original impact parameter of the given particle (Stamm & Voslamber 1979;Gigosos & Cardenoso 1987;Djurović et al. 2009) In this approach, the impact parameter distribution is divided into a large number of bins, and the value is resampled from the same bin or cell that the reinjected particle initially exited the simulation with. In more recent work, Gigosos & Cardeñoso (1996) describe a modified reinjection approach. The impact parameter is resampled according to the same cellbased approach. In the modified approach, velocity is also resampled from the Maxwellian thermal distribution similarly divided into a large number of cells. This avoids correlations between b and v by producing a different bv product while simultaneously avoiding a progressive drifting of the particles toward the edge of the simulation sphere as well as the cooling of the plasma. Tremblay et al. (2020) also choose to implement this approach for reinjection. While this approach does properly preserve both distributions, it does so at the expense of the additional randomization afforded by resampling from the full range of possible values. We employ a different approach to reinjection.
Particles with larger impact parameters and larger velocities will leave the sphere more frequently and will therefore need to be reinjected more frequently. If the impact parameter and velocity for the reinjected particles are simply resampled from the original theoretical distributions, then these distributions will drift to lower values: the distribution will skew toward the longest-lived particles, which are those with smaller impact parameters and velocities. For the initial impact parameters and velocities of the particles at the initialization of each simulation iteration, we sample from the original distributions given by Equations (10) and (12), where the value of the impact parameter is normalized to the radius of the simulation sphere. We modify the conventional approach by sampling the b and v of the reinjected particles from distributions that have been adjusted by particle lifetime to properly preserve the steadystate b and v distribution functions. The modified distributions are given by Equations (11) and (13). This allows us to maximize randomness in the simulation which improves the noise levels of the line shapes and avoids pathological, unphysical behavior.
We define τ, the particle lifetime, as We then adjust the probability distribution functions by dividing the original analytical functions by τ and renormalizing. The equilibrium and reinjected impact parameter distributions are given as follows: Similarly, we obtain the following for the equilibrium and reinjected velocity distributions: Both reinjection distributions have been normalized such that the cumulative probability over the range of possible values sums to 1. We find that reinjecting the particles using the modified distribution functions properly preserves the impact parameter and velocity distributions while maximizing the randomness in the simulation (see Figure 4). To our knowledge, this is the first time such a reinjection scheme has been used in simulation-based line-profile calculations.

Additional Considerations: Microfield Distribution and Poisson Statistics
It is necessary to verify that the changes to the reinjection scheme have not introduced disruptions to the electric microfield distribution-the distribution of net close-range electric field magnitudes felt by the radiator (Hooper 1968(Hooper , 1966. We confirm that the Hooper distribution is still properly preserved (see Figure 5).
We also verified that the particle count in a smaller spherical cavity within the larger simulation sphere is Poissonian. Other authors have used a different approach to achieve the proper Poisson statistics, i.e., a larger "simulation sphere" that feeds particles into a smaller, concentric "calculation sphere" with a radius of 3 Debye lengths. Only particles that lie within the bounds of the calculation sphere contribute to the electric field magnitudes (Hegerfeldt & Kesting 1988;Tremblay et al. 2020). While this does save some computational time, the computational expense of the calculation is largely dominated by the integration of the Schrödinger equation and in particular by the diagonalization of the Hamiltonian. Therefore, we chose to simplify the simulation construction and specify only a single sphere of 5 Debye radii, which is large enough to encompass all particles that make a nonnegligible contribution as well as some that will have a negligible effect on the total microfield. Particle counts in Xenomorph within a smaller sphere of radius 3 Debye lengths converge to a Poissonian distribution. Contributions from all particles within the simulation sphere are considered.

Debye Screening Effects
To fully mimic physical reality, the simulation code would have to take into account particle-to-particle Coulomb interactions for all of the perturbers in the plasma. Given that a minimum of thousands of such simulations would be required to build up the ensemble of electric field time histories that are needed to calculate the averaged profiles, full N-body simulations are still computationally prohibitive. Therefore, most simulation codes instead employ the so-called Debyescreened electric potential. First proposed by Debye & Hückel (1923), most codes now use Debye-screened electric potentials to represent the collective effect of the interparticle interactions. The Debye prescription employs straight-line particle trajectories with a modified electric potential. The Debye screening effect is implemented by modifying the electric field from its original Coulombic form, eZ/r 2 : The Debye length calculation (λ D ) assumes screening from electrons only for a plasma of temperature T, density N s , and charge Ze. It is given by Because the Hamiltonian perturbation depends linearly on the electric potential (V D = − F D · r/Ze), the inclusion of screening results in smaller perturbations leading to narrower line profiles. Stambulchik et al. (2007) evaluated the accuracy of the Debye screening formalism using a full N-body simulation, which accounts for electrostatic interactions between the plasma particles for a hydrogen plasma at T = 1 eV and n e = 10 18 cm −3 . They compared the full N-body results to noninteracting Debye-screened simulation results, where the ions were assumed to be screened by both ions and electrons whereas electrons were screened only by electrons, i.e., 2 D,i D,e l l = . They found that the full N-body simulation approach results in larger electric field magnitudes than the Debye shielding approach. In the full N-body simulation calculation, the FWHM of the line profile was ∼6% greater than that of the screening prescription calculation. This suggests that to properly mimic screening experienced by plasma particles, at least for the particular plasma conditions investigated in that work, the screening may need to be slightly smaller (i.e., a larger screening length) than what is produced by doubly screening ions. The computational burden of the full N-body simulation prohibits exploring the full plasma parameter space. The Debye-screened approach with straight-line trajectories is nevertheless generally taken to be roughly valid at lower densities in the regime of weakly coupled plasmas. The ∼6% difference reported by Stambulchik et al. (2007) suggests this is not an egregious approximation.
For this first generation of line profiles, we chose to adopt the standard approach of singly screened ions and electrons, i.e., both ions and electrons are screened only by electrons, so λ D,i = λ D,e . This likely overestimates the broadness of the line profiles; the singly screened line shape is ∼10% broader than the doubly screened one in our calculations. Figure 6 demonstrates the differences in broadening resulting from different screening prescriptions for Hβ. As expected, including screening from more particles results in smaller electric field magnitudes and consequently narrower profiles. Stambulchik et al. (2007) report a ∼6% difference between their full Nbody line shape and their doubly screened line shape, we note that though performed at similar conditions (T = 1 eV and n e = 10 18 cm −3 ), their calculations were for Hα. Our own Xenomorph calculations for Hα return a similar difference of ∼6.9% between the two screening options: singly screened ions and electrons returns a line profile that is ∼6.9% broader than singly screened electrons and doubly screened ions.
Our goal of evaluating the effects of the new input physics in the Stark broadening calculations relative to the assumptions made by Vidal et al. (1970) and Tremblay & Bergeron (2009) motivated us to use the same Debye screening parameters. It is important to note, however, that the importance of screening should not be underestimated. Indeed, screening has a larger effect on the line shapes than any of the new physics described in Section 3. Though changes in the cores of the profiles tend to have less of an impact on the final emergent stellar spectra than changes in the wings because the WD spectral features approach saturation due to the purity of atmospheric composition, the effects need to be definitively parameterized by including them in model atmospheres. In future work, we intend to more closely investigate which screening prescription yields the best approximation of physical reality. This may require a screening approach that is a function of the plasma parameters.

Choice of Numerical Method
Simulation-based line-shape codes employ one of two numerical methods to calculate the line shape from the timeevolved dipole moment operator. We will henceforth refer to these as "autocorrelation" and "power spectrum." The power spectrum method described in early texts relies on a recasting of classical electrodynamics and the Poynting theorem into QM analogs (Schiff 1955). This method treats the electron as a classical oscillator and calculates the resulting Poynting flux propagating out from the current density produced by the oscillating electron. The classical current density (J(r)) is replaced by the quantum mechanical charge density multiplied by velocity, where the velocity is given by the momentum operator divided by mass.
The final expression for the line shape computed using the power spectrum is given by    where α and β represent the initial and final states of the transition respectively, and ρ α corresponds to the density matrix that contains the level populations of the various states considered in the system, the Boltzmann factors. D(t) is the dipole moment between the two states and ω is the frequency corresponding to the displacement in energy (ΔE) away from the line center. The subscript "av" denotes an ensemble average over the plasma perturbations, which we approximate with a large number of iterations to reduce noise. On the other hand, the dipole autocorrelation method first described by Anderson (1949), and subsequently rederived using an alternative method by Bloom & Margenau (1953), takes advantage of the properties of Fourier transforms to recast and simplify the computation. With the autocorrelation method, the characteristic decorrelation of the dipole moment operator reveals the line shape through the Fourier transform of the dipole moment's autocorrelation function (Baranger 1958b): Here, the capital D indicates this is the matrix corresponding to the operator form of the dipole moments. α and β again represent the initial and final states of the transition, and U is again the time-evolution operator. The subscript "av" again refers to an average over a large number of iterations. Historically, simulation-based line-shape codes have chosen to use the dipole autocorrelation function. Though the rate of convergence for both methods is the same and scales as the inverse of the square root of the number of iterations as expected (  N  1 ), the variance of the power spectrum method is much lower (Rosato et al. 2020). Our implementation of the autocorrelation method compared to the power spectrum method corroborates these findings (see Figure 7). We therefore choose to use the power spectrum method; the advantage of this method has been recognized in recent works (Stambulchik et al. 2007;Rosato et al. 2020).

VCS Shortcomings
Critique of the accuracy of VCS goes back to the 1970 s, beginning with Lee (1971). Lee (1971) pointed out that VCS unifies the impact limit in the line core with the static (one-electron) limit in the wings, but this unification of the two limits does not mean the resulting line shapes are correct across the entire profile. While Lee (1971) has many criticisms of VCS, probably the most important is that the one-electron limit has not been established to produce the correct wing behavior. In reality, the line shape is influenced by the combined effect of ion and electron microfields and particularly as densities increase, simultaneous close approaches of more than one plasma particle become more common. Additionally, Godfrey et al. (1971) pointed out that VCS does not include timeordering, which may cause VCS to underpredict the width, as demonstrated by Roszman (1975). Because TB09 profiles are based on the VCS formalism, these points apply to the TB09 profiles as well.
VCS (and by extension TB09) also suffer from calculational limitations. First, VCS/TB09 do not include the effect of the neighboring states on the transition of interest. For example, if one is calculating an Hβ line shape (n = 4 → 2), then the calculation should also include the n = 3 and n = 5 states; ignoring the adjacent states in this way is known as the noquenching approximation. The no-quenching approximation is suitable at low densities, but at high densities, as perturbative effects increase, this approximation breaks down. Second, for a given transition, VCS/TB09 fails to converge on a solution on the high-density part of the grid. The upper bound on density decreases with decreasing temperature and increasing quantum number n. To fill out the grid, Lemke and TB09 simply copy over the profile calculated at the nearest condition for which the calculation did converge. In some cases, the same profile is copied over multiple adjacent grid points.
In Xenomorph, on the other hand, the calculation is performed by considering the total microfield of both ions and electrons together and does not preclude the possibility that multiple plasma particles can approach the radiator simultaneously. Furthermore, by constructing a chronological history of the time-evolution of the electric field, Xenomorph integrates the time-dependent Schrödinger equation in a way that implicitly includes time-ordering. Finally, Xenomorph does not fail to converge at any set of plasma conditions.

New Physics
Xenomorph incorporates three new pieces of input physics (along with time-ordering which is included by default) that represent improvements over the approximations made by VCS. These are ion dynamics, multipole expansion of the Coulomb interaction potential between radiators and perturbers, and an extended basis set. Each new piece of physics can be turned on independently and studied in isolation to examine their effect on WD atmospheres and their emergent stellar spectra.
We computed line profiles on a grid of number density and temperature for the relevant range of plasma parameters-those spanned by the majority of DA WD atmospheres. Specifically, we calculated Hα, Hβ, and Hγ Balmer-line profiles at temperatures of 5000, 10,000, 20,000, and 40,000 K and at electron densities of 10 15 cm −3 −10 19 cm −3 in steps of 0.5 dex.

Agreement with Previous Theory
We begin by establishing correspondence with the VCS calculations to the extent that this is possible. Under the assumption that the VCS theory has been well vetted over the years and is faithful to the physics that it claims to include, we attempt to reproduce the VCS line profiles using a configuration of Xenomorph that is limited to the same physical assumptions. Specifically, we use a configuration of Xenomorph that includes Debye-screened plasma perturbers, a set of random configurations of particles and corresponding spatial evolution of those particles that reproduces the Hooper microfield distribution, the dipole approximation for the Coulomb interaction potential, and isolated lines that include only the quantum states of the two principal quantum numbers involved in the transition. Large spikes of power are occasionally observed with static ions when one or more ions happens to be placed near the central radiator. We therefore include ion dynamics in these calculations by default to avoid spurious noisy behavior in the far wings of the profiles.
However, there is one fundamental difference between Xenomorph and VCS that makes a perfect correspondence between the two theories impossible. Xenomorph implicitly incorporates the effects of time-ordering on the computed line shapes by virtue of the fact that it constructs a time sequence of the electric field history. On the other hand, the standard VCS theory ignores time-ordering because it relies on an analytic microfield distribution and does not employ the time-ordering operator. Therefore, we do not expect the central structure of line-shape calculations made using the two approaches to agree exactly. As Roszman (1975) points out, the inclusion of timeordering yields differences in the cores of the line profiles. In particular, including time-ordering leads to a lower Hα peak, which is accompanied by an increase in the FWHM. Additionally, time-ordering decreases the relative intensity of the Hβ central dip; in other words, the central dip is filled in relative to VCS calculations which ignore time-ordering. The slight disagreement observed between the VCS and Xenomorph profiles is therefore expected and the agreement with VCS is otherwise excellent in the wings. To the extent that including time-ordering leads to a more physically realistic simulation, we claim that the discrepancies in the core represent another improvement over previous calculations. Figures 8, 9, and 10 show comparisons between the Xenomorph and VCS profiles.

Ion Dynamics
The effect of ion dynamics on line shapes has been studied extensively (Stamm & Voslamber 1979;Gigosos & Cardenoso 1987;Stambulchik & Maron 2006). Many of the earliest simulation line-shape codes were constructed to examine the effect of moving ions (as opposed to static ions) on the profiles to resolve long-standing discrepancies between laboratory measurements and analytical calculations. The inclusion of ion dynamics in the simulation has been shown to yield much better agreement with experimental results in the cores of the profiles (Stamm & Voslamber 1979;Gigosos & Cardenoso 1987). In particular, it has been shown theoretically and verified experimentally that ion dynamics broadens the core of Hα and fills in the central dip of Hβ (Dufty 1969(Dufty , 1970Kelleher & Wiese 1973;Oza et al. 1988). Line-shape calculations used for laboratory diagnostics are particularly concerned with accuracy in this regime. In the context of DA WD model atmospheres, however, the lines in the final stellar spectrum are already so broad that changes to the cores have a negligible overall effect.

Multipole Expansion
Inclusion of higher-order multipole terms in the Coulomb interaction potential between radiator and plasma perturbers also manifests in the cores of the profiles. Gomez et al. (2016) demonstrated that the Hβ core asymmetry was more accurately predicted when the quadrupole term is included in the perturbative potential of the time-dependent Hamiltonian. Closer investigation revealed an error in the calculation of the quadrupole moment in Xenomorph. We reproduce Figure 2 of Gomez et al. (2016) with the revised version of the code and find slightly worse agreement with the data. Additional details regarding the difference will be published in an errata. Xenomorph's predictions for the core asymmetry are smaller than measured values. Similar to the changes made to ion dynamics, the inclusion of this feature does not yield large differences in the features of the synthetic stellar spectra.

Expanded Basis Set
The expanded basis set refers to the expansion of the matrix elements of the relevant operators to include those that correspond to neighboring states. In our investigation, we examine the effect of including one additional set of states belonging to the next principal quantum number above the upper level of the transition in question. For example, we include the n = 4 states in calculations of Hα-the transition from n = 3 to n = 2. Expanding the basis set amounts to a more accurate calculation of the Stark effect. In the perturbation theory approach, this equates to including higher-order Stark terms in the perturbative expansion.

Occupation Probability
In addition to the broadening of spectral lines through the processes mentioned in the previous sections, there is another effect in high-density plasmas that leads to these lines decreasing in strength and merging with the continuum. This is believed to be the result of the dissolution of atomic states due to the high electric fields of nearby ions, which results in the optical electron being removed from the emitting/ absorbing atom. Because this removal occurs for states with negative energies that would otherwise be bound, this effect is often referred to as "ionization potential depression" or "continuum lowering" (Inglis & Teller 1939). To compute this effect on a particular atomic level, one must average over a statistical ensemble of nearby ions to determine the fraction of atoms that are affected by this process. This results in an "effective statistical weight" or an "occupation probability" OP for the given state (e.g., Hummer & Mihalas 1988;Fisher & Maron 2002). In stellar astronomy, the term OP is most frequently used, and we employ it in the remainder of this paper.
In their most recent work on line profiles, Tremblay & Bergeron (2009) used the HM88 formalism (Hummer & Mihalas 1988) with the implementation of Nayfonov et al. (1999) for calculating the OP of states. This formalism assumes that for each state there is a critical, uniform microfield above which the state no longer exists. Then, using a model for the distribution of the microfield values in the plasma, the fraction of atoms that are experiencing microfields less than this value can be calculated, and this is taken to be the OP for the level. This approach uses the critical microfield as a hard cutoff and does not consider gradients in this field.
In contrast, Fisher & Maron (2002, hereafter FM02) employ a criterion based on the distance to the nearest ion. They claim this to be more physical as it is the nearest ion that will be able to "share" or "steal" the electron. They find that a smaller microfield is required to remove an electron because the field increases in the direction of the ion. Thus, FM02 typically calculate OPs that are smaller than those of HM88.
We note that FM02 use a hard cutoff based on the distance to the nearest ion, whereas HM88 use a hard cutoff based on the value of the microfield, assumed to be uniform. In either case, the inclusion of an OP formalism results in lines that are narrower because the broadening effects of the highest microfields are excluded (see Figure 11).
This picture of OP is almost certainly too simplistic. The behavior of level dissolution is not binary with the survival of the level on one side of the microfield (or nearest-neighbor) cutoff and dissolution on the other side. In addition, while the liberated electron may no longer be bound to a single ion, it may be shared between two or more ions in a state with energies and transition probabilities close to their bound-state values (Fisher & Maron 2003). Hence, the spectral features of the state may also not go away suddenly as this criterion is reached. In fact, Fisher & Maron (2003) explicitly state that they do not expect the spectral features to be greatly altered as their nearest-neighbor criterion is reached; rather, these features will begin to weaken only as the distance to the nearest ion further decreases.
Nevertheless, the HM88 formalism is what is currently used for WD atmospheres, and we employ it and FM02 independently in the calculations that follow. For the implementation of the HM88 OP, we limit the closest approach of the ions such that the total electric field does not exceed the maximum microfield value. To enforce this, when the total electric field felt by the radiator exceeds this maximum value, we reverse the radial velocity of the closest ion that is approaching the radiator. In the next time step, if the electric field continues to exceed the threshold value, we reverse the radial velocity of the next-nearest ion that is approaching. We continue this procedure until the electric field falls below the critical microfield value. This successfully produces a sharp cutoff in the microfield distribution, which qualitatively matches the analytical distribution with the cutoff at β crit used in TB09. For the FM02 prescription, the velocity reversals only occur at the one radius corresponding to the nearest ion criterion given Calculations are shown on a linear scale (top row) to highlight core behavior and on a log scale (bottom row) for wing behavior. The leftmost column shows calculations performed at an electron density of n e = 10 16 cm −3 , and the electron density increases to the right. Where labeled, each profile includes only one of the additional pieces of input physics. For example, the calculations labeled "Xeno: idyn" use dynamic ions but do not include higher order multipole terms in the Coulomb potential or an expanded basis set. However, all other profiles have been calculated with ion dynamics turned on. For example, the calculations labeled "Xeno: quad" use both dynamic ions and expand the Coulomb interaction potential to include the quadrupole moment but do not use an expanded basis set. This choice follows from the fact that the prolonged presence of a close ion throughout the duration of a given iteration generates large artificial regions of power in the line profile at the frequency corresponding to the E-field generated by the close ion. Static ions therefore generate spiky profiles in the far wings due to these close ion approaches, and we leave ion dynamics on to avoid this type of pathological behavior in the simulation and in the profiles. Comparisons between VCS and Xenomorph exhibit the expected disagreement in the height of the cores of the profiles. The wings show good agreement between the two theories. by FM02; in practice, this leads to a more severe limitation on large microfields, producing a microfield distribution that more heavily favors smaller field values (see Figure 12). In light of the discussion in Fisher & Maron (2003), we note that this is an upper limit (and likely an overestimate) on the degree of dissolution of these states.
We have performed a cursory comparison of our profiles that include the OP with those of TB09 (see Figure 13). As expected, the FWHM of the profiles depends on the chosen prescription of the OP. The more extreme microfield cutoff of FM02 yields narrower profiles than that of TB09, which uses the milder HM88 criteria. The Xenomorph line shape with the HM88 OP has an FWHM that is ∼5% greater than that of TB09 at the same conditions. The larger width of the Xenomorph profiles is likely due to the inclusion of ion dynamics, and the discrepancy in the heights of the peaks in the core comes from time-ordering. Qualitatively, our profiles differ from those of VCS in the same way that TB09 profiles do; profiles with an OP have narrower cores and a steeper decline in the wings.
While the above OP formalisms are designed to account for state dissolution due to the ions, TB09 employ an additional exponential damping factor to the electron-broadening profiles to account for electron collisions. The physical justification they provide is as follows. Most electron collisions occur so quickly that the electron is unlikely to become unbound. If, on the other hand, the time between successive electron collisions is sufficiently small, due to the time/energy uncertainty principle, there is a finite probability for the electron to have sufficient energy to become unbound. After a sufficient amount of time, the bound electron is assumed to undergo ionization if the plasma is in a steady state (Tremblay & Bergeron 2009). Practically, the correction is applied as follows: 5 where f(ν) refers to the original line shape and f n ¢( ) represents the line shape after the electron-damping correction has been applied. ν c − ν 0 is the difference between the frequency of the line center (ν 0 ) and the frequency corresponding to the energy of the ionization threshold of the lower level in the transition (ν c ).
Using E = hν, we can convert this into an expression written in terms of energy: In this example using Hβ, these would be n up = 4 and n low = 2. We applied this damping factor to our own profiles; however, we were unable to establish correspondence with the published profiles of TB09 (see Figures 14). Applying the same exponential damping factor described in TB09 results in similar wing behavior out to roughly 0.2 eV from line center for Hβ at 10,000 K and 10 18 cm −3 . However, there is a noticeable difference in the cores of the profiles, and the drop in intensity is more dramatic in the far wings for the TB09 profiles than allowed by the prescription provided by Equation (19). Figure 14 shows our attempt to independently reproduce the final electron damped line shapes from the original undamped TB09 line shapes; however, we were again unable to reproduce the final published profiles (P. E. Tremblay, 2021, private communication). The differences in the damped versus undamped profiles seem small. As TB09 point out, the inclusion of electron damping results in changes to log g and T eff of ∼0.5%. Varying the electron-damping prescription can change the fundamental parameters by about half of that amount (±∼0.25%). Compared to the changes of roughly ∼4% in log g and T eff resulting from the inclusion of OP, the effect seems minor. Preliminary investigations reveal only small differences in the final emergent WD spectra when we exclude damping effects from the electron corrections for the Balmer lines. Whether the difference in the profiles presented here will yield significant differences in the emergent stellar spectrum will require careful examination of their effect once propagated through a grid of model atmosphere calculations. We will present a more detailed comparison of the Xenomorph and TB09 profiles in our next publication.
To investigate the importance of OP on the line-profile shapes, we have computed two additional grids of line profiles, one each for the HM88 and FM02 formalisms, for the same three Balmer lines, Hα, Hβ, and Hγ, on the same temperature and electron density grid as before (see Figures 15, 16, and 17). Calculations for both grids also include all three new pieces of input physics described in the preceding sections.
It is clear that continuum lowering/ionization potential depression and line merging are physically observable phenomena, but it is not clear which of the above-mentioned models are accurate, if any. It is therefore unknown which is the best prescription to use: HM88, FM02, or none at all. This is an area of research in atomic physics that is hotly debated, and there is not a consensus as to which models yield the best comparison with experiment. Figure 10. Same as Figure 8 but for Hγ. The VCS code fails to converge at higher densities and lower temperatures as discussed in Section 2.6. This includes Hγ at n e = 10 18 cm −3 and T 10,000 K (right panels). The VCS line profile for T = 20,000 K, n e = 10 18 cm −3 is shown instead because this is the profile copied to lower temperatures in the tables of Lemke.

Comparison to Experiment
We adopt the standard mode for validation of theoretical hydrogen Balmer profiles at high densities by conducting a preliminary and limited comparison of our profiles against the Wiese et al. (1972) experimental data. We use the tabulated TB09 profiles available in the model atmosphere code TLUSTY (e.g., Hubeny & Lanz 1995;Hubeny et al. 2021) to fit for number density with temperature fixed to the value reported in the TB09 paper. We compare these fits against those performed with our own simulation line-profile calculations made with Xenomorph (see Figures 18 and 19). Appropriate subsections of the Wiese et al. (1972) data, which contain a single line, were isolated with the left and right cutoff wavelength points corresponding to the points at which the continuum reaches a minimum between the two adjacent lines. The fits were performed using least-squares minimization. Electron density, scaling factor, and the continuum polynomial were permitted to float as free parameters. Because the theoretical line profiles are area-normalized following convention, we applied a scaling factor to match the intensities of the Wiese data. The continuum in the Wiese data across a given line was fit using a first-order polynomial for Hβ and a second-order polynomial for Hγ and added to the line profiles to bring the continuum levels into agreement.
Xenomorph shows improvements in the fit to the Wiese data in two aspects. First, the Hβ profile fit in the core is vastly improved owing to the combined contribution from all three new pieces of input physics. In particular, ion dynamics decreases the depth of the core and the expanded basis set produces the bulk of the core asymmetry, both of which yield Figure 11. Comparison of Hβ line shapes resulting from different occupation probability prescriptions. Incorporating occupation probability removes the most strongly perturbing ions from the simulation, resulting in narrower line shapes. The FM02 prescription (blue) has a more severe microfield cutoff than that of HM88 (green) at the same plasma conditions. The TB09 line profiles (orange) apply the HM88 prescription to the standard VCS theory (red). The wings of the line profiles calculated with FM02 and HM88 using Xenomorph are almost identical.  . Comparison of Xenomorph profiles calculated with the HM88 OP formalism against TB09. The comparison of the Xenomorph profile against that of TB09 without electron corrections applied represents the fairest comparison because the physical assumptions are as close as can be achieved. The lines labeled "TB09" and "Xeno with HM88 OP" implement the same OP prescription. However, as described above, Xenomorph implicitly includes time-ordering and was calculated with ion dynamics to avoid spurious noisy behavior in the far wings. These differences account for the depressed peak heights and greater FWHM, respectively. The close match between these two lends credence to our simulation-based implementation of occupation probability. improvement in the fit to the data in the core of Hβ. The number densities inferred by independent fits to Hβ and Hγ match much more closely between the Xenomorph fits. While the number densities inferred using the TB09 profiles differ by ∼12% between the two profiles, those inferred using the Xenomorph profiles differ by less than 1%.
It is important to note that the lack of detailed error measurements reported in Wiese et al. (1972) prevents us from making a proper estimation of uncertainties. Using the value of the average ∼1% statistical uncertainty quoted in Wiese yields very high estimates for the value of the reduced chi square. This leads us to conclude that none of the models produces a good fit to the data if the quoted level of uncertainty is to be believed. For these reasons, we choose not to report any uncertainty estimates. We caution against drawing the conclusion that the Xenomorph calculations provide better fits in general. The discrepancy in fit parameters observed for the TB09 calculations might be diminished if the fits were performed just using the wings of the profiles and excluding the cores. For WD atmospheres, it is likely that the differences in the core have a negligible effect on the final effective temperature and g log determinations. Nonetheless, the more accurate reproduction of the cores builds confidence that the extra physics implemented in Xenomorph produces more accurate line profiles in general.

Preliminary Application to Model Atmospheres
We incorporated the new Xenomorph line-shape calculations into DA model atmospheres, which were calculated using the newest released version of TLUSTY. The corresponding highresolution synthetic spectra were also produced using the newest version of SYNSPEC (Hubeny & Lanz 1995). Because the Xenomorph line profiles are not symmetric across the line center when we use an expanded basis set and/or include the quadrupole contribution to the Coulomb interactions, TLUSTY and SYN-SPEC were modified to accept separate tabulations for the profiles to account for the differences between the red and blue sides of the profile. We tabulated the Xenomorph line-profile calculations according to the format described in Lemke (1997). The profiles are area-normalized such that S 1 ò a a D = -¥ ¥ ( ) . We calculated each set of DA WD models at four different grid points for effective temperature and surface gravity: T eff = 17,000 and 20,000 K and log g = 7.0, 8.0, and 9.0 (see Figures 20 and 21). We chose to show the synthetic spectrum calculations for a range that includes only the three lines for which we have calculated new Balmer-line shapes (Hα, Hβ, and Hγ).
In the first set of calculations shown in the left panel of both plots, it is apparent that the model spectra exhibit virtually no differences. Including all three new pieces of physics (ion dynamics, an expanded basis set, and quadrupole interactions) does not produce any appreciable change in the model atmospheres or the corresponding synthetic spectra. On the other hand, including the HM88 prescription of OP as implemented in our simulation context produces narrower lines in the emergent spectra than VCS. Similar to TB09, the differences are negligible for Hα and Hβ. However, the magnitude of the difference is potentially significant for Hγ. Because we decided not to include the prescription for electron damping used in TB09 in our profiles, the drop-off in intensity in the wings of the Xenomorph profiles with HM88 OP is less precipitous. This explains why our calculations differ from the VCS calculations less drastically than those of TB09.
The treatment of OP in the synthetic spectrum calculation calls for further scrutiny. TLUSTY uses the HM88 OP formalism to treat the pseudo-continuum or transitions between bound levels and the set of "dissolved" atomic levels. A self-consistent model atmosphere calculation therefore incorporates line profiles calculated with the HM88 OP formalism as well. The TLUSTY implementation of OP uses cutoffs for each line defining empirical frequency bounds beyond which there is no contribution to the pseudo-continuum from transitions involving the dissolved levels. The location of these empirical frequency bounds is defined somewhat arbitrarily and should be reevaluated to find a firstprinciples definition for their locations.
The model atmosphere calculations demonstrate the viability of using our new simulation-based hydrogen Balmer-line-shape calculations in DA WD model atmosphere and synthetic spectrum calculations. However, we emphasize that these are only preliminary results and, as mentioned previously, there are several issues regarding screening and different OP formalisms that we intend to investigate further. For now, this precludes any mass or temperature determinations using the new profiles.

Conclusions
We present new line-shape calculations of the first three Balmer spectral lines, Hα, Hβ, and Hγ, intended for use in DA Figure 14. Comparison of exponential damping prescription with the published TB09 profile. The published profile (red) corresponds to the VCS theory with occupation probability and electron damping incorporated. When the published electron-damping formula is applied to the original undamped profile with HM88 OP included (solid blue) using the functional form provided in TB09, the result (dashed blue) is different from the published (damped) profile in the core and in the wings. The same procedure applied to the Xenomorph calculation (green) also does not yield good correspondence with the published TB09 profile.  WD model atmosphere calculations. The line shapes were calculated using Xenomorph, a simulation-based, line-shape code. We implement a novel approach to reinjection that uses modified distribution functions that take into account the higher probability that particles with a large impact parameter and/or velocity will preferentially leave the bounds of the simulation more quickly. This allows us to preserve the correct impact parameter and velocity distributions throughout the duration of the simulation while maximizing randomness.
We found that using a different numerical method than most simulation codes yielded much cleaner line profiles with far better noise properties. We choose to use the power spectrum method employed by Stambulchik & Maron (2006) instead of the autocorrelation method to achieve much cleaner lines, eliminating the possibility of encountering any potential spurious effects when interpolating throughout the grids of profiles for Figure 17. Hγ calculations comparing different microfield prescriptions. Again, the VCS line profile at T = 10,000 K, n e = 10 18 cm −3 is copied over from the next higher temperature grid point as discussed in the caption to Figure 10. Figure 18. Fits to high-density Wiese data: Hβ. The Xenomorph fits were performed with profiles that include all of the new physics discussed above but without OP. Fit parameters obtained with profiles calculated including OP did not exhibit significant differences to Xenomorph profiles without OP. Qualitatively, the core asymmetry is more accurately predicted by Xenomorph than by either TB09 or VCS. Figure 19. Fits to high-density Wiese data: Hγ. The best-fit number density obtained using Xenomorph is almost identical to that obtained for Hβ. The shape of the core is again more accurately predicted by Xenomorph than by either TB09 or VCS.
Figure 20. DA model atmospheres calculated at six different grid points corresponding to T eff = 17,000 K and log g = 7, 8, and 9. The left panel compares the VCS line shapes and Xenomorph line shapes with all three new pieces of physics included (but excluding occupation probability). The log g = 7 and 8 models have been shifted up for clarity. The spectra exhibit virtually no differences when we include ion dynamics, expanded basis set, and quadrupole Coulomb interactions. The right panel shows model calculations at the same grid points but with the HM88 occupation probability included in the Xenomorph line shapes. We also overplot models calculated with the TB09 Balmer-line profiles for comparison. The inclusion of occupation probability has the largest impact on the emergent stellar spectra as expected based on the differences observed in the line-shape calculations. The differences between all three sets of models VCS, TB09, and Xenomorph are minimal for Hα and Hβ. However, there are differences for Hγ that might produce appreciable differences in the log g and effective temperature determinations. This can likely be attributed to the fact that the Xenomorph line profiles do not use the prescription for electron damping included in TB09. The specifics of screening and OP formalism call for further scrutiny.  Figure 20 but for T eff = 20,000 K. Again, the log g = 7 and 8 models have been shifted up for clarity. model atmosphere calculations. The cleaner profiles also eliminate the need to perform any smoothing in particularly noisy regions of parameter space or out in the wings of the profiles where noise increases.
Xenomorph incorporates three improvements in the physical treatment of Stark broadening over previous codes: ion dynamics, the inclusion of the quadrupole moment in the Coulomb interactions between radiators and perturbers, and an expanded basis set that includes states from neighboring principal quantum numbers. Additionally, by virtue of the fact that a time history of the electric field is constructed and used to integrate the Schrödinger equation, time-ordering is inherently included in the simulation. Xenomorph also has the flexibility to include a simulation-based approach to OP under two different frameworks that limit either the closest approach of particles to exclude particles within a given radius (FM02) or the maximum value of the electric field (HM88). The HM88 implementation in our simulation was tuned to produce an electric microfield distribution similar to the analytical distribution implemented by Tremblay & Bergeron (2009).
We constructed four separate grids of line profiles corresponding to each of the three new pieces of input physics and one grid that simultaneously incorporates all of the new features. We also performed two sets of calculations with all new features turned on in addition to the HM88 and FM02 OP formalism. Ion dynamics and quadrupole terms in the Coulomb interaction produce minor changes relative to VCS. Calculations performed with an expanded basis set generate small features that mimic the next higher Balmer line in the blue wings of the profiles. Though interesting, it is unlikely that these features will turn out to produce an appreciable difference in the synthetic stellar emergent spectra. The automatic inclusion of time-ordering leads to larger differences than any of the individual pieces of new physics, particularly in the core of Hα; however, this change is also unlikely to make a difference in the final emergent spectra.
We also perform initial investigations of the effect of two different permutations of the Debye screening prescription on Hβ. Screening ions by both ions and electrons, and electrons by electrons only yields an FWHM for Hβ that is ∼15% smaller than screening both ions and electrons only by electrons. The choice of screening prescription has the potential to induce larger changes in the final emergent stellar spectra than any of the three new pieces of input physics described above. Similarly, the effect of OP outweighs that of any of the new input physics, particularly in the wings of the profiles. The particular choices made in implementation affect the final calculations, and the literature offers more than one possible formalism. Future publications will systematically investigate the different choices for OP implementation and their effects on model atmospheres.
Preliminary model atmosphere calculations have been performed using TLUSTY (model atmosphere) and SYNSPEC (synthetic spectrum) at six grid points corresponding to T eff = 17,000 K and 20,000 K and log g = 7.0, 8.0, and 9.0. While these are only preliminary calculations, they demonstrate that Xenomorph line shapes can be successfully used in model atmosphere and synthetic spectrum calculations of DA WDs. The spectra also demonstrate that the three new pieces of physics presented in this paper (ion dynamics, expanded basis set, and quadrupole interactions) are unlikely to have a measurable impact for the lower principal quantum number lines Hα, Hβ, and Hγ on either the synthetic spectra or subsequent mass and temperature determinations. On the other hand, the specifics of the implementation of the OP prescription and screening are likely to have possibly consequential effects. All of the new physical effects incorporated into Xenomorph are likely to be more consequential for lines of higher principal quantum numbers. We will reserve judgment and speculation on whether any of these new physical effects will impact final mass and temperature determinations until we examine higher Balmer lines as well.
The questions and topics investigated in this paper have long been debated in the atomic physics community. Many of them remain open sources of debate. We present this collection of topics in the context of a simulation code developed specifically for the relevant parameter space of DA WD model atmospheres, though the range of possible calculations extends somewhat beyond this. Many of the topics reviewed in this publication require further extensive examination and their implementation in Xenomorph and other line-profile codes will benefit from further refinement, in particular, OP. Despite the remaining outstanding questions, the combination of all of the new physics and time-ordering in the line-shape calculations produces significant improvements in fits to laboratory data and in particular to the standard for the hydrogen Balmer series: the Wiese et al. (1972) data. Though we will continue to investigate OP and refine our implementation, we intend to tabulate grids of line profiles for use in DA WD model atmosphere calculations and for comparison against existing standards for line profiles. These databases will be provided in a forthcoming publication.