# A Novel Solution for Resonant Scattering Using Self-Consistent Boundary
Conditions

B. Connor McClellan Shane Davis Phil Arras Department of Astronomy, University
of Virginia, Charlottesville, VA 22904, USA

###### Abstract

We present two novel additions to the semi-analytic solution of Lyman $\alpha$
(Ly$\alpha$) radiative transfer in spherical geometry: (1) implementation of
the correct boundary condition for a steady source, and (2) solution of the
time-dependent problem for an impulsive source. For the steady-state problem,
the solution can be represented as a sum of two terms: a previously-known
analytic solution of the equation with mean intensity $J=0$ at the surface,
and a novel, semi-analytic solution which enforces the correct boundary
condition of zero-ingoing intensity at the surface. This solution is compared
to that of the Monte Carlo method, which is valid at arbitrary optical depth.
It is shown that the size of the correction is of order unity when the
spectral peaks approach the Doppler core and decreases slowly with line center
optical depth, specifically as $(a\tau_{0})^{-1/3}$, which may explain
discrepancies seen in previous studies. For the impulsive problem, the time,
spatial, and frequency dependence of the solution are expressed using an
eigenfunction expansion in order to characterize the escape time distribution
and emergent spectra of photons. It is shown that the lowest-order
eigenfrequency agrees well with the decay rate found in the Monte Carlo escape
time distribution at sufficiently large line-center optical depths. The
characterization of the escape-time distribution highlights the potential for
a Monte Carlo acceleration method, which would sample photon escape properties
from distributions rather than calculating every photon scattering, thereby
reducing computational demand.

††software: numpy (Harris et al., 2020), scipy (Virtanen et al., 2020), Coblis
- Color Blindness Simulator (color-blindness.com)

##  1 Introduction

Given the abundance of hydrogen in the universe, the Lyman $\alpha$
(Ly$\alpha$) line is an important component of radiation fields in a wide
range of astrophysical settings. Ly$\alpha$ radiation transport is an active
area of research in the study of planets, stars, galaxies, and cosmology
(Dijkstra, 2019). An example application motivating our work is the role of
Ly$\alpha$ in planetary atmospheres. The outer layers of the atmosphere are
central to a planet’s evolution, since they can shelter the lower atmosphere
from high energy radiation as well as regulate the escape of gas into space.
There are two sources of Ly$\alpha$: the star, and recombinations in the
planet’s atmosphere. Ly$\alpha$ may ionize atoms and dissociate molecules, as
well as exert pressure forces that drive an outflow (Bourrier et al., 2018).
Ly$\alpha$ can also excite H atoms to the 2p state, creating a population of
Balmer-line absorbers that can be observed via transmission spectroscopy
(Huang et al., 2017; Yan et al., 2021). Due to the low gas densities in the
upper atmosphere, collisional de-excitation and broadening are of secondary
importance and Ly$\alpha$ may undergo “resonant scattering”.

Hubble Space Telescope (HST) observations with the STIS have found large
Ly$\alpha$ transit depths around a handful of exoplanets (Vidal-Madjar et al.,
2003; Lecavelier des Etangs et al., 2012; Ehrenreich et al., 2012, 2015;
Bourrier et al., 2017a, b, c, 2018; Waalkes et al., 2019; Lavie et al., 2019;
García Muñoz et al., 2020; Bourrier et al., 2021). These observations have
revealed a population of atoms extending out to distances of order a few
planetary radii or more for several planets around bright, nearby stars,
motivating a study of the physics of Ly$\alpha$ interactions with the H atom
population. The transition from the atomic to the molecular layer in these hot
upper atmospheres may take place at pressures of order ${\sim}\ 10\ \mu$bar
(see discussion in Huang et al. 2017 for details). This suggests the presence
of a thick layer of atomic H which can have a line center optical depth of
$\tau_{0}\ {\sim}\ 10^{8}$. A careful treatment of resonant scattering is
necessary in order to construct accurate models of H atom excitation, heating,
and radiative forces.

Due to the technical challenge of including resonant scattering, the fully
three-dimensional geometry, and the presence of an outflow, numerical
simulations may be required to fully understand the dynamics of these
irradiated exoplanet atmospheres. The large optical depths at Ly$\alpha$ line
center impose a steep computational cost for solving radiative transfer with
Monte Carlo methods directly coupled with fluid dynamics (Smith et al., 2017).
The number of scatterings a photon undergoes is proportional to the line
center optical depth, $\tau_{0}$, of the domain (Adams, 1972). Near the base
of the atomic layer, the line center optical depth is ${\sim}10^{8}$, so most
of the time is spent following photons in these cells. A method that can
accurately characterize transfer through these zones without following every
photon scattering has the potential to greatly accelerate the calculation (see
e.g. Auer 1968; Ahn et al. 2002).

Approximate analytic solutions for resonant scattering exist in certain
limits. Harrington (1973) showed that when most of the radiation is in the
damping wings, the transfer equation reduces to the Poisson equation. However,
their solution uses an ansatz to handle the boundary condition. To our
knowledge, the errors introduced by this treatment have never been quantified.
They attempt a separation of variables as
$J(\tau,\sigma)=\theta(\tau)j(\sigma)$ in spatial variable $\tau$ and
frequency variable $\sigma$ (their Equations 16 and 23). The solutions for the
eigenfunctions $\theta(\tau)$ and $j(\sigma)$ then depend explicitly on the
separation constant $\lambda$. In order to satisfy the boundary conditions,
the separation constant is shown to satisfy an eigenvalue equation of the form

$\displaystyle\lambda\tan(\lambda B)$ $\displaystyle=$
$\displaystyle\frac{3}{2}\phi\Delta,$ (1)

where $2B$ is the slab optical depth at line center, $\phi$ is the line
profile, and $\Delta$ is the Doppler width. The key point is that the line
profile depends on one of the coordinates: frequency. This causes the
eigenvalues of the separation constant to be frequency-dependent. Thus, the
separation “constant” is not constant, and the function does not satisfy the
Poisson equation since the frequency derivatives will act on the separation
“constant”, giving extra terms. In the limit of large optical depth $B$, they
approximate the eigenvalues as $\lambda_{n}B\simeq\pi(n-1/2)$, which gives
zero mean intensity at the surface. Their Equation 34 subsequently allows
$\lambda$ to have a small deviation from the above expression, which is
explicitly frequency dependent. This allows a nonzero intensity at the
surface, but at the cost of rendering the separation of variables assumption
invalid. Our treatment, using the correct boundary condition, quantifies the
errors in this ansatz.

Several other works have followed Harrington (1973). Neufeld (1990) extends
the solution to media of intermediate optical depth, including the effects of
scattering in the Doppler core of the line. Dijkstra et al. (2006) generalize
the same problem to spherical geometry, as is used here. Lao & Smith (2020)
generalize both the slab and sphere solutions to arbitrary power-law density
and emissivity profiles. Each of these works, and several others (Seon & Kim,
2020; Tomaselli & Ferrara, 2021), use either the same surface boundary
condition and ansatz as Harrington (1973), or use a solution that does not
handle the frequency-dependence of the boundary condition. Our novel steady-
state solution involves a frequency-dependent correction to the solution that
fixes an observed excess at the spectral peaks as compared with Monte Carlo,
which is present in many of the works cited above.

The motivation for including time-dependence in the transfer equation is to
characterize the distribution of photon escape times, which is needed to
calculate the radiation moments in the Monte Carlo simulation. Additionally,
steady-state solutions to this problem are not always sufficient to describe
all the physics of Ly$\alpha$ transport. Time-variable, optically-thick
environments necessitate a time-dependent solution to include the dynamic
effects of Ly$\alpha$ transfer. These include the optical afterglow of gamma-
ray bursts (Roy et al., 2010) and Ly$\alpha$ sources redshifted by
cosmological expansion (Xu et al., 2011), among others.

##  2 STEADY-STATE SOLUTION

Consider a sphere of radius $R$ with uniform density $n_{\rm sc}$, luminosity
$L$, and line-center optical depth $\tau_{0}$, containing a point source of
photons. We aim to find the intensity within the sphere as a function of
radius and photon frequency. The point source is assumed to be a delta
function in space and photon frequency. Photons of frequency $\nu$ near the
line center frequency $\nu_{0}$ are considered. The photon frequency of the
source is $\nu_{s}$. The Doppler width is $\Delta=\nu_{0}v_{\rm th}/c$, where
$v_{\rm th}=\sqrt{2k_{\rm B}T/m_{\rm H}}$ is the thermal speed of hydrogen
atoms of mass $m_{\rm H}$ and temperature $T$, and $c$ is the speed of light.
The photon frequency in Doppler units is $x=(\nu-\nu_{0})/\Delta$, and
$x_{s}=(\nu_{s}-\nu_{0})/\Delta$ is the corresponding source frequency. For
upper-state de-excitation rate $\Gamma$, the ratio of natural to Doppler
broadening is $a=\Gamma/(4\pi\Delta)$. For the Ly$\alpha$ transition and
T=$10^{4}$ K, $a=4.72\times 10^{-4}$. $\mathcal{H}(x,a)$ is the Voigt
function, and the Voigt line profile is
$\phi=\mathcal{H}(x,a)/(\sqrt{\pi}\Delta)$, which is normalized as $\int
d\nu\,\phi(\nu)=1$. The line center optical depth is
$\tau_{0}=kR/(\sqrt{\pi}\Delta)$, where $k=n_{\rm sc}\pi e^{2}f/(m_{e}c)$.
Here, $e$ and $m_{e}$ are the charge and mass of the electron, and $f$ is the
oscillator strength of the transition, which is 0.4162 for Ly$\alpha$ (Rybicki
& Lightman, 1986).

Appendix A contains a derivation of the transfer equation for convenience.
Starting with the full transfer equation, Equation (A24), ignoring photon
destruction and including a photon emission term given by Equation (A25), the
steady-state transfer equation is

$\displaystyle\nabla^{2}J+\left(\frac{k}{\Delta}\right)^{2}\frac{\partial^{2}J}{\partial\sigma^{2}}$
$\displaystyle=$
$\displaystyle-\frac{\sqrt{6}kL}{4\pi\Delta^{2}}\delta^{3}(\mathbf{x}-\mathbf{x}_{s})\delta(\sigma-\sigma_{\rm
s}).$ (2)

where $J$ is the mean intensity, the spatial variable is $\mathbf{x}$, and
$\mathbf{x}_{s}$ is the position of the source. We will consider only the case
where $\mathbf{x}_{s}=0$. Following Harrington (1973), we have used a change
of variables in photon frequency from $x$ to $\sigma$,

$\displaystyle\sigma(x)=\sqrt{\frac{2}{3}}\int_{0}^{x}\frac{dx}{\phi(x)\Delta}\approx\sqrt{\frac{2}{3}}\frac{\pi}{a}\frac{x^{3}}{3},$
(3)

where the approximation is applicable in the damping wing. From Equation
(A17), the line profile is then approximately

$\displaystyle\phi\approx\frac{a}{\pi
x^{2}\Delta}\approx\frac{1}{3\Delta}\left(\frac{2a}{\pi}\right)^{1/3}|\sigma|^{-2/3}.$
(4)

In Equation (2), $\sigma_{s}\equiv\sigma(x_{s})$ is the photon frequency of
the source. $\sigma_{s}$ is interchangeable with $x_{s}$ and $\nu_{s}$ in
Doppler widths or Hz, respectively. Balancing the two terms on the left-hand
side of Equation (2) gives $\sigma\ {\sim}\ \tau_{0}$, or $x_{\rm peak}\
{\sim}\ (a\tau_{0})^{1/3}$. The boundary condition of no incoming intensity at
the surface (Rybicki & Lightman, 1986) is

$\displaystyle J$ $\displaystyle=$ $\displaystyle\sqrt{3}H$ (5)

at $r=R$.

A solution for the mean intensity $J_{d}$ which is divergent at the origin and
$\sigma=\sigma_{s}$ and is zero at infinity is presented in Neufeld (1990).
Here it is extended to spherical geometry and generalized to allow emission
frequencies away from line center:

$\displaystyle J_{\rm d}$ $\displaystyle=$
$\displaystyle\left(\frac{\sqrt{6}k^{2}L}{16\pi^{3}\Delta^{3}}\right)\left(\frac{1}{(kr/\Delta)^{2}+(\sigma-\sigma_{\rm
s})^{2}}\right)$ (6)

$\displaystyle H_{\rm d}$ $\displaystyle=$
$\displaystyle-\frac{1}{3k\phi}\frac{\partial J_{d}}{\partial
r}=\left(\frac{1}{3k\phi}\right)\left(\frac{\sqrt{6}k^{3}L}{8\pi^{3}\Delta^{4}}\right)\left(\frac{kr/\Delta}{\left[(kr/\Delta)^{2}+(\sigma-\sigma_{\rm
s})^{2}\right]^{2}}\right).$ (7)

This solution is useful as a simple analytic formula. However, it is not a
good approximation to the true solution, as it is too large at $r=R$ by a
factor of $J_{\rm d}(R,\sigma)/H_{\rm d}(R,\sigma)\sim
a\tau_{0}/x^{2}\sim(a\tau_{0})^{1/3}\gg 1$ and does not adhere to the correct
boundary condition. This solution is included in Figure 1 for illustration.

A better approximation to the true solution has been derived by Dijkstra et
al. (2006), who generalized the closed-form solution in slab geometry found in
Harrington (1973). It satisfies a $J=0$ boundary condition at $r=R$. Again, we
generalize their solution to allow emission at frequency $\sigma_{\rm s}$ away
from line center. The result can be written as a sum over spatial modes,

$\displaystyle
J_{0}=\frac{\sqrt{6}L}{16\pi\Delta}\frac{1}{R^{2}}\sum_{n=1}^{\infty}n\frac{\sin{\kappa_{n}r}}{\kappa_{n}r}\exp{\left(\frac{-\kappa_{n}\Delta}{k}|\sigma-\sigma_{s}|\right)},$
(8)

and

$\displaystyle H_{0}=-\frac{1}{3k\phi}\frac{\partial J_{0}}{\partial
r}=-\frac{1}{3k\phi}\frac{\sqrt{6}L}{16\pi\Delta}\frac{1}{R^{2}}\sum_{n=1}^{\infty}n\left(\frac{\cos{\kappa_{n}r}}{r}-\frac{\sin{\kappa_{n}r}}{\kappa_{n}r^{2}}\right)\exp{\left(\frac{-\kappa_{n}\Delta}{k}|\sigma-\sigma_{s}|\right)},$
(9)

where $\kappa_{n}=n\pi/R$. These can be summed to give the closed form
expressions

$\displaystyle J_{0}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{6}L}{32\pi^{2}\Delta}\frac{1}{Rr}\left(\frac{\sin(\pi
r/R)}{\cosh\left[\frac{\pi\Delta}{kR}(\sigma-\sigma_{s})\right]-\cos(\pi
r/R)}\right)$ (10)

and

$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{3k\phi}\frac{\sqrt{6}L}{32\pi^{2}\Delta}\frac{1}{Rr^{2}}\left(\frac{\sin(\pi
r/R)}{\cosh\left[\frac{\pi\Delta}{kR}(\sigma-\sigma_{s})\right]-\cos(\pi
r/R)}\right.$ (11) $\displaystyle\left.-\left(\frac{\pi
r}{R}\right)\frac{\cos(\pi
r/R)}{\cosh\left[\frac{\pi\Delta}{kR}(\sigma-\sigma_{s})\right]-\cos(\pi
r/R)}+\left(\frac{\pi r}{R}\right)\frac{\sin^{2}(\pi
r/R)}{\left[\cosh\left[\frac{\pi\Delta}{kR}(\sigma-\sigma_{s})\right]-\cos(\pi
r/R)\right]^{2}}\right).$

These solutions agree with Equations (6) and (7) when the arguments of the
trigonometric and hyperbolic functions are small. Again $J_{0}\gg H_{0}$,
except near $r=R$, where it goes to zero. The flux at $r=R$ can be written

$\displaystyle H_{0}(R,\sigma)$ $\displaystyle=$
$\displaystyle-\frac{1}{3k\phi}\frac{\sqrt{6}L}{16\pi\Delta}\frac{1}{R^{3}}\sum_{n=1}^{\infty}n(-1)^{n}\exp{\left(\frac{-\kappa_{n}\Delta}{k}|\sigma-\sigma_{s}|\right)}$
(12) $\displaystyle=$
$\displaystyle\frac{1}{3k\phi}\frac{\sqrt{6}L}{32\pi\Delta}\frac{1}{R^{3}}\left(\frac{1}{\cosh\left[\frac{\pi\Delta}{kR}(\sigma-\sigma_{s})\right]+1}\right).$

Equation (11) will be shown to be a better approximation to the solution than
Equation (7). It is still valid near the delta function at $r=0$, but is also
a better approximation at $r=R$. $J_{0}$ decreases exponentially, rather than
as a power-law in frequency as it does for $J_{d}$, giving a much smaller flux
in the line wings as compared to the divergent solution.

In order to enforce the boundary conditions, a different solution method is
attempted here, namely a continuous Fourier expansion in the frequency
variable $\sigma$. The solution of this problem is split into two pieces:
$J_{0}$ which includes the delta function source and satisfies $J=0$ at $r=R$,
and $J_{\rm bc}$ which allows the boundary condition $J=\sqrt{3}H$ to be
satisfied at $r=R$. The total solution is

$\displaystyle J(r,\sigma)$ $\displaystyle=$ $\displaystyle
J_{0}(r,\sigma)+J_{\rm bc}(r,\sigma)$ (13)

and

$\displaystyle H(r,\sigma)$ $\displaystyle=$ $\displaystyle
H_{0}(r,\sigma)+H_{\rm bc}(r,\sigma).$ (14)

The additional term $J_{\rm bc}$ must then be a solution of the homogeneous
equation

$\displaystyle\frac{\partial^{2}J_{\rm bc}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial J_{\rm bc}}{\partial
r}+\left(\frac{k}{\Delta}\right)^{2}\frac{\partial^{2}J_{\rm
bc}}{\partial\sigma^{2}}$ $\displaystyle=$ $\displaystyle 0$ (15)

with no delta function source term, and it must allow the boundary conditions
to be satisfied at the surface. Since $J_{0}(R,\sigma)=0$, the surface
boundary condition becomes

$\displaystyle J_{\rm bc}(R,\sigma)-\sqrt{3}H_{\rm bc}(R,\sigma)$
$\displaystyle=\sqrt{3}H_{0}(R,\sigma).$ (16)

Inserting a frequency dependence $J_{\rm bc}\propto e^{is\sigma}$, for
“wavenumber” $s$, gives the equation for modified spherical Bessel functions
of the first kind, $i_{0}(z)=\sinh(z)/z$ for the radial dependence. The
solution can then be represented as

$\displaystyle J_{\rm bc}(r,\sigma)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\frac{ds}{2\pi}e^{is\sigma}A(s)\frac{i_{0}(krs/\Delta)}{i_{0}(kRs/\Delta)},$
(17)

where $A(s)$ is the Fourier amplitude. Inserting Equation (17) into Equation
(16) leads to the following equation for the Fourier amplitudes,

$\displaystyle\int_{-\infty}^{\infty}\frac{ds}{2\pi}e^{is\sigma}A(s)\left[1+\left(\frac{s}{\sqrt{3}\Delta\phi}\right)\left(\frac{i_{0}^{\prime}(kRs/\Delta)}{i_{0}(kRs/\Delta)}\right)\right]$
$\displaystyle=$ $\displaystyle\sqrt{3}H_{0}(R,\sigma).$ (18)

Discretization of Equation (18) for frequency variables $\sigma_{i}$ and
wavenumbers $s_{j}$ leads to a set of coupled linear equations for the
$A(s_{j})$. We use equally-spaced points $\delta\sigma=2\sigma_{\rm
max}/(N-1)$ and $\delta s=2\pi/(N\delta\sigma)$, where $N$ is the number of
points for each grid. The maximum frequency is set as $\sigma_{\rm max}={\rm
constant}\times\tau_{0}$, for a large enough constant that the end of the
frequency grid is at such small intensities that it does not affect the
solution except close to the boundaries. The number of points was increased
until the solution was well-resolved near line center, and only became
inaccurate close to the boundaries. We found that values of $N=4097$ and
$\sigma_{\rm max}=60\tau_{0}$ were sufficient. Given the Fourier amplitudes
$A(s)$, $J_{\rm bc}$ is computed using Equation (17), and the flux is given by

$\displaystyle H_{\rm bc}(r,\sigma)$ $\displaystyle=$
$\displaystyle-\frac{1}{3k\phi}\frac{\partial J_{\rm bc}(r,\sigma)}{\partial
r}=-\frac{1}{3k\phi}\int_{-\infty}^{\infty}\frac{ds}{2\pi}e^{is\sigma}A(s)\left(\frac{ks}{\Delta}\right)\left(\frac{i_{0}^{\prime}(krs/\Delta)}{i_{0}(kRs/\Delta)}\right).$
(19)

The Bessel functions are finite at the center and rise steeply toward the
surface when $kRs/\Delta\gg 1$.

###  2.1 Scaling with Line Center Optical Depth $\tau_{0}$

We now estimate the scaling of $H_{\rm bc}$ with $\tau_{0}$. In the limit
$J_{\rm bc}\gg H_{\rm bc}$, we find that $J_{\rm bc}\approx\sqrt{3}H_{0}$ from
Equation (16). We estimate $H_{\rm bc}$ from $J_{\rm bc}$ using Equation (19)
as

$\displaystyle H_{\rm
bc}(R,\sigma)\approx\frac{1}{\sqrt{3}k\phi}\frac{ks}{\Delta}H_{0}\ {\sim}\
H_{0}s\frac{x^{2}}{a}\ {\sim}\
H_{0}\frac{1}{\tau_{0}}\frac{(a\tau_{0})^{2/3}}{a}\ {\sim}\
H_{0}(a\tau_{0})^{-1/3},$ (20)

where we have used $s\ {\sim}\ 1/\sigma\ {\sim}\ 1/\tau_{0}$ so that

$\displaystyle\frac{H_{\rm
bc}(R,\sigma)}{H_{0}(R,\sigma)}\propto(a\tau_{0})^{-1/3}.$ (21)

At large $\tau_{0}$, it is expected that the correction term will be small,
but it will become increasingly important as $\tau_{0}$ decreases. Our
solution of the transfer equation is only valid when the peaks of the spectral
energy distribution lie well outside of the Doppler core, i.e., for large
$\tau_{0}$. The value of $x$ at which the Doppler and Lorentzian components of
the line profile are equal is $x_{\rm cw}=3.3$. Setting $x_{\rm cw}=x_{\rm
peak}$ and solving for $\tau_{0}$ gives the value at which the peak of the
spectrum falls at the Doppler core boundary, which is $\tau_{\rm cp}\approx
10^{5}$ (“core-peak” optical depth). Hence $H_{\rm bc}/H_{0}{\sim}(\tau_{\rm
cp}/\tau_{0})^{1/3}$ is large at $\tau_{0}\leq\tau_{\rm cp}$ and decreases
relatively slowly as $\tau_{0}$ increases. Additionally, the optical depth at
$x_{\rm peak}$ is proportional to $(a\tau_{0})^{1/3}$, so photons here become
optically thin when $a\tau_{0}{\sim}1$.

Figure 1: Spectrum $P(x)$ vs. frequency $x$ for $\tau_{0}=10^{7}$ with
$x_{s}=0$. At $T=10^{4}$ K, this corresponds to $a\tau_{0}=4.7\times 10^{3}$.
The legend to the right describes each line style. For reference, $\rm H_{0}$:
fiducial solution, $\rm H_{d}$: divergent solution, $\rm H_{bc}$: our boundary
condition correction to the fiducial solution. The top panel is linear scale,
the middle panel is log scale with $|\rm H_{bc}|$ shown instead of $\rm
H_{bc}$, and the bottom panel is the residual of each solution with Monte
Carlo.

###  2.2 Comparison to Monte Carlo

The Monte Carlo method is used to solve the transfer equation numerically in
order to compare the analytic approximation to an “exact” solution. This
method is valid at all $\tau_{0}$, being restricted only by the computational
demand, which grows proportionally to the number of photons used and
$\tau_{0}$. For each simulation, a total of ${\sim}10^{6}$ photon packets are
initialized at a monochromatic source frequency $x_{s}$ and are allowed to
propagate through the sphere until escaping, at which point their positions,
outgoing angles, and escape frequencies are tabulated to obtain the spectrum
at the surface of the spherical simulation domain. A constant temperature of
$T=10^{4}\ \rm K$ is set for the gas. Frequency redistribution is calculated
at each scattering. In the comparisons shown in this section, the raw photon
data is binned in frequency to obtain spectra. Further details of the Monte
Carlo implementation are discussed in Huang et al. (2017).

We now compare each of the previously-discussed solutions for surface flux to
the Monte Carlo results. The spectrum $P(x)$ is defined as the specific
luminosity at the surface divided by the source luminosity, or

$\displaystyle P(x)=\frac{16\pi^{2}R^{2}H(R,x)\Delta}{L}.$ (22)

This is normalized so that $\int P(x)dx=1$. Since $H(R,x)$ is per $d\nu$, a
factor of $\Delta$ gives the expression the correct units.

In Figure 1, the Monte Carlo spectrum is shown along with that of the
solutions $H_{\rm d}$, $H_{0}$, and $H_{0}+H_{\rm bc}$ for an optical depth of
$\tau_{0}=10^{7}$ and photons emitted at line center $\rm x_{s}=0$. Note that
the errorbars shown on the Monte Carlo data points are proportional to
$\sqrt{N}$, with $N$ being the photon count in each frequency bin, since the
photons are all equally weighted. The $H_{\rm bc}$ term is negative at the
peak of the spectrum and positive in the line wing such that, when added to
$H_{0}$, it corrects for the apparent excess of flux in the peaks of the
spectrum. The solution with the correct frequency-dependent boundary condition
enforced, $H_{0}+H_{\rm bc}$, has lower residuals to Monte Carlo results than
the other solutions, especially in the line wing. The boundary term corrects
the deficit of $H_{0}$ in the line wings, further improving agreement with the
numerical result. The residuals to the $H_{0}$ solution are a close match to
the $H_{\rm bc}$ term, since the Monte Carlo represents the “true” solution,
$H$, and $H_{\rm bc}=H-H_{0}$. It is evident that the divergent solution
$H_{\rm d}$ fails in the line wings. Also note that the “V” shape of the
solution in the line core is due to the low number of points plotted, as the
analytic solutions are not valid in this frequency regime since they utilize
the damping wing approximation of the Voigt line profile.

Figure 2: The same as Figure 1, but for $\tau_{0}=10^{5}$ (top panel),
$10^{6}$ (middle panel), and $10^{7}$ (bottom panel). The $x$ and $y$ axes are
scaled by $(a\tau_{0})^{1/3}$.

The size of $H_{\rm bc}$ is dependent on $\tau_{0}$. $H_{\rm bc}$ is
significant even at $\tau_{0}{\sim}10^{7}$ where the $H_{0}$ solution is
expected to perform well, i.e., photons are pushed further out into the wing
where the simplifying assumptions made in the derivation of the differential
equation are a better approximation.

In Figure 2 we show the solutions alongside Monte Carlo, now for three
different optical depths $\tau_{0}=10^{5},10^{6}$, and $10^{7}$. From Equation
(21), the size of the term $H_{\rm bc}$ should become smaller with larger
optical depths, following a $(a\tau_{0})^{-1/3}$ scaling. Indeed, agreement
between Equation (12) and the Monte Carlo points in Figure 2 improves as
$\tau_{0}$ increases, with $H_{\rm bc}$ providing a fractionally smaller
correction to $H_{0}$. One factor of $(a\tau_{0})^{1/3}$ has been scaled out
of the x-axis such that the peaks of the distributions are horizontally
aligned. This scaling has also been applied to the y-axis to preserve
normalization of the escape probability. At lower $\tau_{0}$, the scattering
of photons within the Doppler core of the line becomes important, but our
analytic solution does not include this effect. The effects of line core
scattering can be seen in the Monte Carlo data for $\tau_{0}=10^{5}$ and, to a
lesser extent, $\tau_{0}=10^{6}$.

Figure 3: The same as Figure 1, but for $\tau_{0}=10^{7}$ and $\rm x_{s}=0$
(top panel), $6$ (middle panel), and $12$ (bottom panel). The optical depth at
each of these source frequencies is $\tau_{s}=10^{7},77,$ and $19$,
respectively. Figure 4: The same as Figure 3, but at a lower optical depth
$\tau_{0}=10^{6}$. The shift $x_{s}$ is a much larger fraction of the distance
to the spectral peak $(a\tau_{0})^{1/3}$, and thus the asymmetry in the
spectrum is much larger. The optical depth at the source frequency is
$\tau_{s}=10^{6}$, $7.7$, and $1.9$ for $x_{s}=0,6,$ and $12$, respectively.

Next, we show $P(x)$ for $x_{s}\neq 0$. Photons initialized further out in the
line wing have larger mean free paths. The larger spatial diffusion implies
greater escape probability for these photons. In the limit that $|\rm x_{s}|$
becomes large, the distribution becomes a delta function at $x_{s}$ as all
photons escape the sphere without scattering.

Figure 3 shows calculations performed for $x_{s}=0,6$, and $12$ and
$\tau_{0}=10^{7}$. The asymmetry of the spectrum is slight for $\rm x_{s}=6$
where $\tau(x_{s})=77$, but is larger for $\rm x_{s}=12$ outside the line core
where $\tau(x_{s})=19$. It is seen here that the difference between the Monte
Carlo data and $H_{0}$ becomes larger as $\rm x_{s}$ increases. Thus, for
large $|x_{s}|$, inclusion of $H_{\rm bc}$ is more important.

Figure 4 shows emission away from line center at the same values of $x_{s}$ as
in Figure 3, but for $\tau_{0}=10^{6}$ rather than $10^{7}$. The difference
between the left and right side of the escaping spectrum is now substantial
since $(a\tau_{0})^{-1/3}$ increased by a factor of ${\sim}2$. It is clear
from the figure that as $x_{s}$ extends further into the wing, the spectrum
becomes more strongly peaked in frequency. Additionally, since the sphere is
increasingly optically thin in the wing, we expect there to be a stronger
disagreement with the Monte Carlo as the analytic solution assumed large
optical depths.

Figure 5: $P_{\rm bc}(x)=16\pi^{2}R^{2}H_{\rm bc}(R,x)\Delta/L$ (Equation (22)
with $H\rightarrow H_{\rm bc}$) vs. frequency $x$. Each panel shows
$\sigma_{s}$ further from line center, labeled by
$\sigma_{s}=0,\tau_{0},2\tau_{0}$. The solution is shown at a range of line
center optical depths between $\tau_{0}=10^{5}$ and $10^{9}$ with solid lines.
The location of the core-wing boundary $x_{cw}$ for each $\tau_{0}$ is shown
with dashed vertical lines which match the color of the corresponding
solution. The axes are scaled to show that the size of the correction factor
agrees with the predicted $(a\tau_{0})^{-1/3}$ scaling.

Figure 5 shows how the correction $P_{\rm bc}(x)$, Equation (22) with
$H\rightarrow H_{\rm bc}$, scales with both $\tau_{0}$ and $x_{s}$. In this
figure, $\sigma_{s}$ is shifted by integer multiples of $\tau_{0}$ (Equation
A23) in each panel such that the source falls near the peak of the spectrum
for each $\tau_{0}$. For clarity, only the $x>0$ side of the spectrum is
shown. From Equation (21), it is expected that the fractional size of $H_{\rm
bc}$ relative to $H_{0}$ should become smaller with larger optical depths
following $(a\tau_{0})^{-1/3}$. This factor has been scaled out of the figure
such that solutions for different $\tau_{0}$ and the same $\sigma_{s}$ should
show close agreement in scale on the figure’s vertical axis if the relation
holds. Indeed, the scaled solutions converge as $\tau_{0}$ becomes larger,
indicating agreement with the $(a\tau_{0})^{-1/3}$ scaling. The remaining
discrepancy present in the vertical axis for fixed $\sigma_{s}$ results from
$x_{\rm peak}$ becoming close to $x_{\rm cw}$; at lower $\tau_{0}$, this
causes the line profile approximation in the wing, Equation (A17), to break
down. From this, we conclude that the errors introduced by the incorrect
separation of variables in Harrington (1973), Neufeld (1990), Dijkstra et al.
(2006) and others are indeed proportional to $(a\tau_{0})^{-1/3}$.

##  3 Time-Dependent Diffusion

In order to understand how long it takes for the photons to escape the uniform
sphere of gas we must reintroduce the time dependence of the diffusion
equation, which was ignored in the steady-state calculations in Section 2. To
obtain the radiative intensity $I=dE/(dAdtd\Omega d\nu)$ on timescales
comparable to the light-crossing time $t_{\rm lc}=R/c$, the time-dependent
response to a delta function impulse is found. This allows the distribution of
photon escape times (the “wait time distribution”) to be characterized. For
simplicity, a $J=0$ boundary condition will be used in the following
derivations, which is a rough approximation for $a\tau_{0}\gg 1$.

###  3.1 Derivation of the time-dependent solution

The emissivity for an impulsive source with energy E, source position
$\mathbf{x}_{s}$, and frequency $\nu_{s}$ is derived in Appendix A.
Considering a photon source at $\mathbf{x}_{s}=0$, we have (Equation A25)

$\displaystyle j_{\rm em}$ $\displaystyle=$
$\displaystyle\frac{E}{4\pi}\delta^{3}(\mathbf{x})\delta(\nu-\nu_{\rm
s})\delta(t),$ (23)

The resulting equation for $J(r,\sigma,t)$ is

$\displaystyle\frac{-3k\phi}{c}\frac{\partial J}{\partial
t}+\nabla^{2}J+\left(\frac{k}{\Delta}\right)^{2}\frac{\partial^{2}J}{\partial\sigma^{2}}$
$\displaystyle=$
$\displaystyle-\frac{\sqrt{6}kE}{4\pi\Delta^{2}}\delta^{3}(\mathbf{x})\delta(\sigma-\sigma_{s})\delta(t).$
(24)

We employ an expansion in terms of spherical Bessel functions in $r$ and
Fourier transform in time. The zeroth spherical Bessel function is
$j_{0}(x)=\sin{x}/x$. The expansion for $J(r,\sigma,t)$ is then

$\displaystyle
J(r,\sigma,t)=\sum_{n=1}^{\infty}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}e^{-i\omega
t}j_{0}\left(\kappa_{n}r\right)J(n,\sigma,\omega),$ (25)

with

$\displaystyle J(n,\sigma,\omega)=\frac{2\kappa_{n}^{2}}{R}\int_{0}^{R}dr\
r^{2}j_{0}(\kappa_{n}r)\int_{-\infty}^{\infty}dt\ e^{i\omega t}J(r,\sigma,t).$
(26)

Here, $\kappa_{n}=n\pi/R$ and $\omega$ describes the time-dependence of $J$.
Though it is written as a function of the photon frequency variable $\sigma$,
$J(r,\sigma,t)$ is the specific mean intensity $dE/(dAdtd\nu)$ and is a
distribution in $\nu$. The Fourier coefficient $J(n,\sigma,\omega)$ has units
$dE/(dAd\nu)$. Using Equation (26), we obtain

$\displaystyle\left(\frac{3k\phi}{c}i\omega-\kappa_{n}^{2}\right)J(n,\sigma,\omega)$
$\displaystyle+$
$\displaystyle\left(\frac{k}{\Delta}\right)^{2}\frac{\partial^{2}J(n,\sigma,\omega)}{\partial\sigma^{2}}=-\frac{2\kappa_{n}^{2}}{R}\frac{\sqrt{6}}{4\pi}\frac{kE}{\Delta^{2}}\frac{1}{4\pi}\delta(\sigma-\sigma_{s}).$
(27)

At $\sigma=\sigma_{s}$, continuity must be enforced,

$\displaystyle J(n,\sigma^{-},\omega)=J(n,\sigma^{+},\omega),$ (28)

and the discontinuity in $dJ/d\sigma$ due to the source is

$\displaystyle\frac{\partial
J(n,\sigma^{+},\omega)}{\partial\sigma}-\frac{\partial
J(n,\sigma^{-},\omega)}{\partial\sigma}$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{6}}{8}n^{2}\frac{E}{kR^{3}}.$ (29)

At large values of $\sigma$ the line profile $\phi$ is small and Equation (27)
becomes

$\displaystyle\frac{\partial^{2}J}{\partial\sigma^{2}}\approx\frac{\Delta^{2}\kappa_{n}^{2}}{k^{2}}J,$
(30)

which has solutions

$\displaystyle J(n,\sigma,\omega)\ {\sim}\ e^{\pm\kappa_{n}\sigma\Delta/k}.$
(31)

This approximate solution implies a boundary condition at large $|\sigma|$

$\displaystyle\frac{\partial
J}{\partial\sigma}=\mp\frac{\kappa_{n}\Delta}{k}J,$ (32)

where a negative sign is taken for large $+\sigma$ and a positive sign is
taken for large $-\sigma$ to choose the finite solution as
$|\sigma|\to\infty$. Numerical integrations are performed inward toward
$\sigma_{s}$ over several domains: from large $|\sigma|$ to $\sigma_{s}$, from
large $|\sigma|$ to 0, and from 0 to $\sigma_{s}$, depending on whether
$\sigma_{s}$ is positive or negative. If $\sigma_{s}=0$, just two integrations
are performed inward from large $|\sigma|$ to 0. Initial values for
integration are obtained either by setting $J=1$ and $dJ/d\sigma$ from
Equation (32) at large $|\sigma|$ or by matching $J$ and $dJ/d\sigma$ at 0.
This gives $J$ and $J^{\prime}$ on either side of $\sigma_{s}$, where a prime
indicates the derivative $\partial/\partial\sigma$. By enforcing the matching
conditions, Equations (28) and (29), the eigenfunctions $J(n,\sigma,\omega)$
are obtained over the domain of photon frequencies $\sigma$. Since the
solutions are linear in the starting conditions, only two integrations with
different starting values are necessary.

We now wish to reconstruct the specific mean intensity $J(r,\sigma,t)$. While
one might expect this could be expressed as a sum over eigenmodes, the
analysis presented in Appendix B suggests this treatment is incomplete in the
case where $x_{s}\neq 0$ and the solution is asymmetric about the line center.
This ansatz does, however, roughly agree with Monte Carlo results for
$x_{s}=0$ based on numerical calculations of this result.

Let us define the damping rate to be $\gamma\equiv i\omega$, which is real and
positive for damped solutions. At the eigenvalues $\gamma=\gamma_{nm}$, the
response $J(n,\sigma,\omega)$ is resonant. We find that near these
$\gamma_{nm}$ poles an approximate expression for the resonant response of the
eigenfunctions is

$\displaystyle J(n,\sigma,-i\gamma)$
$\displaystyle\simeq\frac{J_{nm}(\sigma)}{\gamma_{nm}-\gamma}+C(\gamma,\sigma),$
(33)

where $C(\gamma,\sigma)$ varies slowly in $\gamma$. If the $\omega$-integral
in Equation (25) could be closed at infinity and evaluated using the residue
theorem, the result would be

$\displaystyle J(r,\sigma,t)$ $\displaystyle\simeq$ $\displaystyle
j_{0}(\kappa_{n}r)J_{nm}(\sigma)e^{-\gamma_{nm}t}.$ (34)

Summing over all spatial modes $n$ and over all eigenmodes $m$ for a given
$n$, we obtain

$\displaystyle J(r,\sigma,t)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}j_{0}(\kappa_{n}r)\sum_{m=1}^{\infty}J_{nm}(\sigma)e^{-\gamma_{nm}t}.$
(35)

This ansatz captures the contributions from $n\times m$ simple poles. Taking a
derivative with respect to $r$ and evaluating at the surface $r=R$, we use

$\displaystyle\frac{dj_{0}(\kappa_{n}r)}{dr}\bigg{\rvert}_{R}$
$\displaystyle=$
$\displaystyle\frac{d}{dr}\left[\frac{\sin(\kappa_{n}r)}{\kappa_{n}r}\right]\bigg{\rvert}_{R}=\left(\frac{\cos(\kappa_{n}R)}{R}-\frac{\sin(\kappa_{n}R)}{\kappa_{n}R^{2}}\right)=\frac{(-1)^{n}}{R}$
(36)

to obtain the flux, which is

$\displaystyle F(R,\sigma,t)$ $\displaystyle=$
$\displaystyle-\frac{4\pi}{3k\phi}\frac{dJ(R,\sigma,t)}{dr}=-\frac{4\pi}{3k\phi
R}\sum_{nm}(-1)^{n}J_{nm}(\sigma)e^{-\gamma_{nm}t}.$ (37)

Multiplying by $4\pi R^{2}$ gives the energy per time per frequency emerging
from the sphere to be

$\displaystyle\frac{dE}{dtd\nu}$ $\displaystyle=$
$\displaystyle-\frac{16\pi^{2}R}{3k\phi}\sum_{nm}(-1)^{n}J_{nm}(\sigma)e^{-\gamma_{nm}t}.$
(38)

Integrating over time yields a factor $1/\gamma_{nm}$, and by integrating over
$d\nu$ we find

$\displaystyle E$ $\displaystyle=$
$\displaystyle\sqrt{\frac{3}{2}}\frac{16\pi^{2}R\Delta^{2}}{3k}\sum_{nm}(-1)^{n+1}\gamma_{nm}^{-1}\int
d\sigma J_{nm}(\sigma).$ (39)

This non-trivial “sum rule” provides a check on the values of $\gamma_{nm}$
and $J_{nm}(\sigma)$. This expression can also be written as

$\displaystyle 1$ $\displaystyle=$ $\displaystyle\sum_{nm}P_{nm},$ (40)

where the contribution of each mode is

$\displaystyle P_{nm}$ $\displaystyle\equiv$
$\displaystyle\sqrt{\frac{3}{2}}\frac{16\pi^{2}R\Delta^{2}}{3kE}(-1)^{n+1}\gamma_{nm}^{-1}\int
d\sigma J_{nm}(\sigma).$ (41)

These coefficients $P_{nm}$ are negative for odd values of $n$ and positive
for even $n$. The size of each contribution scales roughly as
$0.5/(m-7/8)^{2/3}$, with a weak dependence on $n$. This indicates the need
for a large number of $n$ and $m$ to converge, in that it takes roughly ten
times as many $m$ modes for a given $n$ to reduce the size of $P_{nm}$ by a
factor of ${\sim}$5\. The physical intuition for the convergence of these
terms is that the $n$ spatial terms must provide sufficient spatial resolution
to resolve the steep falloff in intensity at the surface of the sphere.
Additionally, the function falls off steeply in frequency in the line wing,
which requires more $m$ terms in the Fourier sum to resolve (also see the
discussion of Figure 9 in Section 3.3).

###  3.2 Numerical calculation

We seek now to calculate the eigenmodes $J_{nm}(\sigma)$ and eigenfrequencies
$\gamma_{nm}$ for a given spatial $n$. These will be labelled by an index
$m=1,2,...$. To measure the size of the response to detect where resonances
occur, we sum the absolute value of $J(n,\sigma,-i\gamma)$ over the array
$\sigma$. We call this response $f$, and use the index $j$ to represent the
value of the response at discrete points $\gamma_{j}$ over a range of
$\gamma$. In places where $f_{j}>f_{j-1}$ and $f_{j}>f_{j+1}$, we have
bracketed a resonance that occurs in the interval
$(\gamma_{j-1},\gamma_{j+1})$. To refine the value of the eigenfrequency
before continuing the sweep in $\gamma$, we evaluate $f_{j-1}$, $f_{j}$, and
$f_{j+1}$ at the points $(\gamma_{j-1},\gamma_{j},\gamma_{j+1})$. Assuming the
form in Equation (33), a guess at the correct eigenvalue $\gamma_{nm}$ can be
calculated by linear interpolation from

$\displaystyle\gamma_{\rm guess}$ $\displaystyle=$
$\displaystyle\frac{b\gamma_{j-1}-\gamma_{j+1}}{b-1},$ (42)

where

$\displaystyle b$ $\displaystyle=$
$\displaystyle\left(\frac{f_{j}-f_{j-1}}{f_{j}-f_{j+1}}\right)\left(\frac{\gamma_{j+1}-\gamma_{j}}{\gamma_{j-1}-\gamma_{j}}\right).$
(43)

The error of the current guess is $|\gamma_{\rm guess}-\gamma_{j}|$. This
error is reduced iteratively by replacing initial points
$(\gamma_{j-1},\gamma_{j},\gamma_{j+1})$ with closer estimates while the size
of the response grows as the resonance is approached. After iterating an
eigenvalue $\gamma_{nm}$ to convergence, we now find the corresponding
eigenfunction $J_{nm}(\sigma)$. We evaluate Equation 33 at two points
$\gamma_{1}$ and $\gamma_{2}$ near the resonance, subtracting them and solving
for $J_{nm}(\sigma)$ to find

$\displaystyle J_{nm}(\sigma)$ $\displaystyle\simeq$
$\displaystyle\frac{J(n,\sigma,-i\gamma_{1})-J(n,\sigma,-i\gamma_{2})}{1/(\gamma_{nm}-\gamma_{1})-1/(\gamma_{nm}-\gamma_{2})},$
(44)

where $C(\gamma,\sigma)$ has cancelled in the difference.

Figure 6: Eigenfunctions $J_{nm}(x)$ for the lowest-order spatial eigenmode
$n=1$, and $m=1,...,100$ with $x_{s}=0$ and $\tau_{0}=10^{7}$. The scale of
the $J_{nm}(\sigma)$ are set by the factor $E/(kR^{3})$, which here is
${\sim}10^{-37}$ in ergs/(cm2 Hz).

The form of a single eigenmode $J_{nm}(\sigma)$ is oscillatory out to some
turning point, $\sigma_{\rm tp}$, at which point the function becomes
evanescent. The location of the turning point can be found by ignoring the
delta-function discontinuity at the source frequency $\sigma_{s}$ in Equation
(27) and examining the resulting homogeneous differential equation. We obtain

$\displaystyle\frac{d^{2}J}{d\sigma^{2}}$ $\displaystyle=$
$\displaystyle\left[\left(\frac{\kappa_{n}\Delta}{k}\right)^{2}-\frac{3\phi\gamma\Delta^{2}}{ck}\right]J,$
(45)

where the line profile is approximated as in Equation (A17). When the
coefficient on the right hand side is positive, exponential growth or decaying
evanescent solutions are found. This occurs in the line wings. When the
coefficient on the right hand side is negative, oscillatory solutions are
found (propagation), which occurs near the line core. The boundary between
propagation and evanescence occurs at the turning point, given by

$\displaystyle\sigma_{\rm tp}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2a}{\pi}}\left(\frac{k\gamma}{\kappa_{n}^{2}c\Delta}\right)^{3/2}$
(46)

Thus, to ensure accuracy in each term of Equation (35), the bounds of $\sigma$
must be set sufficiently far outside of $\sigma_{\rm tp}$ such that the
function is small at the edges. The scale of an $e$-folding in
$J_{nm}(\sigma)$ is $k/(\kappa_{n}\Delta)=\tau_{0}/(\sqrt{\pi}n)$, so a grid
of $\sigma$ is chosen that spans a large enough number of $e$-foldings that no
oscillatory behavior is present at the boundaries of the domain.

The eigenfunction’s oscillatory forms have varying amplitudes which sum in
Equation (35) to create the final form of the mean intensity. The largest
contribution at late times always comes from the ($n=1,m=1$) lowest-order
eigenfunction. Figure 6 shows a set of eigenfunctions $J_{nm}(\sigma)$ to
illustrate their relative scales for different $m$ at a fixed spatial
eigenmode $n$. The overall scale of the $J_{nm}(\sigma)$ are set by the factor
$E/(kR^{3})$ with $E$ arbitrarily set to 1. For H atoms with $T=10^{4}$ K and
$\tau_{0}=10^{7}$, an eigenfunction has typical size
${\sim}Ea/\left(R^{2}\Delta\right)=10^{-37}$ in units of specific mean
intensity times time. Additional terms add smaller-magnitude, faster-
oscillating components that lead to higher accuracy upon summation with the
lower-order terms in Equation (35). The oscillations of various modes must
cancel in the Fourier sum, so many modes $m$ and $n$ are required for
convergence to the solution.

Figure 7: Dimensionless decay time vs. mode number. Overlapping solid lines
are plotted for each $n$. The thickness of the lines decrease with $n$ in
order to clearly show where they are in agreement. These numerically-obtained
resonant frequencies are compared with the analytic expression in Equation
(47), shown as a dashed line. The parameters for the calculation are
$\tau_{0}=10^{7}$ and $x_{s}=0$. The value of $n$ is given by the color bar on
the right-hand side.

The values of the $\gamma_{nm}$ can be described approximately with Equation
(47). Their values depend on $m$, $n$, and other physical parameters according
to

$\displaystyle\gamma_{nm}=2^{-1/3}\pi^{13/6}n^{4/3}\left(m-\frac{7}{8}\right)^{2/3}\frac{c}{R}(a\tau_{0})^{-1/3}$
(47)

as shown by setting the denominator of Equation B5 to zero in WKB
approximation of Appendix B. The power law in $m$ is weak, requiring up to
$m=1000$ to reduce the scale of $\gamma_{nm}^{\ \ -1}$ by two orders of
magnitude. When sweeping through to find resonances, Equation (47) is used to
set the scale of the sweep points $\gamma_{j}$ to ensure no $\gamma_{nm}$ are
missed. The close agreement with the analytic expression shown in Figure 7
indicate that the numerical solutions are accurate.

###  3.3 Comparison with Steady State and Monte Carlo

We now calculate the wait time distribution for escape from the sphere. This
is obtained by integrating Equation (38) over all frequencies. We find

$\displaystyle P(t)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{3}{2}}\frac{16\pi^{2}R\Delta^{2}}{3kE}\sum_{nm}(-1)^{n+1}e^{-\gamma_{nm}t}\int
d\sigma J_{nm}(\sigma)$ (48) $\displaystyle=$
$\displaystyle\sum_{nm}P_{nm}\gamma_{nm}e^{-\gamma_{nm}t},$

which is normalized to unity. For a sufficiently large number of spatial modes
$n$ and frequency modes $m$, the result of this sum can agree with Monte Carlo
escape time distributions when $x_{s}=0$. The late-time distribution is simply
an exponential falloff. The rate constant of the falloff is the lowest-order
eigenfrequency, $\gamma_{11}$, and its scale is determined by the coefficient
$P_{11}$ as in Equation (41). Thus, an approximate “fitting function” that
captures both the peak of the escape time distribution and the exponential
falloff is

$\displaystyle P(t)=\exp{\left[-\left(\frac{t_{\rm
diff}}{t}\right)^{2}\right]}\times\gamma_{11}P_{11}e^{-\gamma_{11}t}.$ (49)

The first term represents the early-time distribution, which then transitions
to an exponential falloff past a point $ct_{\rm diff}/R=(a\tau_{0})^{1/3}$,
where $t_{\rm diff}$ is the characteristic diffusion timescale.

In Figure 8, the late-time decay timescale of the wait time distribution is
shown as a function of $\tau_{0}$. It is shown that the time constant of
exponential decay in fitted Monte Carlo escape time distributions converges
with $\gamma_{11}^{-1}$ at sufficiently high $\tau_{0}$, following a
$t\propto(a\tau_{0})^{1/3}$ scaling. The coefficient of this scaling ($0.51$)
is within a factor of 2 of the approximate “light-trapping time” defined in
Lao & Smith (2020), which predicts $ct/R=0.901(a\tau_{0})^{1/3}$. At lower
$\tau_{0}$, the effects of line core scattering are most important, leading to
a larger discrepancy in the characteristic escape timescale. Here, the Monte
Carlo accurately includes the photons which scatter in the core many times
before escaping, while the semi-analytic solution does not capture this
behavior as it uses only the Lorentzian piece of the line profile, and also
does not use enough spatial modes to accurately model the frequency regime
near line center. However, as $\tau_{0}$ grows, the effect of core scattering
becomes smaller and the approximations hold, agreeing better with the expected
$(a\tau_{0})^{1/3}$ scaling (Adams, 1975) for the rate constant at late times.
The excess in the Monte Carlo data points due to core scattering decreases
exponentially at higher $\tau_{0}$, and though these points are not shown at
$\tau_{0}=10^{8}$ and $10^{9}$ due to computational expense, it is expected
that the fractional error between the Monte Carlo and the $(a\tau_{0})^{1/3}$
scaling would be less than 2% at $\tau_{0}=10^{9}$.

Figure 8: Late-time exponential decay timescale as a function of $\tau_{0}$.
The Monte Carlo points on this figure are obtained by fitting only the
exponential piece of the escape time distribution to obtain the rate constant.

We now evaluate the time-integrated spectrum (fluence) of the response to an
impulse and compare it with the solution for the $H_{0}$ steady-state spectrum
(Equation 12). Integrating Equation (38) over all times and dividing by the
energy $E$, we find the fluence

$\displaystyle P(x)$ $\displaystyle=$
$\displaystyle\frac{16\pi^{2}R\Delta}{3k\phi
E}\sum_{nm}(-1)^{n+1}\gamma_{nm}^{-1}J_{nm}(\sigma).$ (50)

Integrating over $\nu$ then gives unity as required by the sum rule in
Equation (40).

In Figure 9, the fluence for $x_{s}=0$ and $\tau_{0}=10^{7}$ is shown for a
sum up to $n=20$ and $m=500$, labelled “Time-integrated”, and is compared with
two analytic solutions: the steady-state $H_{0}$ solution (Equation 12),
labelled “Steady State”, and the result for summing a finite number of spatial
modes in the steady-state eigenfunction expansion as in the first line of
Equation (12), labelled “Partial Sum”. Additional spatial modes $n$ increase
the solutions’ accuracy in the core of the line. If more spatial modes are
included, the agreement with the steady-state spectrum extends further toward
the line core. If additional frequency modes are included, faster-oscillating
terms are incorporated into the Fourier sum over eigenmodes which create more
perfect cancellations with the lower-order terms, reducing the “ringing” seen
in the time-integrated spectrum. Extending the calculation deep into the line
core by adding additional spatial modes could have an impact on the accuracy
of the escape time distribution, but this would primarily affect the
distribution at early times since the late time distribution is determined by
the lowest order modes. This was the motivation for choosing a comparatively
low number of spatial eigenmodes with respect to the number of frequency
eigenmodes calculated.

Figure 9: Fluence $P(x)$ vs. $x$ (see Equation 50). Fluence is the radiation
flux integrated over time. Steady-state and time-integrated spectra for
$n=1,...,20$ and $m=1,...,500$ are shown with $x_{s}=0$ and $\tau_{0}=10^{7}$.
Note that the x-axis begins near the edge of the line core, as we are only
concerned with the solutions’ accuracy near the line wing.

In Figure 10, the escape time distributions calculated from Equation (48) are
shown alongside Monte Carlo and the fitting function Equation (49) for
$\tau_{0}=10^{6},10^{7}$ with $x_{s}=0$. The disagreement between the tail of
the distribution and the Monte Carlo data is due to line core scattering which
is not modeled by the eigenfunction solution, but improves for larger optical
depth as seen in the figure. A large number of scatterings in the Doppler core
affects the tail of the escape-time distribution, since photons with
frequencies near line center will take longer to escape. Thus, the rate
constant for the exponential falloff is overestimated slightly in the
eigenfunction solution as compared with the Monte Carlo. The error in this
rate constant is a function of $\tau_{0}$ since the effect from the Doppler
core is greatest when it extends into the peak of the spectrum.

Figure 10: Wait-time distribution of escaping photons $P(t)$ vs. t for
$\tau_{0}=10^{6}$ and $10^{7}$, including the fitting function from Equation
(49). A sum over 20 spatial eigenmodes and 500 frequency eigenmodes is labeled
“Eigenfunctions”. All calculations were performed with a monochromatic source
of photons at line center ($x_{s}=0$).

##  4 Discussion

###  4.1 Steady-State Source

A primary goal of this work is to present a solution for resonant scattering
of photons near the line-center frequency $\nu_{0}$ in a uniform sphere. We
have generalized a spherically symmetric solution derived by Dijkstra et al.
(2006) (called $H_{0}$ here) to allow a monochromatic source of photons with
frequencies away from line center. We introduce a new term to this solution,
$J_{\rm bc}$, which allows the boundary condition $J=\sqrt{3}H$ to be
satisfied at the surface of the sphere. This is solved using a continuous
Fourier expansion in frequency. The integrals are discretized and the Fourier
coefficients solved for numerically. The resulting flux correction, $H_{\rm
bc}$, scales as $H_{0}(a\tau_{0})^{-1/3}$. Thus, for large $a\tau_{0}$, only a
small correction to $H_{0}$ is needed, while larger errors are present in
calculations performed at lower $a\tau_{0}$. Since the Laplacian form for
frequency redistribution in the differential equation is only correct for
photons in the wing where the line profile is $\phi\approx a/(\pi
x^{2}\Delta)$, our solutions do not accurately model the Doppler core of the
Ly$\alpha$ line. Because the peak of the spectral energy distribution of
escaping photons is $x_{\rm peak}{\sim}(a\tau_{0})^{1/3}$, calculations
performed at small $a\tau_{0}$ are inaccurate due to the close proximity of
the spectral peak and the Doppler core of the line.

By comparison with Monte Carlo simulations, we have shown that the enforcement
of the correct frequency-dependent boundary conditions improves the accuracy
of these analytic solutions for $a\tau_{0}\gg 1$. Specifically, this solution
shows improvement over previous solutions that utilized a $J=0$ surface
boundary condition presented in Harrington (1973), Neufeld (1990), and
Dijkstra et al. (2006). Several papers have previously compared these analytic
models to Monte Carlo and seen discrepancies on the order of this correction.
For example, in the top-left panel of Figure 1 from Dijkstra et al. (2006),
the Ly$\alpha$ spectrum emergent from a sphere of uniform optical depth is
shown for $\tau_{0}=10^{5},10^{6},$ and $10^{7}$ at a temperature of $T=10$ K,
corresponding to $a=1.5\times 10^{-2}$. The dotted line showing their
theoretically-derived spectrum ($H_{0}$) displays an excess at the peak of at
least 5-10 percent as compared with the Monte Carlo for $\tau_{0}=10^{5}$ and
$10^{6}$. Another example is in Smith et al. (2015), where the peak excess in
the Ly$\alpha$ spectrum is particularly noticeable for line center optical
depths of $\tau_{0}=10^{6}$ and $10^{7}$ in the top panel of their Figure 5,
which used slab geometry and a gas temperature of $T=10^{4}$ K. Again, the
error in their solution is of order 5-10 percent. Both of these solutions are
too large at the spectral peaks and too small further out in the wing, and the
error scales approximately as $(a\tau_{0})^{-1/3}$. We show in our Figure 2
that the error present in $H_{0}$ is corrected by our treatment of the
boundary condition at $\tau_{0}=10^{7}$ for T$=10^{4}$ K, corresponding to
$a=4.72\times 10^{-4}$. We note that our correction term $H_{\rm bc}$ is
positive in the line wing and negative at the peak of the spectrum, which
matches with the discrepancies noted in the aforementioned solutions.

###  4.2 Impulsive Source

The time-dependent transfer equation is solved in order to characterize the
distribution of photon escape times. A semi-analytic approach is used,
utilizing an expansion in space, time, and photon frequency. This boundary
value problem in frequency $\sigma$ is solved to find the flux at the surface
of the sphere as a function of $t$ and $\nu$. This solution is expressed as a
sum over spatial and frequency modes $n$ and $m$, respectively. Calculating
additional spatial eigenmodes increases the accuracy nearer to line center,
but convergence is slow due to each eigenmode’s weak dependence on $n$.
Additional frequency eigenmodes introduce fast-oscillating terms that improve
the accuracy of the Fourier sum, as their contributions cancel with components
of lower-order terms to better represent the true solution. Integrating the
solution over time produces a fluence that is shown to broadly agree with the
steady-state calculations in Section 2, provided a sufficient number of terms
in the sum and emission at the line-center frequency $\nu_{0}$. Integrating
the solution over frequency leads to a distribution of photon escape time,
which can be compared directly with Monte Carlo simulations. The sum over
eigenmodes produces an escape-time distribution that broadly captures the
behavior shown by Monte Carlo data—a rise at early times, transitioning to
exponential decay in the tail of the distribution. It is expected that the
accuracy of the rate constant for the tail of the distribution is limited by
the effect of the Doppler core, which can trap photons at high optical depths
until they diffuse outward in frequency, weighting the distribution toward
later times. This physics is not modeled by our solution for two reasons: 1)
our calculations ignored the Gaussian component of the Voigt line profile,
leaving the Lorentzian piece which is accurate only in the line wing, and 2)
knowing the core is not modeled accurately, we do not include a large enough
number of spatial eigenmodes in the sum to resolve it. However, an approximate
fitting function dependent on parameters $a$ and $\tau_{0}$ is found that
adequately represents the escape time distribution of the Monte Carlo results
within these constraints.

Our characterization of the escape time distribution leads to a possible
application of this work. Models of the interaction of stellar Ly$\alpha$ with
the upper atmosphere of exoplanets and the associated transmission spectrum
can be constructed with a treatment of resonant scattering in spherical
geometry (Huang et al., 2017; Yan et al., 2021). The Monte Carlo method can be
used for this problem, but is limited by its high computational demand for
large $\tau_{0}$ where there are many photon scatterings before escape. We
seek to develop a method to accelerate the radiative transfer calculation.

There are several methods that are commonly used to accelerate Monte Carlo
radiation transfer calculations, including core skipping methods (Auer, 1968;
Ahn et al., 2002) and hybrid diffusion methods (Smith et al., 2018). Another
approach with wide application is modified random walk methods, such as those
discussed in Fleck & Canfield (1984); Min et al. (2009); Robitaille (2010). In
this approach, an outgoing photon is randomly sampled on the surface of the
outgoing sphere by drawing its properties from distributions in outgoing
frequencies, directions, and escape times, based on solutions to the diffusion
equation. A method similar to this has been applied by Tasitsiomi (2006) to
Lyman $\alpha$ transfer using the Neufeld (1990) solution, but this solution
of course does not utilize the frequency-dependent boundary condition at the
surface of the sphere. Furthermore, to perform a full radiation hydrodynamic
simulation with Monte Carlo acceleration, it will be necessary to calculate
radiation forces within each cell due to Ly$\alpha$ transfer. Similar
calculations have been done in Weymann (1976) in plane-parallel geometry.
However, these solutions are limited to optical depths below $2.5\times
10^{3}$. For this work, it would be necessary to model line center optical
depths of up to 1 million or more.

##  5 Summary

We have examined previous solutions to Ly$\alpha$ transfer including resonant
scattering in the limit of large optical depth, noting that the separation of
variables and treatment of the boundary condition in Harrington (1973),
Neufeld (1990), Dijkstra et al. (2006) and others produces a discrepancy in
the outgoing spectrum as compared with Monte Carlo. Here, we have derived the
solution in spherical geometry with an appropriate treatment of the surface
boundary condition. The key result is that the errors in the previously-cited
works have been quantified via a correction term, $H_{\rm bc}$, which explains
an excess in flux at the spectral peak and a deficit in the line wing of the
calculated spectrum of Ly$\alpha$ radiation as compared with Monte Carlo. The
size of $H_{\rm bc}/H_{0}$ is of order unity when the spectral peaks are near
the Doppler core, and diminishes at larger $\tau_{0}$ following a
$(a\tau_{0})^{-1/3}$ scaling.

The time-dependent transfer equation for the impulsive source is solved
numerically with an eigenfunction expansion. We demonstrate that it agrees
with the steady-state spectrum for $x_{s}=0$ when integrated over time, though
its rate of numerical convergence is slow and requires a sum over many modes
to become accurate. The time-dependent solution is utilized to create wait-
time distributions for photons escaping the sphere of optically-thick hydrogen
gas. We compare the calculations from the time-dependent solution with Monte
Carlo for a sample of $\tau_{0}$, noting general agreement in the resulting
escape time distributions. The solution derived in our work here may be used
as the basis for a novel implementation of the modified random walk method,
which would accelerate Monte Carlo Ly$\alpha$ transfer at large optical depths
with potential applications in radiation hydrodynamic simulations of the
atmospheres of exoplanets.

This research was funded by NASA ATP grant 80NSSC18K0696, “Exoplanetary MHD
Outflows Driven by EUV Heating, Lyman alpha Radiation Forces and Stellar
Tides”. Resources supporting this work were provided by the NASA High-End
Computing (HEC) Program through the NASA Advanced Supercomputing (NAS)
Division at Ames Research Center. We thank the referee for a detailed report
that helped improve the presentation of our work.

##  Appendix A  derivation of the transfer equation

The problem is as follows. The radiative intensity $I=dE/(dAdtd\Omega d\nu)$
is the energy per perpendicular area $dA$, per time $dt$, per solid angle
$d\Omega$ and per frequency $d\nu$ (Rybicki & Lightman, 1986). The intensity
$I=I(\mathbf{x},t,\mathbf{n},\nu)$ will be considered a function of position
$\mathbf{x}$, time $t$, photon (unit) direction vector $\mathbf{n}$, and
cyclic frequency $\nu$. In the Eddington and two-stream approximations,
$I(\mathbf{x},\nu)\simeq
J(\mathbf{x},\nu)+3\mathbf{n}\cdot\mathbf{H}(\mathbf{x},\nu)$, where
$J=(1/4\pi)\int d\Omega I$ is the mean intensity and
$\mathbf{F}=4\pi\mathbf{H}=\int d\Omega\mathbf{n}I$ is the flux.

The transfer equation is (Rybicki & Lightman, 1986)

$\displaystyle\frac{1}{c}\frac{\partial I}{\partial
t}+\mathbf{n}\cdot\mathbf{\nabla}I$ $\displaystyle=$
$\displaystyle-\left(\alpha_{\rm sc}+\alpha_{\rm abs}\right)I+(1-p)j_{\rm
sc}+j_{\rm em}.$ (A1)

The scattering coefficient, or inverse mean free path to scattering, is

$\displaystyle\alpha_{\rm sc}$ $\displaystyle=$ $\displaystyle n_{\rm
sc}\,\frac{\pi
e^{2}}{m_{e}c}\,f\,\frac{\mathcal{H}(x,a)}{\sqrt{\pi}\Delta}=k\phi.$ (A2)

The absorption coefficient $\alpha_{\rm abs}$, or inverse mean free path to
true absorption, is a sum over species number density times absorption cross
section. Once the incoming photon has promoted the electron to an excited
state, the collisional de-excitation probability is $p$, and hence only a
fraction $1-p$ of the excitations lead to re-emission of photons.

Harrington (1973) first showed that the transfer equation for the mean
intensity $J$ will satisfy a Poisson equation involving space and frequency.
In this section we will briefly review the derivation of this equation
including photon destruction terms and an emission term.

“Hummer Case II-b” (Hummer, 1962) will be used for the redistribution
function, for which the incoming photon is absorbed by the atom according to
the natural broadening profile in the rest frame, then is re-emitted with a
dipole phase function
$g(\mathbf{n},\mathbf{n}^{\prime})=(3/16\pi)(1+[\mathbf{n}\cdot\mathbf{n}^{\prime}]^{2})$,
which is appropriate for a 1s-2p transition (Berestetskii et al., 1982), and
then is averaged over a Maxwell-Boltzmann distribution of speeds for the atom.
The result can be written

$\displaystyle j_{\rm sc}(\mathbf{x},\mathbf{n},\nu)$ $\displaystyle=$
$\displaystyle k\int\frac{d^{3}v}{\pi^{3/2}v_{\rm th}^{3}}e^{-v^{2}/v_{\rm
th}^{2}}\,\int d\Omega^{\prime}\int
d\nu^{\prime}\,g(\mathbf{n},\mathbf{n}^{\prime})$ (A3) $\displaystyle\times$
$\displaystyle\delta\left(\nu-\nu^{\prime}-\nu_{0}\mathbf{v}\cdot(\mathbf{n}-\mathbf{n}^{\prime})/c\right)\left(\frac{\Gamma/4\pi^{2}}{\left(\nu^{\prime}-\nu_{0}-\nu_{0}\mathbf{v}\cdot\mathbf{n}^{\prime}/c\right)^{2}+(\Gamma/4\pi)^{2}}\right)\,I(\mathbf{x},\mathbf{n}^{\prime},\nu^{\prime})$
$\displaystyle=$ $\displaystyle 4\pi k\int d\Omega^{\prime}\int
d\nu^{\prime}R(\mathbf{n},\nu;\mathbf{n}^{\prime},\nu^{\prime})I(\mathbf{x},\mathbf{n}^{\prime},\nu^{\prime}),$

which defines the Case II-b redistribution function

$\displaystyle R(\mathbf{n},\nu;\mathbf{n}^{\prime},\nu^{\prime})$
$\displaystyle=$
$\displaystyle\frac{g(\mathbf{n},\mathbf{n}^{\prime})}{4\pi}\int\frac{d^{3}v}{\pi^{3/2}v_{\rm
th}^{3}}e^{-v^{2}/v_{\rm
th}^{2}}\,\delta\left(\nu-\nu^{\prime}-\nu_{0}\mathbf{v}\cdot(\mathbf{n}-\mathbf{n}^{\prime})/c\right)\left(\frac{\Gamma/4\pi^{2}}{\left(\nu^{\prime}-\nu_{0}-\nu_{0}\mathbf{v}\cdot\mathbf{n}^{\prime}/c\right)^{2}+(\Gamma/4\pi)^{2}}\right),$
(A4)

(Unno, 1952; Hummer, 1962).

The integral of the redistribution function over outgoing and incoming
frequency are

$\displaystyle\int d\nu\ R(\mathbf{n},\nu;\mathbf{n}^{\prime},\nu^{\prime})$
$\displaystyle=$
$\displaystyle\frac{1}{4\pi}g(\mathbf{n},\mathbf{n}^{\prime})\phi(\nu^{\prime})$
(A5)

and

$\displaystyle\int d\nu^{\prime}\
R(\mathbf{n},\nu;\mathbf{n}^{\prime},\nu^{\prime})$ $\displaystyle=$
$\displaystyle\frac{1}{4\pi}g(\mathbf{n},\mathbf{n}^{\prime})\phi(\nu)$ (A6)

where the right hand side is the usual Voigt function, the thermal average of
the Lorentzian. The former result implies that the integrated source and sink
terms for scattering cancel for $p=0$. In addition, $4\pi
R(\mathbf{n},\nu;\mathbf{n}^{\prime},\nu^{\prime})/\phi(\nu^{\prime})$ is the
normalized distribution for the outgoing $\mathbf{n}$ and $\nu$ given the
incoming $\mathbf{n}^{\prime}$ and $\nu^{\prime}$.

This probability distribution can be used to define the moments of the
frequency shift (Osterbrock, 1962)

$\displaystyle\langle\delta\nu^{n}\rangle$ $\displaystyle=$
$\displaystyle\frac{\int d\nu^{\prime}(\nu^{\prime}-\nu)^{n}R}{\int
d\nu^{\prime}R}=\frac{1}{\phi(\nu)}\int\frac{d^{3}v}{\pi^{3/2}v_{\rm
th}^{3}}e^{-v^{2}/v_{\rm
th}^{2}}\,\left(\frac{\nu_{0}\mathbf{v}\cdot(\mathbf{n}^{\prime}-\mathbf{n})}{c}\right)^{n}\left(\frac{\Gamma/4\pi^{2}}{\left(\nu-\nu_{0}-\nu_{0}\mathbf{v}\cdot\mathbf{n}/c\right)^{2}+(\Gamma/4\pi)^{2}}\right),$
(A7)

which are functions of $\nu$, $\mathbf{n}$ and $\mathbf{n}^{\prime}$. These
integrals can be evaluated in terms of the dimensionless moments of the
parallel velocity distribution, defined as

$\displaystyle\langle u_{\parallel}^{n}\rangle(x,a)$ $\displaystyle=$
$\displaystyle\frac{a/\pi}{\mathcal{H}(x,a)}\int\frac{du_{\parallel}u_{\parallel}^{n}e^{-u_{\parallel}^{2}}}{(x-u_{\parallel})^{2}+a^{2}}.$
(A8)

The end results for the first and second moments are

$\displaystyle\langle\delta\nu\rangle$ $\displaystyle=$ $\displaystyle-\langle
u_{\parallel}\rangle\left(1-\mathbf{n}\cdot\mathbf{n}^{\prime}\right)\Delta$
(A9) $\displaystyle\langle\delta\nu^{2}\rangle$ $\displaystyle=$
$\displaystyle\left[\langle
u_{\parallel}^{2}\rangle\left(1-\mathbf{n}\cdot\mathbf{n}^{\prime}\right)^{2}+\frac{1}{2}\left(1-\left(\mathbf{n}\cdot\mathbf{n}^{\prime}\right)^{2}\right)\right]\Delta^{2}.$
(A10)

For small frequency shifts $\nu^{\prime}-\nu$, the incoming intensity may be
expanded as

$\displaystyle I(\mathbf{x},\mathbf{n}^{\prime},\nu^{\prime})$
$\displaystyle\simeq$ $\displaystyle
I(\mathbf{x},\mathbf{n}^{\prime},\nu)+\frac{\partial
I(\mathbf{x},\mathbf{n}^{\prime},\nu)}{\partial\nu}(\nu^{\prime}-\nu)+\frac{1}{2}\frac{\partial^{2}I(\mathbf{x},\mathbf{n}^{\prime},\nu)}{\partial\nu^{2}}(\nu^{\prime}-\nu)^{2}+...$
(A11)

and the Fokker-Planck expansion of $j_{\rm sc}$ is (Rybicki & dell’Antonio,
1994)

$\displaystyle j_{\rm sc}(\mathbf{x},\mathbf{n},\nu)$ $\displaystyle\simeq$
$\displaystyle 4\pi k\int d\Omega^{\prime}\int
d\nu^{\prime}R(\mathbf{n},\nu;\mathbf{n}^{\prime},\nu^{\prime})\left[I(\mathbf{x},\mathbf{n}^{\prime},\nu)+\frac{\partial
I(\mathbf{x},\mathbf{n}^{\prime},\nu)}{\partial\nu}(\nu^{\prime}-\nu)+\frac{1}{2}\frac{\partial^{2}I(\mathbf{x},\mathbf{n}^{\prime},\nu)}{\partial\nu^{2}}(\nu^{\prime}-\nu)^{2}\right]$
(A12) $\displaystyle=$ $\displaystyle k\phi(\nu)\int
d\Omega^{\prime}g\left[I(\mathbf{x},\mathbf{n}^{\prime},\nu)+\frac{\partial
I(\mathbf{x},\mathbf{n}^{\prime},\nu)}{\partial\nu}\langle\delta\nu\rangle+\frac{1}{2}\frac{\partial^{2}I(\mathbf{x},\mathbf{n}^{\prime},\nu)}{\partial\nu^{2}}\langle\delta\nu^{2}\rangle\right].$

To perform the angular integrals, the Eddington approximation for the angular
dependence is inserted with the following result

$\displaystyle j_{\rm sc}$ $\displaystyle=$ $\displaystyle k\phi
J-k\phi\Delta\langle u_{\parallel}\rangle\left(\frac{\partial
J}{\partial\nu}-\frac{6}{5}\mathbf{n}\cdot\frac{\partial\mathbf{H}}{\partial\nu}\right)+\frac{1}{2}\Delta^{2}k\phi\left[\frac{\partial^{2}J}{\partial\nu^{2}}\left(\frac{7}{5}\langle
u_{\parallel}^{2}\rangle+\frac{3}{10}\right)-\frac{12}{5}\langle
u_{\parallel}^{2}\rangle\mathbf{n}\cdot\frac{\partial^{2}\mathbf{H}}{\partial\nu^{2}}\right]$
(A13) $\displaystyle\simeq$ $\displaystyle k\phi J-k\phi\Delta\langle
u_{\parallel}\rangle\frac{\partial
J}{\partial\nu}+\frac{1}{2}\Delta^{2}k\phi\left(\frac{7}{5}\langle
u_{\parallel}^{2}\rangle+\frac{3}{10}\right)\frac{\partial^{2}J}{\partial\nu^{2}}.$

The first term in Equation (A13), $k\phi J$, represents re-emission of the
photon through de-excitation of the atom. It cancels the $-k\phi J$ term in
Equation (A1) that corresponds to excitation of the atom. The terms involving
frequency derivatives of $\mathbf{H}$, if carried through the calculation, end
up giving terms smaller than the largest terms by a factor of $1/x^{2}$, which
is small in the line wing. These terms are ignored from here onward.

If only scattering is included, the transfer equation becomes

$\displaystyle\frac{1}{c}\frac{\partial}{\partial
t}\left(J+3\mathbf{n}\cdot\mathbf{H}\right)+\mathbf{n}\cdot\mathbf{\nabla}\left(J+3\mathbf{n}\cdot\mathbf{H}\right)$
$\displaystyle=$
$\displaystyle-3k\phi\mathbf{n}\cdot\mathbf{H}-k\phi\Delta\langle
u_{\parallel}\rangle\frac{\partial
J}{\partial\nu}+\frac{1}{2}\Delta^{2}k\phi\left(\frac{7}{5}\langle
u_{\parallel}^{2}\rangle+\frac{3}{10}\right)\frac{\partial^{2}J}{\partial\nu^{2}}.$
(A14)

Integrating over angle and frequency then gives

$\displaystyle\frac{1}{c}\frac{\partial J(\mathbf{x})}{\partial
t}+\mathbf{\nabla}\cdot\mathbf{H}(\mathbf{x})$ $\displaystyle=$
$\displaystyle\int d\nu\left(-k\phi\Delta\langle
u_{\parallel}\rangle\frac{\partial
J}{\partial\nu}+\frac{1}{2}\Delta^{2}k\phi\left(\frac{7}{5}\langle
u_{\parallel}^{2}\rangle+\frac{3}{10}\right)\frac{\partial^{2}J}{\partial\nu^{2}}\right)$
(A15) $\displaystyle=$ $\displaystyle k\int d\nu
J\frac{\partial}{\partial\nu}\left(\phi\Delta\langle
u_{\parallel}\rangle+\frac{\partial}{\partial\nu}\left[\frac{1}{2}\phi\Delta^{2}\left(\frac{7}{5}\langle
u_{\parallel}^{2}\rangle+\frac{3}{10}\right)\right]\right),$

where $\mathbf{H}(\mathbf{x})$ is the frequency integrated flux, and
$\mathbf{\nabla}\cdot\mathbf{H}(\mathbf{x})=0$ if there are no sources or
sinks of radiation. Integration by parts has been used to factor $J$ out,
assuming each term goes to zero at infinity. The quantity inside parentheses
must be constant, and since each term should go to zero at infinity, that
constant is zero. Hence the first and second moments of the frequency shift
are related by

$\displaystyle\phi\Delta\langle u_{\parallel}\rangle$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial\nu}\left[\frac{1}{2}\phi\Delta^{2}\left(\frac{7}{5}\langle
u_{\parallel}^{2}\rangle+\frac{3}{10}\right)\right].$ (A16)

As an example, in the damping wing, the line profile can be approximated as

$\displaystyle\phi\simeq\frac{a}{\pi x^{2}\Delta}$ (A17)

with $\langle u_{\parallel}\rangle\simeq 1/x$ and $\langle
u_{\parallel}^{2}\rangle\simeq 1/2$, and so this identity is satisfied. The
scattering source function can then be rewritten

$\displaystyle j_{\rm sc}$ $\displaystyle\simeq$ $\displaystyle k\phi
J+\frac{1}{2}k\Delta^{2}\frac{\partial}{\partial\nu}\left[\phi\left(\frac{7}{5}\langle
u_{\parallel}^{2}\rangle+\frac{3}{10}\right)\frac{\partial
J}{\partial\nu}\right]\simeq k\phi
J+\frac{1}{2}k\Delta^{2}\frac{\partial}{\partial\nu}\left(\phi\frac{\partial
J}{\partial\nu}\right),$ (A18)

The following equations will use the approximations for the damping wing. Thus
far the transfer equation is

$\displaystyle\frac{1}{c}\frac{\partial}{\partial
t}\left(J+3\mathbf{n}\cdot\mathbf{H}\right)+\mathbf{n}\cdot\mathbf{\nabla}\left(J+3\mathbf{n}\cdot\mathbf{H}\right)$
$\displaystyle=$ $\displaystyle j_{\rm em}-\left(k\phi+\alpha_{\rm
abs}\right)\left(J+3\mathbf{n}\cdot\mathbf{H}\right)+(1-p)\left[k\phi
J+\frac{1}{2}k\Delta^{2}\frac{\partial}{\partial\nu}\left(\phi\frac{\partial
J}{\partial\nu}\right)\right]$ (A19) $\displaystyle\simeq$ $\displaystyle
j_{\rm em}-3(k\phi+\alpha_{\rm
abs})\mathbf{n}\cdot\mathbf{H}-\left(pk\phi+\alpha_{\rm
abs}\right)J+\frac{1}{2}k\Delta^{2}\frac{\partial}{\partial\nu}\left(\phi\frac{\partial
J}{\partial\nu}\right),$

where leading order dissipative terms were kept in the second equality. The
moment equations are

$\displaystyle\frac{1}{c}\frac{\partial J}{\partial
t}+\mathbf{\nabla}\cdot\mathbf{H}$ $\displaystyle=$ $\displaystyle j_{\rm
em}-\left(pk\phi+\alpha_{\rm
abs}\right)J+\frac{1}{2}k\Delta^{2}\frac{\partial}{\partial\nu}\left(\phi\frac{\partial
J}{\partial\nu}\right)$ (A20)

and

$\displaystyle\frac{1}{c}\frac{\partial\mathbf{H}}{\partial
t}+\frac{1}{3}\mathbf{\nabla}J$ $\displaystyle=$
$\displaystyle-\left(k\phi+\alpha_{\rm abs}\right)\mathbf{H}.$ (A21)

The $\partial\mathbf{H}/\partial t$ term may be dropped for slowly changing
sources. Assuming the coefficients are constant in space, these two equations
can be combined together to find

$\displaystyle\frac{1}{c}\frac{\partial J}{\partial
t}-\frac{1}{3(k\phi+\alpha_{\rm abs})}\nabla^{2}J$ $\displaystyle=$
$\displaystyle j_{\rm em}-\left(pk\phi+\alpha_{\rm
abs}\right)J+\frac{1}{2}k\Delta^{2}\frac{\partial}{\partial\nu}\left(\phi\frac{\partial
J}{\partial\nu}\right).$ (A22)

Making the change of variables to $d\sigma$ using

$\displaystyle d\sigma=\sqrt{\frac{2}{3}}\frac{d\nu}{\phi\Delta^{2}},$ (A23)

the equation can be rewritten in the standard form (Harrington, 1973)

$\displaystyle-3\left(\frac{k\phi+\alpha_{\rm abs}}{c}\right)\frac{\partial
J}{\partial
t}+\nabla^{2}J+\left(\frac{k}{\Delta}\right)^{2}\left(1+\frac{\alpha_{\rm
abs}}{k\phi}\right)\frac{\partial^{2}J}{\partial\sigma^{2}}$ $\displaystyle=$
$\displaystyle-3\left(k\phi+\alpha_{\rm abs}\right)j_{\rm
em}+3\left(k\phi+\alpha_{\rm abs}\right)\left(pk\phi+\alpha_{\rm
abs}\right)J.$ (A24)

For emission with luminosity $L$, with a delta function in space
$\delta^{3}(\mathbf{x}-\mathbf{x}_{s})$ at source position $\mathbf{x}_{s}$, a
delta function in frequency $\delta(\nu-\nu_{\rm s})$ at source frequency
$\nu_{\rm s}$, and isotropic in angles, the emission coefficient is

$\displaystyle j_{\rm em}$ $\displaystyle=$
$\displaystyle\frac{L}{4\pi}\delta^{3}(\mathbf{x}-\mathbf{x}_{s})\delta(\nu-\nu_{\rm
s}).$ (A25)

Multiplying by a factor $-3k\phi$, as appears in Equation (A24), the delta
function in $\nu$ becomes a delta function in $\sigma$ of the form

$\displaystyle-3k\phi j_{\rm em}$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{6}kL}{4\pi\Delta^{2}}\delta^{3}(\mathbf{x}-\mathbf{x}_{s})\delta(\sigma-\sigma_{\rm
s}),$ (A26)

where $\sigma_{s}=\sigma(\nu_{s})$.

##  Appendix B WKB Approximation

In this appendix, we will find a WKB solution for the differential equation
Equation (45) subject to the boundary conditions in Equations (28) and (29).
We will find this solution for the case $\sigma_{s}=0$, assuming $\sigma\geq
0$ for simplicity since the solution is symmetric about $0$.

Near the line center, there are two independent solutions which can be written

$\displaystyle J(n,\sigma,s)$ $\displaystyle=$ $\displaystyle
a_{1}\left(\frac{\sigma}{\sigma_{\rm
tp}}\right)^{1/2}J_{-3/4}\left(k_{c}\left[\frac{\sigma}{\sigma_{\rm
tp}}\right]^{2/3}\right)$ (B1) $\displaystyle+$ $\displaystyle
a_{2}\left(\frac{\sigma}{\sigma_{\rm
tp}}\right)^{1/2}J_{3/4}\left(k_{c}\left[\frac{\sigma}{\sigma_{\rm
tp}}\right]^{2/3}\right),$

where $J_{\alpha}(x)$ are Bessel functions of the first kind, $a_{1}$ and
$a_{2}$ are normalization constants, $\sigma_{tp}$ is given in Equation (46),
and $k_{c}=3\kappa_{n}\Delta\sigma_{\rm tp}/2k$. The solution far from line
center is

$\displaystyle J(n,\sigma,s)$ $\displaystyle=$ $\displaystyle
a_{3}\left(\frac{\sigma}{\sigma_{\rm
tp}}\right)^{1/6}Ai\left(k_{c}^{2/3}\left[\left(\frac{\sigma}{\sigma_{\rm
tp}}\right)^{2/3}-1\right]\right),$ (B2)

where $Ai(x)$ is the Airy function of the first kind. Again, $a_{3}$ is a
normalization constant. We can enforce the discontinuity in Equation (29) to
find

$\displaystyle a_{2}$ $\displaystyle=$
$\displaystyle-\frac{\Gamma(7/4)\sigma_{\rm
tp}}{2^{1/4}k_{c}^{3/4}}\frac{6^{1/2}}{8}n^{2}\frac{E}{kR^{3}}.$ (B3)

The discontinuity sets the value of $a_{2}$, and the corresponding term in
Equation (B1) is non-resonant. The discontinuity does not set the value of
$a_{1}$ as that solution does not have a term linear in $\sigma$. Matched
asymptotic expansions (Bender & Orszag, 1999) allows us to relate the
normalization coefficients $a_{1}$ and $a_{3}$. The matching conditions are

$\displaystyle\frac{a_{1}}{a_{2}}$ $\displaystyle=$
$\displaystyle\frac{\sin\left(\frac{\kappa_{n}\sigma_{\rm
tp}\Delta}{k}+\frac{\pi}{8}\right)}{\sin\left(\frac{\kappa_{n}\sigma_{\rm
tp}\Delta}{k}-\frac{\pi}{8}\right)}$ (B4) $\displaystyle\frac{a_{3}}{a_{2}}$
$\displaystyle=$
$\displaystyle\frac{k_{c}^{-1/3}}{\sin\left(\frac{\kappa_{n}\sigma_{\rm
tp}\Delta}{k}-\frac{\pi}{8}\right)}.$ (B5)

By setting the denominator of this expression equal to zero and solving for
the eigenfrequency $\gamma_{nm}$ contained in $\sigma_{\rm tp}$, we find the
dispersion relation in Equation (47).

While the solution above is specific to the case $\sigma_{s}=0$, we can extend
this approach to understand the case where $\sigma_{s}\neq 0$. In the interval
$\sigma\in(0,\sigma_{s})$, both the $Ai$ and $Bi$ Airy function solutions must
be included, while for $\sigma\in(-\infty,0)$ only $Ai$ is finite. The $Bi$
term causes the asymmetry. As $\gamma\rightarrow\gamma_{nm}$, however, it is
small compared to the $Ai$ term, and hence the eigenfunctions are symmetric.

## References

  * Adams (1972) Adams, T. F. 1972, ApJ, 174, 439, doi: 10.1086/151503
  * Adams (1975) —. 1975, ApJ, 201, 350, doi: 10.1086/153891
  * Ahn et al. (2002) Ahn, S.-H., Lee, H.-W., & Lee, H. M. 2002, ApJ, 567, 922, doi: 10.1086/338497
  * Auer (1968) Auer, L. H. 1968, ApJ, 153, 783, doi: 10.1086/149705
  * Bender & Orszag (1999) Bender, C. M., & Orszag, S. A. 1999, Advanced mathematical methods for scientists and engineers I: asymptotic methods and perturbation theory.
  * Berestetskii et al. (1982) Berestetskii, V. B., Lifshitz, E. M., & Pitaevskii, V. B. 1982, Quantum Electrodynamics, 2nd ed. (Butterworth-Heinemann Press)
  * Bourrier et al. (2017a) Bourrier, V., Ehrenreich, D., King, G., et al. 2017a, A&A, 597, A26, doi: 10.1051/0004-6361/201629253
  * Bourrier et al. (2017b) Bourrier, V., Ehrenreich, D., Wheatley, P. J., et al. 2017b, A&A, 599, L3, doi: 10.1051/0004-6361/201630238
  * Bourrier et al. (2017c) Bourrier, V., Ehrenreich, D., Allart, R., et al. 2017c, A&A, 602, A106, doi: 10.1051/0004-6361/201730542
  * Bourrier et al. (2018) Bourrier, V., Lecavelier des Etangs, A., Ehrenreich, D., et al. 2018, A&A, 620, A147, doi: 10.1051/0004-6361/201833675
  * Bourrier et al. (2021) Bourrier, V., dos Santos, L. A., Sanz-Forcada, J., et al. 2021, arXiv e-prints, arXiv:2103.09864.  https://arxiv.org/abs/2103.09864
  * Dijkstra (2019) Dijkstra, M. 2019, Saas-Fee Advanced Course, 46, 1, doi: 10.1007/978-3-662-59623-4_1
  * Dijkstra et al. (2006) Dijkstra, M., Haiman, Z., & Spaans, M. 2006, ApJ, 649, 14, doi: 10.1086/506243
  * Ehrenreich et al. (2012) Ehrenreich, D., Bourrier, V., Bonfils, X., et al. 2012, A&A, 547, A18, doi: 10.1051/0004-6361/201219981
  * Ehrenreich et al. (2015) Ehrenreich, D., Bourrier, V., Wheatley, P. J., et al. 2015, Nature, 522, 459, doi: 10.1038/nature14501
  * Fleck & Canfield (1984) Fleck, J. A., J., & Canfield, E. H. 1984, Journal of Computational Physics, 54, 508, doi: 10.1016/0021-9991(84)90130-X
  * García Muñoz et al. (2020) García Muñoz, A., Youngblood, A., Fossati, L., et al. 2020, ApJ, 888, L21, doi: 10.3847/2041-8213/ab61ff
  * Harrington (1973) Harrington, J. P. 1973, MNRAS, 162, 43, doi: 10.1093/mnras/162.1.43
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362, doi: 10.1038/s41586-020-2649-2
  * Huang et al. (2017) Huang, C., Arras, P., Christie, D., & Li, Z.-Y. 2017, ApJ, 851, 150, doi: 10.3847/1538-4357/aa9b32
  * Hummer (1962) Hummer, D. G. 1962, MNRAS, 125, 21, doi: 10.1093/mnras/125.1.21
  * Lao & Smith (2020) Lao, B.-X., & Smith, A. 2020, MNRAS, 497, 3925, doi: 10.1093/mnras/staa2198
  * Lavie et al. (2019) Lavie, B., Bourrier, V., & Ehrenreich, D. 2019, in EPSC-DPS Joint Meeting 2019, Vol. 2019, EPSC–DPS2019–1928
  * Lecavelier des Etangs et al. (2012) Lecavelier des Etangs, A., Bourrier, V., Wheatley, P. J., et al. 2012, A&A, 543, L4, doi: 10.1051/0004-6361/201219363
  * Min et al. (2009) Min, M., Dullemond, C. P., Dominik, C., de Koter, A., & Hovenier, J. W. 2009, A&A, 497, 155, doi: 10.1051/0004-6361/200811470
  * Neufeld (1990) Neufeld, D. A. 1990, ApJ, 350, 216, doi: 10.1086/168375
  * Osterbrock (1962) Osterbrock, D. E. 1962, ApJ, 135, 195, doi: 10.1086/147258
  * Robitaille (2010) Robitaille, T. P. 2010, A&A, 520, A70, doi: 10.1051/0004-6361/201015025
  * Roy et al. (2010) Roy, I., Shu, C.-W., & Fang, L.-Z. 2010, ApJ, 716, 604, doi: 10.1088/0004-637X/716/1/604
  * Rybicki & dell’Antonio (1994) Rybicki, G. B., & dell’Antonio, I. P. 1994, ApJ, 427, 603, doi: 10.1086/174170
  * Rybicki & Lightman (1986) Rybicki, G. B., & Lightman, A. P. 1986, Radiative Processes in Astrophysics
  * Seon & Kim (2020) Seon, K.-i., & Kim, C.-G. 2020, ApJS, 250, 9, doi: 10.3847/1538-4365/aba2d6
  * Smith et al. (2017) Smith, A., Bromm, V., & Loeb, A. 2017, MNRAS, 464, 2963, doi: 10.1093/mnras/stw2591
  * Smith et al. (2015) Smith, A., Safranek-Shrader, C., Bromm, V., & Milosavljević, M. 2015, MNRAS, 449, 4336, doi: 10.1093/mnras/stv565
  * Smith et al. (2018) Smith, A., Tsang, B. T. H., Bromm, V., & Milosavljević, M. 2018, MNRAS, 479, 2065, doi: 10.1093/mnras/sty1509
  * Tasitsiomi (2006) Tasitsiomi, A. 2006, ApJ, 645, 792, doi: 10.1086/504460
  * Tomaselli & Ferrara (2021) Tomaselli, G. M., & Ferrara, A. 2021, MNRAS, 504, 89, doi: 10.1093/mnras/stab876
  * Unno (1952) Unno, W. 1952, PASJ, 4, 100
  * Vidal-Madjar et al. (2003) Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J. M., et al. 2003, Nature, 422, 143, doi: 10.1038/nature01448
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Waalkes et al. (2019) Waalkes, W. C., Berta-Thompson, Z., Bourrier, V., et al. 2019, AJ, 158, 50, doi: 10.3847/1538-3881/ab24c2
  * Weymann (1976) Weymann, R. 1976, ApJ, 208, 286, doi: 10.1086/154605
  * Xu et al. (2011) Xu, W., Wu, X.-P., & Fang, L.-Z. 2011, MNRAS, 418, 853, doi: 10.1111/j.1365-2966.2011.19539.x
  * Yan et al. (2021) Yan, D., Guo, J., Huang, C., & Xing, L. 2021, ApJ, 907, L47, doi: 10.3847/2041-8213/abda41