Stark-Broadened Profiles for Ionized Helium Lines Using Computer Simulations

1. Introduction

White dwarf stars represent the endpoint of more than 97% of the stars in the Galaxy, and the determination of their properties thus provides a variety of useful information on the star formation history in our galaxy, and on stellar evolution in general [1]. White dwarfs are found in a variety of flavors: The DA stars, whose spectra are dominated by hydrogen lines, comprise about 80% of the white dwarf population, and their atmospheres are presumably hydrogen-rich. The remaining 20% are generally referred to as non-DA stars. These include the DB (DO) stars, whose spectra are dominated by neutral (ionized) helium lines, the DQ stars that show molecular carbon Swan bands, the DZ stars with metallic absorption features, and the DC stars, which are completely featureless. A significant fraction of white dwarf stars enter the cooling sequence with helium-dominated atmospheres (DO-type) but eventually develop hydrogen-dominated atmospheres (DA-type) due to the upward diffusion of residual hydrogen in their envelope – the so-called float-up model (see [2] for a review).

[3] empirically studied this phenomenon by conducting a detailed spectroscopic analysis of 1806 hot white dwarfs ($T_{\mathrm{eff}}\ge 30,000$ K) selected from the Sloan Digital Sky Survey (SDSS). Their study demonstrated that approximately 24% of white dwarfs begin their evolution as DO stars, of which about two-thirds subsequently transform into DA stars. The atmospheric parameters derived from their model atmosphere analysis are summarized in the $logg-T_{\mathrm{eff}}$ diagram shown in Figure 13 of [3], where several peculiar behaviors can be observed. On the one hand, at low effective temperatures, most DA and DO/DB white dwarfs form a continuous sequence, essentially parallel to constant-mass curves and tightly centered around $M\sim 0.55$$M_{\odot }$, in agreement with our understanding of white dwarf evolution. On the other hand, at high temperatures, DA stars with $T_{\mathrm{eff}}\gtrsim 60,000$ K exhibit lower-than-average surface gravities, while DO stars with $T_{\mathrm{eff}}\gtrsim 50,000$ K display higher-than-average values. Consequently, the hydrogen-dominated sequence bends upward toward lower masses, whereas the helium-dominated sequence sharply drops toward higher masses. These features are in clear contradiction with the well-established facts that white dwarfs evolve at a constant mass and that the DA and DB spectral classes at lower temperatures have indeed similar mean masses [4,5]. These mass discrepancies are not unique to the study of [3]. Indeed, based on their compilation of spectroscopic parameters of DO stars from the SDSS, [6][see their Figure 5] also obtained a mass distribution peaking at a larger value of $M\sim 0.675$$M_{\odot }$. They argue that this high-mass peak is physically real, but the striking contrast between the mass distributions of DO and DB stars suggests otherwise. Based on these observations, [3] thus concluded that the spectroscopic masses of very hot SDSS white dwarfs are affected by widespread issues that are not limited to atmospheric models, cooling sequences, or analysis techniques. One possible issue is the line broadening theory of ionized helium (He ii) used in model spectra for white dwarfs.

All model spectra for hot white dwarfs published in the literature rely on the Stark broadening calculations of [7][see also [8]] for the He ii lines, in particular for He ii $\lambda$4686, which is the strongest helium line observed in DO white dwarfs, and DAO stars as well. These calculations are based on the Unified Theory of [9], but including the effects of the ionic charge on the paths of the electrons. Recently, we have used computer simulations to study the Stark broadening of He i lines [10,11], in order to replace the semi-analytical calculations of [12] commonly used in spectroscopic analyses of DB white dwarfs (see, e.g., [13]). It is in this context that we have undertaken a similar investigation of the line broadening of ionized helium using numerical simulations to verify the validity of the profiles of [7] employed in previous spectroscopic studies. In this particular work, we focus only on He ii $\lambda$4686. In contrast with our previous computer simulations [11,14], the perturbed emitting particle here is He ii, which, unlike neutral helium, is a charged particle. In this case, the perturbing particles (electrons and He ii) no longer follow straight-line trajectories, a complex problem that we address in this paper.

We first introduce in Section 2 the fundamental concepts underlying the calculation of Stark-broadened line profiles using computer simulations. In Section 3, we review the method employed for these calculations when the perturbing particles follow straight-line trajectories and extend this approach to hyperbolic trajectories in Section 4. The results of our He ii $\lambda$4686 calculations are presented in Section 5, followed by our conclusions in Section 6.

2. Computer Simulations of Stark-Broadened Line Profiles

2.1. Extending Trivial Molecular Dynamics to Charged Emitters and Hyperbolic Trajectories

The aim of this work is to generalize the computer simulation method developed by several authors [10,15,16,17,18,19,20,21,22] to calculate Stark-broadened profiles for the He ii $\lambda$4686 line. Computer simulation is the state-of-the-art approach for detailed Stark-broadened profile calculations, and tables of line profiles of hydrogen and helium atoms have been calculated using a version of the method called Trivial Molecular Dynamics (which is described below). In this method, the charged particles interacting with the emitting (neutral) ion follow straight-line trajectories at constant velocity, and adapting it to perturbers with more general trajectories in the presence of a charged emitter requires targeted modifications.

[23] were the first to examine the applicability of the simulation method to charged emitters, adopting a hyperbolic orbit for the perturbers. They derived the general equation of the trajectory and proposed a numerical method to determine the position of a particle as a function of time along its trajectory. This approach involves inverting the non-linear equation that defines the hyperbolic eccentric anomaly as a function of time. However, their work does not provide a statistical distribution of phase space coordinates allowing one to generate particles within the simulation volume at initialization, nor to handle incoming particles during the course of the simulation.

Approximately two decades ago, numerical simulation methods diverged into two main families. The first generation of simulation methods, the so-called Trivial Molecular Dynamics (TMD) models treat the perturbers as if they experience no electrostatic interaction other than with the emitter located at the center of the simulation volume, and the electric field at the position of the emitter is computed with a screened-debye potentiel in order to approximate correlation effects in the particle’s motion.

In the case of an uncharged emitter, the trajectory reduces trivially to a straight line within TMD. However, when the central particle carries a charge, two levels of approximation can be employed, depending on the form of the adopted interaction potential. A pure Coulomb potential results in a hyperbolic orbit, whereas the more rigorous alternative of using a screened Debye potential—the same potential employed to compute the electric field in the vicinity of the emitter—leads to a more complex trajectory that cannot be expressed in analytical form. [24] followed this second, i.e., more rigourous, approach to generate line shapes of some H- and He-like ions. While the TMD approach has become increasingly realistic for hydrogenic emitters—benefiting from major advances in the underlying quantum calculations, such as the use of a full Coulomb interaction rather than a multipole-expansion approximation and improved treatment for penetrating collisions [25]—no simulation-based line-profile tables for ionized helium are yet available for astrophysical applications.

The second family of simulations, Full Molecular Dynamics (FMD), reproduces the physical processes more faithfully by adopting a Coulomb interaction between all particles, whether perturbers or the emitter, and by solving the equations of motion simultaneously for the entire system [26]. This constitutes an N-body problem that requires significantly more computing ressources than the TMD approach.

With the aim of generating line-profile tables for ionized helium, we investigate a TMD approach whose level of physical realism lies between the semi-analytical Stark-broadening calculations [7] and computer simulations such as those implemented by [24]. In accordance with the TMD methodology, the motion of the perturbers is determined solely by their interaction with the emitter, which in this work is described by a Coulomb potential rather than a Debye-screened one. This choice is less rigorous physically, but it allows the displacement of a particle along the resulting hyperbolic orbit to be modeled using a strategy similar to that proposed by [23].

2.2. Quantum Treatment of the Emitting Ion

The simulation method developed here preserves the essential elements of the TMD framework, including the correct initialization of perturbers within the simulation volume and the handling of new particles entering the domain over time. It differs, however, by replacing the traditionally employed autocorrelation function with the power spectrum, as proposed by [27] and [22].

The power spectrum is defined as

$$P\left(\omega \right)\propto {\omega }^4\underset{T\to +\infty }{lim}{\displaystyle \frac{1}{T}}\sum _{ab}{\rho }_a{〈|{\int }_{-T/2}^{T/2}dt\langle b\left|\mathbf{d}\left(t\right)\right|a\rangle e^{i\omega t}{|}^2〉}_{\mathrm{av}}$$

(1)

where $|a\rangle$ and $|b\rangle$ are the upper and lower states of the emitting ion for the considered transition, ${\rho }_a$ is the statistical weight of the upper states, and $\mathbf{d}\left(t\right)$ is the dipole operator of the emitting ion in the Heisenberg picture. The above equation includes a thermal average over all possible configurations of moving perturbers.

As usual, the slowly varying ${\omega }^4$ term is excluded from the detailed calculations, yielding the line profile

$$L\left(\omega \right)\propto \underset{T\to +\infty }{lim}{\displaystyle \frac{1}{T}}\sum _{ab}{\rho }_a{〈|{\int }_{-T/2}^{T/2}dt\langle b\left|\mathbf{d}\left(t\right)\right|a\rangle e^{i\omega t}{|}^2〉}_{\mathrm{av}}$$

(2)

which is then normalized to unit area. In the simulation approach, the thermal average is approximated by an average over a large but finite set of N dynamic configurations. The line profile is thus computed as the mean of N individual profiles $L_i\left(\omega \right)$, each corresponding to a specific dynamic configuration i,

$$L\left(\omega \right)={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}{L}_i\left(\omega \right)$$

(3)

where

$$L_i\left(\omega \right)\propto \underset{T\to +\infty }{lim}{\displaystyle \frac{1}{T}}\sum _{ab}{\rho }_a|{\int }_{-T/2}^{T/2}dt\langle b\left|\mathbf{d}\left(t\right)\right|a\rangle e^{i\omega t}{|}^2$$

(4)

is computed from the ith temporal sequence of dipole moments.

The no-quenching approximation is used in this work. Consequently, the matrix element of the dipole operator is expressed as an expansion over two disjoint sets of upper and lower states, $|a^{\prime }\rangle$ and $|b^{\prime }\rangle$, which include all states sharing the same principal quantum number (and spin) as the $|a\rangle$ and $|b\rangle$ states involved in the transition:

$$\langle b\left|\mathbf{d}\left(t\right)\right|a\rangle =\sum _{a^{\prime }{b}^{\prime }}\langle b|U_b^{†}\left(t\right)|b^{\prime }\rangle \langle b^{\prime }\left|\mathbf{d}\right|a^{\prime }\rangle \langle a^{\prime }|U_a\left(t\right)|a\rangle$$

(5)

where $\mathbf{d}$ is now the dipole operator in the Schrödinger picture, and $U_{a,b}$ are the evolution operators projected onto the upper and lower subspaces of states, respectively.

The projected evolution operators satisfy

$$\begin{array}{ccc}\hfill i\hslash {\displaystyle \frac{d{U}_a\left(t\right)}{dt}}& =& H_a\left(t\right)U_a\left(t\right)\hfill \\ \hfill i\hslash {\displaystyle \frac{d{U}_b\left(t\right)}{dt}}& =& H_b\left(t\right)U_b\left(t\right)\hfill \end{array}$$

(6)

where $H_{a,b}$ are projections of the Hamiltonian

$$H=H_0+V\left(t\right)$$

(7)

onto the two relevant subspaces. Here, $H_0$ is the Hamiltonian of the unperturbed ion, while $V\left(t\right)$ represents the time-dependent interaction between the emitting ion and the moving perturbers. In this work, we consider only the first nonvanishing term in the multipole expansion of the interaction—namely, the dipole term:

$$V\left(t\right)=e\mathbf{F}\left(t\right)·\mathbf{R}$$

(8)

where $\mathbf{R}$ is the position operator of the radiating electron of the emitting ion.

In each temporal sequence, the ion emitter remains at the center of a sphere, referred to as the simulation volume, with a radius on the order of the Debye radius. The analysis of the variations in the electric field $\mathbf{F}\left(t\right)$ experienced by the emitter focuses on contributions from the charged particles—ions and electrons—within the simulation volume. The moving perturbers are treated as independent classical particles, following trajectories determined by their charge and that of the central emitter.

The line profile $L_i\left(\omega \right)$ associated with the ith simulation volume is computed in three steps. First, a time sequence of electric fields, ${\mathbf{F}}_i\left(t_k\right)$, at the origin of the simulation volume is generated based on the spatial distribution of moving perturbers for a series of $N_t$ discrete times, $t_k\equiv k{\Delta }_t$. Next, the Schrödinger equation (6) is numerically integrated. The integration scheme ensures that the evolution operator matrices remain unitary, and the time step ${\Delta }_t$ is chosen to be sufficiently small to resolve variations in the field ([18], [19], [28], [10], [22]). Finally, the time sequence of dipole moment matrix elements (equation 5) is substituted into equation (4), where the integral over t is performed for different frequencies.

A key challenge of the computer simulation method is ensuring that, across all simulation volumes, the statistical properties of the system remain unchanged both over time and across parallel simulations. For straight-line trajectories at constant velocity, this is achieved by choosing an appropriate coordinate system in six-dimensional phase space and deriving the corresponding joint statistical distribution. This distribution is used for both the initialization and reinjection phases of the simulation. Particle positions are then generated within this coordinate system, called the impact coordinate system [29], using a uniform random number generator. However, this system is not suitable for Stark broadening calculations involving charged emitters.

To address this limitation, the proposed method introduces a coordinate system, referred to as the hyperbolic impact coordinate system. This framework ensures that generated particles follow hyperbolic trajectories while preserving the statistical properties of the system over time. Section 3 examines the case of an uncharged central particle with perturbers following straight-line trajectories, following the standard approach described in the literature on TMD. This involves replacing Cartesian or radial coordinates ($\mathbf{r},\mathbf{v}$) in phase space with a six-dimensional coordinate vector that accounts for the trajectory shape. As will be shown, this second coordinate system (the hyperbolic impact system) is particularly well-suited for initializing particles within a sphere and randomly generating additional particles that enter and cross the sphere during the simulation. The complete set of equations governing particle generation with hyperbolic trajectories is detailed in Section 4.

3. Straight-Line Trajectory: Standard Method of the TMD

This section presents the particle generating process relevant for the well-known case of a neutral emitter and charged perturbers following straight-line trajectories at constant velocity in the TMD model. This topic has been frequently addressed in the literature ([15], [16], [17], [18], [28], [10], [22]), but a complete review of the steps involved in constructing the straight line trajectory and deriving the relevant joint statistical distributions in the impact coordinate system is necessary to familiarize oneself with the notation and more complex calculations required for hyperbolic trajectories.

3.1. Parameterization of a Straight-Line Trajectory

For a neutral emitter positioned at the center of the simulation volume, each perturbing particle follows a straight-line trajectory given by

$$\mathbf{r}\left(t\right)=\mathbf{b}+\mathbf{v}(t-{\tau }_b)$$

(9)

where $\mathbf{v}$ is the constant velocity vector of the particle, ${\tau }_b$ is the time of closest approach, and $\mathbf{b}$ is the trajectory’s closest point to the origin. The vector $\mathbf{b}$ is perpendicular to $\mathbf{v}$, and its norm, b, represents the impact parameter. The velocity vector $\mathbf{v}$ is defined as

$$\mathbf{v}\equiv v{\mathbf{v}}_{\mathbf{3}}(\theta ,\varphi ),$$

(10)

where v is the magnitude of $\mathbf{v}$ and ${\mathbf{v}}_{\mathbf{3}}$ is a unit vector. The angles ($\theta ,\varphi$) define the orthonormal basis:

$$\begin{array}{ccc}\hfill {\mathbf{v}}_{\mathbf{1}}& =& {(cos\varphi cos\theta ,sin\varphi cos\theta ,-sin\theta )}^{\mathrm{T}}\hfill \\ \hfill {\mathbf{v}}_{\mathbf{2}}& =& {(-sin\varphi ,cos\varphi ,0)}^{\mathrm{T}}\hfill \\ \hfill {\mathbf{v}}_{\mathbf{3}}& =& {(cos\varphi sin\theta ,sin\varphi sin\theta ,cos\theta )}^{\mathrm{T}},\hfill \end{array}$$

(11)

where $0\le \varphi <2\pi$, $0\le \theta \le \pi$, and T denotes the transpose operation. When $\theta =\varphi =0$, these vectors align with the positive x, y, and z axes, respectively. The perpendicularity condition between $\mathbf{b}$ and $\mathbf{v}$ implies that $\mathbf{b}$ can be expressed as a linear combination of ${\mathbf{v}}_{\mathbf{1}}$ and ${\mathbf{v}}_{\mathbf{2}}$:

$$\mathbf{b}=b(cos\alpha {\mathbf{v}}_{\mathbf{1}}+sin\alpha {\mathbf{v}}_{\mathbf{2}}),$$

(12)

introducing a new angle $\alpha$ in the range $0\le \alpha <2\pi$.

The six parameters $\mathbf{x}\equiv (\theta ,\varphi ,\alpha ,v,b,{\tau }_b)$ define a valid coordinate system, as each point ($\mathbf{r},\mathbf{v}$) in phase space corresponds uniquely to a set of values $\mathbf{x}$ for any constant t. The only exception occurs at the origin ($\left|\mathbf{r}\right|=0$), where the angular coordinates become degenerate. However, such points are not considered valid positions for perturbing particles in the simulation.

The probability density function $p\left(\mathbf{r},\mathbf{v}\right)$, associated with an infinitesimal phase-space volume ${\mathbf{d}}^{\mathbf{3}}\mathbf{r}{\mathbf{d}}^{\mathbf{3}}\mathbf{v}$, must align with statistical mechanics:

$$p(\mathbf{r},\mathbf{v})\propto e^{-E\left(\mathbf{r},\mathbf{v}\right)/kT},$$

(13)

where E is the total energy of the particle (kinetic plus potential), k is Boltzmann’s constant, and T is the temperature.

To express this probability density in terms of the new coordinate system, we use the jacobian J of the transformation:

$$p\left(\mathbf{x}\right)\propto e^{-E\left(\mathbf{x}\right)/kT}J\left(\mathbf{x}\right),$$

(14)

where E and

$$J\equiv \left|{\displaystyle \frac{{\partial }^{\mathbf{3}}\mathbf{r}{\partial }^{\mathbf{3}}\mathbf{v}}{{\partial }^{\mathbf{6}}\mathbf{x}}}\right|$$

(15)

are functions of $\mathbf{x}$. The jacobian is given by the determinant of the jacobian matrix:

$$\mathbf{J}\equiv \left(\begin{array}{cc}{\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial v}}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial v}}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial \theta }}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial \theta }}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial \varphi }}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial \varphi }}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial {\tau }_b}}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial {\tau }_b}}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial b}}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial b}}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial \alpha }}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial \alpha }}\end{array}\right).$$

(16)

Because each partial derivative can be expanded as a linear combination of ${\mathbf{v}}_{\mathbf{1},\mathbf{2},\mathbf{3}}$, the jacobian matrix factorizes into two simpler matrices:

$$\mathbf{J}=\left(\begin{array}{cccccc}0& 0& (t-{\tau }_b)& 0& 0& 1\\ v(t-{\tau }_b)& 0& -bcos\alpha & v& 0& 0\\ -bsin\alpha & bcos\alpha cos\theta +vsin\theta (t-{\tau }_b)& 0& 0& vsin\theta & 0\\ 0& 0& -v& 0& 0& 0\\ cos\alpha & sin\alpha & 0& 0& 0& 0\\ -bsin\alpha & bcos\alpha & 0& 0& 0& 0\end{array}\right)\left(\begin{array}{cc}{\mathbf{v}}_{\mathbf{1}}^{\mathrm{T}}& 0\\ {\mathbf{v}}_{\mathbf{2}}^{\mathrm{T}}& 0\\ {\mathbf{v}}_{\mathbf{3}}^{\mathrm{T}}& 0\\ 0& {\mathbf{v}}_{\mathbf{1}}^{\mathrm{T}}\\ 0& {\mathbf{v}}_{\mathbf{2}}^{\mathrm{T}}\\ 0& {\mathbf{v}}_{\mathbf{3}}^{\mathrm{T}}\end{array}\right).$$

(17)

Because the determinant of the second matrix is unity, the jacobian simplifies to the determinant of the first matrix:

$$J=b{v}^{3}sin\theta .$$

(18)

Substituting this expression into equation (14) yields the probability density in the new coordinate system:

$$p(\theta ,\varphi ,\alpha ,v,b,{\tau }_b)\propto b{v}^3{e}^{-v^2/v_T^2}sin\theta ,$$

(19)

where $v_T=\sqrt{2kT/m}$ is the thermal velocity, with m denoting the reduced mass. Because there is no interaction with the neutral emitter, the Boltzmann factor accounts only for the kinetic energy.

This probability distribution serves as the foundation for deriving the joint statistical distributions necessary for particle initialization and the injection of new particles during the simulation.

3.2. Statistical Distribution at Initialization of the Simulation

The purpose of the initialization procedure is to generate a set of perturbers that follow the joint statistical distribution of coordinates described in equation (19), while ensuring their presence within the simulation volume, which is a sphere of radius R.

Because

$$\mathbf{r}(t=0)=\mathbf{b}-\mathbf{v}{\tau }_b,$$

(20)

the constraint that particles reside inside the sphere implies

$$|\mathbf{b}-\mathbf{v}{\tau }_b{|}^2=b^2+v^2{\tau }_b^2\le R^2,$$

(21)

or, equivalently,

$$|{\tau }_b|\le t_{\max },$$

(22)

where

$$t_{\max }\equiv {\displaystyle \frac{1}{v}}\sqrt{R^2-b^2}$$

(23)

represents half the particle’s traversal time. This constraint on ${\tau }_b$ introduces statistical dependence among the random variables $(v,b,{\tau }_b)$.

By combining the joint statistical distribution (equation 19) with the constraint above, the probability distributions for each coordinate can be derived. The distributions for the three angular coordinates are straightforward:

$$\begin{array}{ccc}\hfill p\left(\theta \right)& =& {\displaystyle \frac{1}{2}}sin\theta \hfill \\ \hfill p\left(\varphi \right)& =& {\displaystyle \frac{1}{2\pi }}\hfill \\ \hfill p\left(\alpha \right)& =& {\displaystyle \frac{1}{2\pi }}\hfill \end{array}$$

(24)

which are properly normalized. Substituting these distributions back into equation (19) gives the joint probability distribution for the remaining three coordinates, which factorizes as

$$p(v,b,{\tau }_b)=p\left(v\right)p\left(b\right|v)p\left({\tau }_b\right|b,v)\propto b{v}^3{e}^{-v^2/v_T^2}.$$

(25)

Because this expression has no explicit dependence on ${\tau }_b$, the distribution of ${\tau }_b$ (given b and v) is uniform within its allowed range:

$$p\left({\tau }_b\right|b,v)={\displaystyle \frac{1}{2t_{\max }}},\mathrm{if}\left|{\tau }_b\right|\le t_{\max }.$$

(26)

Substituting this result into equation (25) yields the joint distribution of v and b:

$$\begin{array}{ccc}\hfill p\left(v\right)p\left(b\right|v)& \propto & b{v}^3{e}^{-v^2/v_T^2}\times t_{\max }\hfill \\ & \propto & b\sqrt{R^2-b^2}{v}^2{e}^{-v^2/v_T^2}.\hfill \end{array}$$

(27)

From this, the normalized probability distribution of b (given v) is

$$p\left(b\right|v)={\displaystyle \frac{3}{R^3}}b\sqrt{R^2-b^2}\equiv p\left(b\right),$$

(28)

indicating that b and v are statistically independent. Finally, substituting this result into equation (27) and normalizing again gives the Maxwellian distribution of velocity:

$$p\left(v\right)=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{1}{v_T^3}}{v}^2{e}^{-v^2/v_T^2}.$$

(29)

The particle generation process at initialization follows these steps. For each particle: (1) the angles $\varphi$ and $\alpha$ are drawn from a uniform distribution between 0 and $2\pi$, (2) the angle $\theta$, the impact parameter b, and the velocity v are sampled from the distributions given by equations (24), (28), and (29), respectively, and (3) the coordinate ${\tau }_b$ is drawn from a uniform distribution over $\pm t_{\max }$, where $t_{\max }$ is determined by b and v (equation 23).

Once these coordinates have been generated for a set of $N_p$ particles of each type p, the simulation begins. At each time step $t_k\equiv k{\Delta }_t$, where k is an integer and ${\Delta }_t$ is the time step, a particle with coordinates $(\theta ,\varphi ,\alpha ,v,b,{\tau }_b)$ is moved according to

$$\begin{array}{ccc}\hfill \mathbf{r}\left(t_k\right)& =& b(cos\alpha {\mathbf{v}}_{\mathbf{1}}+sin\alpha {\mathbf{v}}_{\mathbf{2}})+v{\mathbf{v}}_{\mathbf{3}}(t_k-{\tau }_b)\hfill \\ & =& \mathbf{r}\left(t_{k-1}\right)+v{\mathbf{v}}_{\mathbf{3}}{\Delta }_t.\hfill \end{array}$$

(30)

As the simulation progresses, some particles will exit the simulation sphere and must be replaced with new incoming particles of the same type.

3.3. Incoming Particles Generating Process

Historically, finding a reliable procedure for generating incoming particles has been challenging. The initial approach involved creating a new particle with the same b and v values as the outgoing particle, which led to nonphysical behaviors. [18] proposed an alternative method based on discrete bins, ensuring that the system’s statistical properties remained constant over time. However, this approach imposed a rigid distribution for b and v, eliminating random fluctuations. Moreover, it was only applicable when these coordinates were statistically independent— a condition satisfied for straight-line trajectories but, as will be shown, not for hyperbolic ones.

[22] introduced a more general procedure that accounts for statistical dependence between b and v while allowing random fluctuations in both coordinates. The statistical distribution of the incoming particles’ b and v coordinates must match that of the outgoing particles. However, particles with longer traversal times through the sphere (i.e., $2t_{\max }$) leave less frequently. Consequently, the distribution $p_{\mathrm{in}}(v,b)$ of incoming particles must be proportional to the distribution $p(v,b)$ within the sphere (equation 27), divided by the traversal time:

$$p_{\mathrm{in}}(v,b)=p_{\mathrm{in}}\left(v\right)p_{\mathrm{in}}\left(b\right|v)\propto p\left(v\right)p\left(b\right|v)\times {\displaystyle \frac{1}{t_{\max }}},$$

(31)

which yields

$$\begin{array}{ccc}\hfill p_{\mathrm{in}}\left(b\right|v)& =& {\displaystyle \frac{2}{R^2}}b\equiv p_{\mathrm{in}}\left(b\right)\hfill \\ \hfill p_{\mathrm{in}}\left(v\right)& =& {\displaystyle \frac{2}{v_T^4}}{v}^3{e}^{-v^2/v_T^2}\hfill \end{array}$$

(32)

after normalization; b and v remain independent variables in the reinjection process for straight-line trajectories, though the reasoning extends to statistically dependent cases.

When a particle of type p exits the simulation volume at time $t_k$, a new incoming particle of the same type is generated. The angles $\theta ,\varphi ,\alpha$ and the coordinates b and v are drawn from their respective distributions (equations 24 and 32). The coordinate ${\tau }_b$ is then determined by requiring that the incoming particle be located on the surface of the simulation volume, or within a thin shell adjacent to it, at time $t_k$.

4. Hyperbolic Trajectory: Generalization of the Standard Method

The equations used to generate particles with hyperbolic trajectories are derived in this section for the TMD model. We draw inspiration from the proven method for straight-line trajectories, for which the six new coordinates serve two roles: they constitute a coordinate system in phase space (at each time t), and they parameterize all possible trajectories as a function of time. The first step is therefore to uniquely define a hyperbolic trajectory with six parameters, and show they adequately cover the entire phase space volume relevant for the simulation. The jacobian of the coordinate transformation from the Cartesian coordinates $\left(\mathbf{r},\mathbf{v}\right)$ in phase space to the set of new coordinates is determined, along with the relevant statistical distributions used to generate particles at initialization and incoming particles.

4.1. Geometry of the Hyperbola

A hyperbola centered at the origin with its principal axes aligned along the x- and y-coordinates is given by:

$$\mathbf{r}=acosh\left(u\right)\mathbf{x}+bsinh\left(u\right)\mathbf{y},$$

(33)

where u is the hyperbolic eccentric anomaly, and a and b are the semi-major and semi-minor axes, respectively. In gravitational trajectory studies, the convention $a<0$ is often used to unify elliptical and hyperbolic orbit equations. However, in this work, we adopt $a>0$ to avoid confusion with the signs that appear in subsequent equations.

Equation (33) must be modified for simulation purposes by shifting the trajectory so that the source, originally at a focus, is moved to the origin:

$$\begin{array}{ccc}\hfill \mathbf{r}& =& acosh\left(u\right)\mathbf{x}+bsinh\left(u\right)\mathbf{y}+\lambda \sqrt{a^2+b^2}\mathbf{x}\hfill \\ & =& a\left((cosh\left(u\right)+\lambda e)\mathbf{x}+\sqrt{e^2-1}sinh\left(u\right)\mathbf{y}\right),\hfill \end{array}$$

(34)

where $\lambda$ represents the sign of the electric force ($-1$ for attraction and $+1$ for repulsion), and e is the eccentricity:

$$e={\displaystyle \frac{\sqrt{a^2+b^2}}{a}}>1.$$

(35)

Given that $b/a=\sqrt{e^2-1}$, the radial distance from the origin follows as:

$$r=a(ecosh(u)+\lambda ),$$

(36)

which is always positive because $ecosh\left(u\right)>1$. The minimum distance, $\sigma$, occurs at $u=0$, where $cosh\left(u\right)$ is minimized:

$$\sigma =a(e+\lambda ).$$

(37)

A key characteristic of a hyperbolic trajectory is its asymptote and associated impact parameter. As $u\to \pm \infty$, the trajectory approaches two straight lines with a minimum distance to the origin of b, meaning the asymptotic impact parameter is precisely the hyperbola’s b parameter.

To obtain the most general trajectory, three rotations are applied. First, the trajectory in equation (34) undergoes a rotation by an angle $\beta$ in the x-y plane. Then, the normal $\mathbf{z}$ of the orbital plane is rotated to align with:

$$\begin{array}{ccc}\hfill {\mathbf{w}}_{\mathbf{3}}& =& {(cos\varphi sin\theta ,sin\varphi sin\theta ,cos\theta )}^{\mathrm{T}},\hfill \end{array}$$

(38)

where $0\le \varphi <2\pi$ and $0\le \theta \le \pi$. After applying these transformations, the general trajectory becomes:

$$\mathbf{r}=a\left((cosh\left(u\right)+\lambda e){\mathbf{w}}_{\mathbf{1}}+\sqrt{e^2-1}sinh\left(u\right){\mathbf{w}}_{\mathbf{2}}\right),$$

(39)

where:

$$\begin{array}{ccc}\hfill {\mathbf{w}}_{\mathbf{1}}& =& cos\beta {\mathbf{v}}_{\mathbf{1}}+sin\beta {\mathbf{v}}_{\mathbf{2}},\hfill \\ \hfill {\mathbf{w}}_{\mathbf{2}}& =& -sin\beta {\mathbf{v}}_{\mathbf{1}}+cos\beta {\mathbf{v}}_{\mathbf{2}},\hfill \end{array}$$

(40)

and:

$$\begin{array}{ccc}\hfill {\mathbf{v}}_{\mathbf{1}}& =& {(cos\varphi cos\theta ,sin\varphi cos\theta ,-sin\theta )}^{\mathrm{T}},\hfill \\ \hfill {\mathbf{v}}_{\mathbf{2}}& =& {(-sin\varphi ,cos\varphi ,0)}^{\mathrm{T}}.\hfill \end{array}$$

(41)

The vectors ${\mathbf{w}}_{\mathbf{1},\mathbf{2},\mathbf{3}}$ form an orthonormal basis that coincides with the x, y, and z axes when $\theta =\varphi =\beta =0$.

4.2. Dynamics of the Hyperbolic Trajectory

The equations in the previous section described the trajectory’s geometry without considering its dynamics. We now incorporate energy conservation to determine how the implicit parameter u depends on time t.

In a coordinate system where the emitter (charge $q_e$) is stationary, the total energy of a perturbing particle with reduced mass m, charge $q_p$, and relative velocity v is the sum of kinetic and potential energies:

$$\begin{array}{ccc}\hfill E& =& {\displaystyle \frac{m{v}^2}{2}}+{\displaystyle \frac{q_e{q}_p}{r}}\hfill \\ & =& m\left({\displaystyle \frac{v^2}{2}}+\lambda {\displaystyle \frac{\kappa }{r}}\right),\hfill \end{array}$$

(42)

where $\lambda$ represents the sign of $q_e{q}_p$, determining whether the interaction is attractive or repulsive. The interaction intensity parameter $\kappa$ is always positive and defined as:

$$\kappa \equiv {\displaystyle \frac{|q_e{q}_p|}{m}}.$$

(43)

Energy is of course conserved and takes the same value (as a function of a) regardless of whether the force is attractive or repulsive:

$$E={\displaystyle \frac{m\kappa }{2a}}.$$

(44)

Combining equations (42) and (44) gives the relation between distance r and velocity v:

$$\kappa \left({\displaystyle \frac{1}{a}}-\lambda {\displaystyle \frac{2}{r}}\right)=v^2,$$

(45)

from which we obtain the velocity at infinity:

$$v_{\infty }^2=\underset{r\to \infty }{lim}\kappa \left({\displaystyle \frac{1}{a}}-\lambda {\displaystyle \frac{2}{r}}\right)={\displaystyle \frac{\kappa }{a}}.$$

(46)

The velocity vector $\mathbf{v}$ is given by the time derivative of equation (39):

$$\mathbf{v}\equiv {\displaystyle \frac{d\mathbf{r}}{dt}}=a{\displaystyle \frac{du}{dt}}(sinh\left(u\right){\mathbf{w}}_{\mathbf{1}}+\sqrt{e^2-1}cosh\left(u\right){\mathbf{w}}_{\mathbf{2}}),$$

(47)

where the standard derivatives of hyperbolic functions have been used. The time derivative of u follows from equations (36), (45), and (47):

$${\displaystyle \frac{du}{dt}}=\sqrt{{\displaystyle \frac{\kappa }{a^3}}}{\displaystyle \frac{1}{(ecosh(u)+\lambda )}},$$

(48)

where we take the positive square root to ensure that u increases with time. Integrating this equation gives the time dependence of u:

$$esinh\left(u\right)+\lambda u=\sqrt{{\displaystyle \frac{\kappa }{a^3}}}(t-{\tau }_b),$$

(49)

where ${\tau }_b$ is an integration constant, identified as the time of closest approach (when $u=0$). Finally, substituting this result into equation (47) yields an alternative expression for the velocity vector:

$$\mathbf{v}=\sqrt{{\displaystyle \frac{\kappa }{a}}}{\displaystyle \frac{1}{ecosh\left(u\right)+\lambda }}(sinh\left(u\right){\mathbf{w}}_{\mathbf{1}}+\sqrt{e^2-1}cosh\left(u\right){\mathbf{w}}_{\mathbf{2}}).$$

(50)

4.3. Jacobian of the Hyperbolic Impact Coordinate Transformation

The next step is to identify the set of phase-space coordinates that simultaneously represent every accessible point in phase space (at any given time t) and describe every hyperbolic trajectory of a perturber. Several parameters can characterize a hyperbolic trajectory in addition to the three angular parameters: the particle’s energy, the velocity at infinity, the eccentricity, the semi-major and the semi-minor axes, among others. We impose the constraint that the particle generation method must converge to the case of straight-line trajectories when the interaction intensity parameter $\kappa$ vanishes.

The most natural generalization of the velocity coordinate v of a straight-line path is the velocity at infinity, $v_{\infty }$. Two possible choices exist for the coordinate corresponding to the impact parameter in straight-line trajectories: the minimum distance to the origin, $\sigma$, and the asymptotic impact parameter, b. It turns out that the latter naturally leads to a joint statistical distribution that converges to the case of straight-line trajectories. The sixth coordinate adopted is the time of closest approach, ${\tau }_b$, which corresponds to its equivalent in the context of straight-line trajectories.

It can be shown that the six proposed parameters, $\mathbf{x}\equiv (\theta ,\varphi ,\alpha ,v_{\infty },b,{\tau }_b)$, provide a valid coordinate system for the hyperbolic trajectory because each point $\left(\mathbf{r},\mathbf{v}\right)$ in phase space corresponds to a unique solution $\mathbf{x}$ at any constant value of t. The only exception is the set of points with negative energy. However, these points were not considered as admissible perturber positions within the numerical simulation.

As discussed for straight-line trajectories (equation 14), the probability density of the particle’s position in the new coordinates is given by the product of a Boltzmann factor and the jacobian of the coordinate transformation. For hyperbolic trajectories, the jacobian matrix $\mathbf{J}$ is expressed as

$$\mathbf{J}\equiv \left(\begin{array}{cc}{\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial v_{\infty }}}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial v_{\infty }}}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial b}}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial b}}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial {\tau }_b}}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial {\tau }_b}}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial \theta }}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial \theta }}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial \varphi }}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial \varphi }}\\ {\displaystyle \frac{\partial {\mathbf{r}}^{\mathrm{T}}}{\partial \beta }}& {\displaystyle \frac{\partial {\mathbf{v}}^{\mathrm{T}}}{\partial \beta }}\end{array}\right).$$

(51)

After systematic though tedious algebra (see Appendix A), evaluating the determinant of $\mathbf{J}$ leads to the jacobian:

$$J=b{v}_{\infty }^{3}sin\theta ,$$

(52)

which, when substituted into equation (14), yields the probability distribution for the particle’s position in the new coordinate system:

$$p(\theta ,\varphi ,\alpha ,v_{\infty },b,{\tau }_b)\propto b{v}_{\infty }^3{e}^{-v_{\infty }^2/v_T^2}sin\theta .$$

(53)

Surprisingly, this probability distribution retains the same form for both straight-line trajectories at constant velocity v and hyperbolic trajectories. This holds as long as a correspondence is made between v and $v_{\infty }$, and between the impact parameter b in both cases. This expression serves as the foundation for deriving the joint statistical distributions needed for both the initialization of particle generation and the injection of new particles during the simulation.

4.4. The Probability Distribution at Initialization

The goal of this section is to determine the joint statistical distribution of the six coordinates ($\theta ,\varphi ,\alpha ,v_{\infty },b,{\tau }_b$) required for generating particles at the initialization ($t=0$) of the simulation. The statistical distribution (equation 53) is subject to the constraint that particles must be inside the sphere of radius R. Additionally, at each step, we allow the interaction intensity parameter $\kappa$ to converge to 0, ensuring that in the absence of interaction (i.e., a neutral emitter), the solution corresponds to a straight-line trajectory. It is important to note that the minimum distance $\sigma$ of the hyperbolic trajectory from the origin is not identical to the impact parameter b.

Given the isotropy of the problem, as inferred from the angular dependence of the statistical distribution (equation 53), the generation of the three angles follows equation (24), derived for straight-line trajectories. However, the generation of the remaining three coordinates deviates from this simple case.

The minimum distance to the origin is expressed in terms of the variables $v_{\infty },b$, and ${\tau }_b$:

$$\sigma =a(e+\lambda )=\sqrt{a^2+b^2}+\lambda a,$$

(54)

where e is the eccentricity (equation 35), $\lambda =\pm 1$, and $a=\kappa /v_{\infty }^2$. The constraint requiring particles to be within the simulation sphere implies that

$$\sigma =\sqrt{a^2+b^2}+\lambda a\le R,$$

(55)

which leads to the condition

$$R\ge 2amax(0,\lambda ).$$

(56)

This result has distinct implications depending on the sign of $\lambda$. When $\lambda =-1$, the condition reduces to the trivial result $R\ge 0$. However, when $\lambda =+1$, it imposes a constraint on $v_{\infty }$, determining the minimum velocity required for a particle to enter the sphere:

$$v_{\infty }\ge \sqrt{{\displaystyle \frac{2\kappa }{R}}}.$$

(57)

Thus, low-energy particles cannot enter the sphere due to the repulsive interaction. Combining these results for both values of $\lambda$, the minimum velocity (at infinity) required for a particle to enter the sphere is

$$v_{\min }\equiv \mathrm{Max}\left(0,\lambda \sqrt{{\displaystyle \frac{2\kappa }{R}}}\right).$$

(58)

Following the same approach used for straight-line trajectories, the distribution of the three dependent variables is factored into a product of distributions:

$$p(v_{\infty },b,{\tau }_b)=p\left(v_{\infty }\right)p\left(b\right|v_{\infty })p\left({\tau }_b\right|b,v_{\infty })\propto e^{-v_{\infty }^2/v_T^2}{v}_{\infty }^{3}b.$$

(59)

The constraint that particles must be present in the sphere (at initialization) is expressed by an inequality involving the implicit variable u of the hyperbolic trajectory, starting from equation (36):

$$\begin{array}{ccc}\hfill r\left(u\right)=a(ecosh(u)+\lambda )\le R\to \left|u\right|& \le & {cosh}^{-1}\left({\displaystyle \frac{R-\lambda a}{\sqrt{a^2+b^2}}}\right)\equiv u_{\max }.\hfill \end{array}$$

(60)

This result is converted into a constraint on the time of closest approach by using the relation between u and ${\tau }_b$ (equation 49) for the particular case $t=0$:

$${\tau }_b=t-\sqrt{{\displaystyle \frac{a^3}{\kappa }}}(esinh\left(u\right)+\lambda u)\to {\displaystyle \frac{v}{a}}{\tau }_b=-(esinh\left(u\right)+\lambda u).$$

(61)

The constraint $\left|u\right|\le u_{\max }$ implies that $|{\tau }_b|\le t_{\max }$, with

$$t_{\max }={\displaystyle \frac{1}{v}}\left(\sqrt{{(R-\lambda a)}^2-\left(a^2+b^2\right)}+\lambda a{cosh}^{-1}\left({\displaystyle \frac{R-\lambda a}{\sqrt{a^2+b^2}}}\right)\right),$$

(62)

which is a function of $v_{\infty }$ and b. For a neutral emitter, the interaction intensity parameter $\kappa$ (and therefore a) vanishes, reducing the half-crossing time to the solution for a straight-line path:

$$t_{\max }={\displaystyle \frac{1}{v}}\sqrt{R^2-b^2}.$$

Because $p(v_{\infty },b,{\tau }_b)$ has no explicit dependence on ${\tau }_b$, the probability distribution $p\left({\tau }_b\right|b,v_{\infty })$ is uniform within its domain:

$$p\left({\tau }_b\right|b,v_{\infty })={\displaystyle \frac{1}{2t_{\max }}}\mathrm{with}\left|{\tau }_b\right|<t_{\max }.$$

(63)

Substituting this result back into equation (59) gives the joint distribution of $v_{\infty }$ and b:

$$\begin{array}{ccc}\hfill p\left(v_{\infty }\right)p\left(b\right|v_{\infty })& \propto & e^{-v_{\infty }^2/v_T^2}{v}_{\infty }^{3}b\times t_{\max }\hfill \\ & =& b\left(\sqrt{{(R-\lambda a)}^2-(a^2+b^2)}+\lambda a{cosh}^{-1}\left({\displaystyle \frac{R-\lambda a}{\sqrt{a^2+b^2}}}\right)\right)e^{-v_{\infty }^2/v_T^2}{v}_{\infty }^2.\hfill \\ & & \hfill \end{array}$$

(64)

The requirement that the square root remains real imposes a constraint on the impact parameter b for a particle to enter the sphere:

$$\begin{array}{ccc}\hfill {(R-\lambda a)}^2-(a^2+b^2)& \ge & 0\hfill \\ \hfill \to b& \le & \sqrt{R(R-2\lambda a)}\equiv b_{\max },\hfill \end{array}$$

(65)

where $b_{\max }$ depends on $v_{\infty }$. The distribution $p\left(b\right|v_{\infty })$ can be determined by isolating terms involving b:

$$p\left(b\right|v_{\infty })\propto b\left(\sqrt{{(R-\lambda a)}^2-(a^2+b^2)}+\lambda a{cosh}^{-1}\left({\displaystyle \frac{R-\lambda a}{\sqrt{a^2+b^2}}}\right)\right).$$

(66)

To normalize this distribution, we compute the normalization constant $K_B$, given by the following integral over the allowed range of b, separated into two terms, $I_1$ and $I_2$:

$$\begin{array}{ccc}\hfill K_B& \equiv & {\int }_0^{b_{\max }}dbb(\sqrt{{(R-\lambda a)}^2-(a^2+b^2)}+\lambda a{cosh}^{-1}\left({\displaystyle \frac{R-\lambda a}{\sqrt{a^2+b^2}}}\right)\hfill \\ & =& I_1+\lambda a{I}_2,\hfill \end{array}$$

(67)

where both terms have analytical solutions:

$$\begin{array}{ccc}\hfill I_1& =& {\displaystyle \frac{1}{3}}{\left({(R-\lambda a)}^2-a^2\right)}^{{\displaystyle \frac{3}{2}}}\hfill \\ \hfill I_2& =& {\displaystyle \frac{1}{2}}\left((R-\lambda a)\sqrt{{(R-\lambda a)}^2-a^2}-a^2{cosh}^{-1}\left({\displaystyle \frac{R-\lambda a}{a}}\right)\right).\hfill \end{array}$$

(68)

Thus, the properly normalized distribution is

$$p\left(b\right|v_{\infty })=K_B^{-1}\times b\left(\sqrt{{(R-\lambda a)}^2-(a^2+b^2)}+\lambda a{cosh}^{-1}\left({\displaystyle \frac{R-\lambda a}{\sqrt{a^2+b^2}}}\right)\right),$$

(69)

which implies that b and $v_{\infty }$ are not statistically independent, in contrast to the case of a straight-line trajectory.

For a neutral emitter, where there is no electrical interaction ($\kappa =0,a=0$), the integral simplifies to $I_1={\displaystyle \frac{1}{3}}{R}^3$ and $\lambda a{I}_2=0$. Consequently, the resulting distribution reduces to

$$\begin{array}{ccc}\hfill p\left(b\right|v)& =& 3{\displaystyle \frac{b}{R^3}}\sqrt{R^2-b^2},\hfill \end{array}$$

which recovers the well-known result for a straight-line trajectory, demonstrating its independence from v. Substituting equation (69) into (64) and normalizing gives the velocity distribution at infinity:

$$\begin{array}{ccc}\hfill p\left(v_{\infty }\right)& =& K_V^{-1}{e}^{-v_{\infty }^2/v_T^2}{v}^2\times \left(I_1(R,\lambda a)+\lambda a{I}_2(R,\lambda a)\right),\hfill \end{array}$$

(70)

where $a=k/v_{\infty }^2$, and the normalization constant $K_V$ is determined through numerical integration over the valid velocity range:

$$K_V\equiv {\int }_{v_{\min }}^{\infty }p\left(v_{\infty }\right)d{v}_{\infty }={\int }_{v_{\min }}^{\infty }dv{e}^{-v^2/v_T^2}{v}^2\times (I_1\left(R,\lambda {\displaystyle \frac{\kappa }{v^2}})+\lambda {\displaystyle \frac{\kappa }{v^2}}{I}_2\left(R,\lambda {\displaystyle \frac{\kappa }{v^2}}\right)\right),$$

(71)

where $v_{\min }$ is given by equation (58).

The initialization of the particle generation process follows these steps. For each particle, the angles $\varphi$ and $\alpha$ are drawn from a uniform distribution between 0 and $2\pi$. The angle $\theta$ and velocity $v_{\infty }$ are then sampled from equations (24) and (70), respectively. Once $v_{\infty }$ is known, the impact parameter b is determined using equation (69). Finally, given $v_{\infty }$ and b, the time of closest approach ${\tau }_b$ is drawn from a uniform distribution between $\pm t_{\max }$, as defined in equation (62). The numerical details of this random generation process are discussed in Appendix B.

Once the coordinates for a set of $N_p$ particles of each type p are generated, the simulation begins. At each time step $t_k\equiv k{\Delta }_t$, where k is an integer and ${\Delta }_t$ is the time increment, the particle positions are updated according to:

$$\begin{array}{ccc}\hfill \mathbf{r}\left(t_k\right)& =& a\left((cosh\left(u_k\right)+\lambda e){\mathbf{w}}_{\mathbf{1}}+\sqrt{e^2-1}sinh\left(u_k\right){\mathbf{w}}_{\mathbf{2}}\right),\hfill \end{array}$$

(72)

where

$$\begin{array}{ccc}\hfill a& =& {\displaystyle \frac{\kappa }{v_{\infty }^2}},\hfill \\ \hfill e& =& {\displaystyle \frac{\sqrt{a^2+b^2}}{a}}.\hfill \end{array}$$

(73)

The vectors ${\mathbf{w}}_{\mathbf{1}}$ and ${\mathbf{w}}_{\mathbf{2}}$ depend on the three angles $(\theta ,\varphi ,\alpha )$, while $u_k$ is the solution to the nonlinear equation:

$$\begin{array}{ccc}\hfill esinh\left(u_k\right)+\lambda u_k& =& \sqrt{{\displaystyle \frac{\kappa }{a^3}}}(t_k-{\tau }_b).\hfill \end{array}$$

(74)

Thus, computing the particle’s position at time $t_k$ requires solving for $u_k$, a procedure detailed in Appendix C.

4.5. The Probability Distribution of Incoming Particles

The probability distribution of incoming particles is obtained by extending the method of [22] to hyperbolic trajectories. The distribution $p_{\mathrm{in}}(v_{\infty },b)$ of incoming particles must be proportional to the distribution $p(v_{\infty },b)$ of particles inside the sphere (equation 64), divided by the crossing time:

$$p_{\mathrm{in}}(v_{\infty },b)=p_{\mathrm{in}}\left(v_{\infty }\right)p_{\mathrm{in}}\left(b\right|v_{\infty })\propto p\left(v_{\infty }\right)p\left(b\right|v_{\infty })\times {\displaystyle \frac{1}{t_{\max }}}.$$

(75)

The result retains the same form as in the case of straight-line trajectories:

$$p_{\mathrm{in}}\left(v_{\infty }\right)p_{\mathrm{in}}\left(b\right|v_{\infty })\propto e^{-v_{\infty }^2/v_T^2}{v}_{\infty }^{3}b,$$

(76)

which implies that

$$p_{\mathrm{in}}\left(b\right|v_{\infty })\propto b.$$

(77)

The normalization constant is obtained by solving the integral

$${\int }_0^{b_{\max }}bdb={\displaystyle \frac{1}{2}}{b}_{\max }^2,$$

(78)

where $b_{\max }$ is defined in equation (65). The normalized distribution is then

$$p_{\mathrm{in}}\left(b\right|v_{\infty })={\displaystyle \frac{2b}{b_{\max }^2}}={\displaystyle \frac{2}{R}}{\displaystyle \frac{b}{(R-2\lambda a)}},$$

(79)

with $a=\kappa /v_{\infty }^2$. Unlike in the case of straight-line trajectories, b and $v_{\infty }$ are not independent variables.

Substituting this result into equation (76) yields the (unnormalized) velocity distribution at infinity for particles entering the sphere:

$$p_{\mathrm{in}}\left(v_{\infty }\right)\propto e^{-v_{\infty }^2/v_T^2}{v}_{\infty }\left(v_{\infty }^2-2\lambda {\displaystyle \frac{\kappa }{R}}\right),$$

(80)

which is valid for $v_{\infty }>v_{\min }$, as defined in equation (58). Normalizing this expression gives the velocity distributions for incoming particles:

$$\begin{array}{ccc}\hfill \lambda =-1,v_{\infty }\ge 0:p_{\mathrm{in}}\left(v_{\infty }\right)& =& {\displaystyle \frac{2}{v_T^4}}{e}^{-(v_{\infty }^2/v_T^2)}{v}_{\infty }^3{\displaystyle \frac{1+2\frac{\kappa }{R{v}_{\infty }^2}}{1+2\frac{\kappa }{R{v}_T^2}}},\hfill \\ \hfill \lambda =+1,v_{\infty }\ge v_{\min }:p_{\mathrm{in}}\left(v_{\infty }\right)& =& {\displaystyle \frac{2}{v_T^4}}{e}^{-(v_{\infty }^2-v_{\min }^2)/v_T^2}\times v_{\infty }\left(v_{\infty }^2-v_{\min }^2\right),\hfill \end{array}$$

(81)

where $v_{\min }=\sqrt{{\displaystyle \frac{2\kappa }{R}}}$ for $\lambda =+1$, and $v_T$ is the thermal velocity. The result for straight-line trajectories is recovered when $\kappa =0$.

When a particle of type p exits the simulation volume at time $t_k$, a new incoming particle must be generated. The angles $\theta ,\varphi ,\alpha$ and the coordinates b and $v_{\infty }$ are drawn from their respective distributions (equations 24, 79, and 81). The numerical details of this random generation process are discussed in Appendix D. The remaining coordinate ${\tau }_b$ is defined so that, at time $t_k$, the particle is located precisely on the surface of the simulation volume.

5. Results

5.1. Statistical Distribution of the Coordinates

The first validation tests of the simulation consist of examining the statistical distributions of the perturber coordinates—ions and electrons—and of the resulting electric-field intensities, both at initialization and as a function of time.

Figure 1 illustrates the distribution of velocities at infinity, averaged over all simulation volumes, at initialization, at mid-time, and at the end of the simulation, shown separately for ions and electrons. The distributions follow the Maxwell distribution, except in the case of ions, where particles with insufficient energy to enter the simulation volume are excluded because of the repulsive interaction with the emitting ion, as defined by equation (58). A similarly good agreement between the theoretical and numerical distributions is obtained for the impact parameter b, as shown in Figure 2.

The distribution of the electric-field intensity produced by the ions is consistent with the Hooper distribution, derived for the environment of an ion with charge $Z=+1$ [30], as illustrated in Figure 3 for three characteristic times. The electron distribution exhibits a more pronounced contribution in the wings, resulting in a reduced peak intensity. This behavior arises from the attractive interaction, which favors the presence of electrons near the emitting ion.

A diagnostic that is sensitive both to the random generation of the particle coordinates and to the history of particle motion is the pair-correlation function $g\left(r\right)$, also known as the radial distribution function. It quantifies the probability of finding a charged perturber at a distance r from the emitting radiator relative to the case of a uniform distribution. In this way, $g\left(r\right)$ measures the local modification of the particle density around the emitter induced by electrostatic interactions.

Assuming an interaction potential $U\left(r\right)$ between the perturbers of a given type and the emitter, $g\left(r\right)$ is defined as

$$g\left(r\right)=e^{-U\left(r\right)/kT}.$$

(82)

Because the perturbers’ trajectories are generated assuming the full Coulomb interaction with the emitter rather than a screened Debye potential,

$$g\left(r\right)=\exp \left(-Z{e}^2/kT\right)$$

(83)

where $Z=\pm 1$ for the He ii ions and electrons. Figure 4 compares the generated $g\left(r\right)$ function for ions and electrons, at three characteristic times, with its theoretical behaviour (equation 83).

The fluctuation level of the generated functions is relatively high, especially at short distances, where the number of particles contained within a thin shell centered at a distance r has a large standard deviation relative to its expected mean. Nevertheless, the three experimental curves remain consistent with one another within the fluctuation range, indicating that the particle-reinjection procedure successfully reproduces the statistical distribution of the initial particle positions.

Moreover, the behavior of the experimental curves is qualitatively consistent with the theoretical prediction, in the sense that the particle densities of both ions and electrons deviate from the uniform density at short distances in a manner consistent with their respective charges. There remains, however, a systematic bias between the theoretical and experimental curves, most notably for the electrons, although the significance of this trend is difficult to assess because of the high noise level in the experimental data.

A way to reduce the noise level—at the cost of losing information at specific distances r—is to evaluate the statistical distribution of particles contained within a volume of radius r. We define a function $G\left(r\right)$ as the number of particles inside this radius, normalized by the number expected for a uniform distribution:

$$G\left(r\right)\equiv {\displaystyle \frac{{\int }_0^{r}dxg\left(x\right)x^2}{{\int }_0^{r}dx{x}^2}}={\displaystyle \frac{1}{3r^3}}{\int }_0^{r}dxg\left(x\right)x^2.$$

(84)

The $G\left(r\right)$ functions obtained for ions and electrons, shown in Figure 5, exhibit a high degree of smoothness, making it possible to clearly identify the differences between the experimental and theoretical distributions.

Figure 5. Statistical distribution of $G\left(r\right)$ for ions and electrons for a temperature of 20,000 K and an electronic density of ${10}^{17}$ cm⁻³ at three time-step indices $k=t_k/{\Delta }_t$ of 0, $N_t/2$, and $N_t-1$. The distance r is normalized by the Debye radius $R_d$. The theoretical ion and electron functions are shown in orange and cyan, respectively, whereas the corresponding generated functions are depicted in red and blue.

Figure 5. Statistical distribution of $G\left(r\right)$ for ions and electrons for a temperature of 20,000 K and an electronic density of ${10}^{17}$ cm⁻³ at three time-step indices $k=t_k/{\Delta }_t$ of 0, $N_t/2$, and $N_t-1$. The distance r is normalized by the Debye radius $R_d$. The theoretical ion and electron functions are shown in orange and cyan, respectively, whereas the corresponding generated functions are depicted in red and blue.

A slight discrepancy therefore remains, particularly for electrons. This effect may result from the minimum cut-off imposed on the point of closest approach $\sigma$ of the trajectories. Excluding trajectories with small $\sigma$ could lead to a lower-than-expected number of electrons located at short distances r from the emitter. This outcome is reminiscent of a similar result reported by [25] for electrons by excluding "bound" electrons from the statistics.

The preceding tests indicate that the statistical distributions agree with the expectations of statistical mechanics from the very beginning of the simulation, aside from a small, still not fully understood deviation in the spatial distribution of the electrons. As these distributions remain stable over time, no thermalization period is required to ensure that the dynamical behavior of the perturbers is statistically stationary.

5.2. Linear Stark Effect

An additional validity test of the simulation environment is the numerical verification that the linear Stark broadening characteristic of hydrogenic ions under the no–quenching approximation is correctly recovered. Owing to the degeneracy of energy levels sharing the same principal quantum number n, together with the particular analytical form of the dipole–transition matrix elements between states of identical n, the Stark effect is known to be linear for hydrogenic ions when mixing between states of different principal quantum numbers is neglected.

Assuming a static uniform electric field of fixed magnitude F directed along the z axis, the eigenstates of the corresponding Hamiltonian for states with principal quantum number n yield Stark components whose energy shifts can be described using the parabolic quantum numbers $(n_1,n_2,m)$, with the constraint

$$n=n_1+n_2+\left|m\right|+1.$$

(85)

Defining $q=n_1-n_2$, the energy shift of a state $|n,n_1,n_2,m\rangle$ varies linearly with the electric field F and is given by

$$\Delta E_{n,q}={\displaystyle \frac{3}{2Z_c}}nqe{a}_{0}F,$$

(86)

where $Z_c$ is the core charge of the emitter. The allowed values of q are the integers in the range $\pm (n-1)$.

Therefore, for the Heii $\lambda 4686$ line ($Z_c=2$), the Stark shifts of the upper and lower levels are

$$\Delta E_{4,q_4}=3q_{4}e{a}_{0}F,q_4=0,\pm 1,\pm 2,\pm 3,$$

$$\Delta E_{3,q_3}={\displaystyle \frac{9}{4}}{q}_{3}e{a}_{0}F,q_3=0,\pm 1,\pm 2.$$

The behavior of the perturbed energies predicted by the simulation is shown in Figure 6 for the upper and lower levels of the Heii $\lambda 4686$ transition as a function of the electric field. All Stark levels display the characteristic linear behavior, and there are 7 and 5 distinct Stark components for the upper and lower states, respectively, as expected.

5.3. Hydrogen Line Profiles

A final test used to validate the numerical simulation as a whole consists in computing a hydrogen line profile and comparing it with profiles from the literature. To this end, we generated 20,000 electric-field sequences for a temperature of $20,000\mathrm{K}$ and an electron density of $N_e={10}^{16}$ cm⁻³, for which the perturbers—electrons and protons—follow the straight-line trajectories described in Section 3. A second series of sequences was also generated in which the ions were kept static.

Figure 7 compares our computed H$\beta$$\lambda$4861 line profiles —obtained using both static and dynamic ion treatments— and the profiles predicted by the semi-analytical models of [31,32], in which the ionic perturbers are assumed to behave quasi-statically. Differences are observed in the line core and far out in the wings between our profile including ion dynamics and those reported in the literature. Ion dynamics contribute directly only to the core of the line profile, but the resulting change in shape affects the total area of the profile prior to normalization; the normalization procedure, which enforces unit area, can therefore produce a global modification of the profile intensity, even in the far wings.

Imposing static ions in the computer simulation significantly improves the agreement. The residual noise in the far wings of our static-ion profile might arise from an insufficient number of electric-field sequences, which does not adequately sample the range of static ion configurations.

5.4. He II Line Profiles

The current work focuses on two first illustrative calculations. Although a full grid of He ii $\lambda 4686$ line profiles spanning a broad range of thermodynamic conditions is planned, only two profiles have been computed so far, for a temperature of $20,000\mathrm{K}$ and electron densities of $N_e={10}^{16}$ and 10¹⁷ cm⁻³. The perturber population consists of electrons and He ii ions in equal proportions.

The principal parameters of the model are the time step ${\Delta }_t$, the total number of time steps $N_t$, the number of independent time sequences of electric fields used to compute the final profile, the minimal allowed distance from the emitter, and the radius R of the simulation volume—conveniently expressed as a multiple $K_D$ of the Debye radius

$$R=K_D\sqrt{{\displaystyle \frac{kT}{4\pi (N_e+N_{\mathrm{He}\mathrm{II}})e^2}}}.$$

For the specific profiles discussed in this section, $20,000$ independent temporal sequences of electric fields were generated, each consisting of $N_t=50,000$ values sampled at regular times $t_k=k{\Delta }_t$. The time step was set to a fraction $ϵ=0.02$ of the characteristic timescale associated with the temporal variation of the microfield,

$${\Delta }_t=ϵ{\displaystyle \frac{r_0}{v_T}},$$

where $r_0$ denotes the typical interparticle distance in a plasma of electron density $N_e$, and $v_T$ is the electron thermal velocity computed using the reduced mass ${\mu }_e$. The resulting duration, $N_t{\Delta }_t$, ensures that the autocorrelation function of the emitter’s dipole moment has vanished, fulfilling one of the parameter-validity criteria. The radius of the simulation volume was chosen such that it contains 600 and 250 perturbers of each species, corresponding to $K_D=3.0$ and $K_D\approx 3.35$, for $N_e={10}^{16}$ and 10¹⁷ cm⁻³, respectively. Trajectories for which the minimum approach distance $\sigma$ (equation 37) falls below the characteristic radius of the upper level — $n^2{a}_0/Z_c$, with $n=4$, $Z_c=2$ the charge of the nucleus, and $a_0$ the Bohr radius, were excluded from the simulation and replaced with admissible trajectories.

Available profile tables for the He ii $\lambda 4686$ line are scarce and originate from calculations performed approximately three decades ago. The profiles published by different authors were all obtained using semi-analytical theories assuming quasi-static ions. However, as in the present work, the electrons were assumed to follow hyperbolic trajectories.

Figure 8 compares the profiles computed by these authors with the one generated in the present work for a temperature of $20,000$K and an electron density of $N_{\mathrm{e}}={10}^{17}{\mathrm{cm}}^{-3}$. The profile from [33] results from our own implementation of the equations provided by the authors, which represent an approximate formulation of the theory developed by Griem in the 1960s. This profile exhibits a discontinuity arising from the transition between the electron–impact regime in the line core and the effects of the Lewis cut-off [34] in the wings. The profile from [31] is based on the theory developed by [35]. [7] rely on the Unified Theory, which models the transition from the impact to the quasi-static regime for electrons while keeping the ions quasi-static. Effects of ion dynamics are therefore neglected in all of these profiles.

Significant differences are observed both among the published profiles and between each of them and our own. The discrepancies with our profile cannot be attributed solely to ion dynamics, which are included only in our simulation approach. Indeed, we generated a second profile in which the ions were held fixed (with only the electrons in motion), also shown in the figure. The contribution of ion dynamics is negligible at this density, consistent with our results for neutral helium [11]. Similar discrepancies are observed at a density of $N_e={10}^{16}$ cm⁻³, as illustrated in Figure 9.

6. Discussion

The Stark-broadened profiles we have computed for the He ii $\lambda$4686 transition differ in several important respects from the profiles currently used in all published model-atmosphere analyses of DO white dwarfs. These differences arise directly from the improved physical treatment implemented in our simulations—most notably the explicit inclusion of a large set of electric-field time sequences, which allows accurate integration of the ion’s time-dependent Hamiltonian.

At this stage, our simulation environment successfully passes all validation tests, including the statistical distribution of particle coordinates, the distribution of the electric field, the expected linear behavior of the Stark component energies, and the calculation of an H$\beta$$\lambda 4861$ line profile at a specific temperature and density that agrees with the literature. Our attempts to explain the remaining discrepancies have been unsuccessful. This contrasts with our results for the neutral-helium profiles [11], which are nearly identical whether obtained from the semi-analytical theory or from the numerical simulation, except for the contribution of ion dynamics at lower densities.

The discrepancies we find between our new profiles and the earlier ones, particularly those of [7], are not merely of academic interest. As discussed in the Introduction, spectroscopic analyses of hot DO white dwarfs have long exhibited puzzling and physically implausible trends: inferred surface gravities that increase sharply at high effective temperatures, mass distributions that differ markedly from those of cooler DB stars, and systematic offsets relative to the DA sequence. These anomalies have resisted explanation despite extensive efforts to improve atmospheric models, cooling sequences, and fitting techniques. Given that He ii $\lambda$4686 is the dominant diagnostic feature in DO spectra, it is reasonable to suspect that inaccuracies in the underlying Stark broadening theory may be a significant—perhaps even primary—contributor to these long-standing issues.

Our exploratory calculations therefore raise the intriguing possibility that the spectroscopic mass discrepancies observed in DO white dwarfs may originate, at least in part, from limitations in the Stark profiles that have been used for decades. However, it is essential to emphasize that this conclusion cannot yet be drawn with confidence. The present work focuses exclusively on a single transition, and only for a limited set of plasma conditions. To assess the true impact of our improved treatment, a complete grid of He ii line profiles must be computed and incorporated into full non-LTE model atmospheres. Only then will it be possible to perform a systematic reanalysis of DO spectra and determine whether the new profiles lead to more physically consistent atmospheric parameters.

In this sense, the results presented here should be viewed as a first step—an encouraging indication that the standard Stark broadening theory for ionized helium may require revision, but not yet a definitive solution to the DO mass problem. The next phase of this project will involve extending our simulations to all relevant He ii transitions and generating a comprehensive set of profiles suitable for model-atmosphere calculations. Once this grid is available, we will be able to test directly whether our improved Stark profiles yield more coherent mass distributions and eliminate the anomalous trends observed in previous spectroscopic studies. A similar study performed with neutral-helium profiles and a large sample of DB white dwarfs has indeed revealed avenues for improving our profile calculations with the semi-analytical theory [36].

If these forthcoming analyses confirm that the new profiles resolve the inconsistencies highlighted in the Introduction, this would represent a significant advance in our understanding of DO white dwarfs and, more broadly, in the modeling of dense helium plasmas. Until then, the results presented here should be interpreted as a strong motivation for further work rather than a final verdict.