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

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.


INTRODUCTION
Given the abundance of hydrogen in the universe, the Lyman α (Lyα) line is an important component of radiation fields in a wide range of astrophysical settings. Lyα 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α 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α: the star, and recombinations in the planet's atmosphere. Lyα may ionize atoms and dissociate molecules, as well as exert pressure forces that drive an outflow (Bourrier et al. 2018). Lyα 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α may undergo "resonant scattering".
Hubble Space Telescope (HST) observations with the STIS have found large Lyα transit depths around a handful of exoplanets (Vidal-Madjar et al. 2003;Lecavelier des Etangs et al. 2012;Ehrenreich et al. 2012Ehrenreich et al. , 2015Bourrier et al. 2017aBourrier et al. ,b,c, 2018Waalkes 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α 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 ∼ 10 µ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 τ 0 ∼ 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α 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, τ 0 , of the domain (Adams 1972). Near the base of the atomic layer, the line center optical depth is ∼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(τ, σ) = θ(τ )j(σ) in spatial variable τ and frequency variable σ (their Equations 16 and 23). The solutions for the eigenfunctions θ(τ ) and j(σ) then depend explicitly on the separation constant λ. In order to satisfy the boundary conditions, the separation constant is shown to satisfy an eigenvalue equation of the form λ tan(λB) = 3 2 φ∆, where 2B is the slab optical depth at line center, φ is the line profile, and ∆ 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 λ n B π(n − 1/2), which gives zero mean intensity at the surface. Their Equation 34 subsequently allows λ 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, steadystate solutions to this problem are not always sufficient to describe all the physics of Lyα transport. Time-variable, optically-thick environments necessitate a time-dependent solution to include the dynamic effects of Lyα transfer. These include the optical afterglow of gamma-ray bursts (Roy et al. 2010) and Lyα sources redshifted by cosmological expansion (Xu et al. 2011), among others.

STEADY-STATE SOLUTION
Consider a sphere of radius R with uniform density n sc , luminosity L, and line-center optical depth τ 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 ν near the line center frequency ν 0 are considered. The photon frequency of the source is ν s . The Doppler width is ∆ = ν 0 v th /c, where v th = 2k B T /m H is the thermal speed of hydrogen atoms of mass m H and temperature T , and c is the speed of light. The photon frequency in Doppler units is x = (ν − ν 0 )/∆, and x s = (ν s − ν 0 )/∆ is the corresponding source frequency. For upper-state de-excitation rate Γ, the ratio of natural to Doppler broadening is a = Γ/(4π∆). For the Lyα transition and T=10 4 K, a = 4.72 × 10 −4 . H(x, a) is the Voigt function, and the Voigt line profile is φ = H(x, a)/( √ π∆), which is normalized as dν φ(ν) = 1. The line center optical depth is τ 0 = kR/( √ π∆), where k = n sc π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α (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 ( 2) where J is the mean intensity, the spatial variable is x, and x s is the position of the source. We will consider only the case where x s = 0. Following Harrington (1973), we have used a change of variables in photon frequency from x to σ, where the approximation is applicable in the damping wing. From Equation (A17), the line profile is then approximately In Equation (2), σ s ≡ σ(x s ) is the photon frequency of the source. σ s is interchangeable with x s and ν s in Doppler widths or Hz, respectively. Balancing the two terms on the left-hand side of Equation (2) gives σ ∼ τ 0 , or x peak ∼ (aτ 0 ) 1/3 . The boundary condition of no incoming intensity at the surface (Rybicki & Lightman 1986) is A solution for the mean intensity J d which is divergent at the origin and σ = σ 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: 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 d (R, σ)/H d (R, σ) ∼ aτ 0 /x 2 ∼ (aτ 0 ) 1/3 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 closedform 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 σ s away from line center. The result can be written as a sum over spatial modes, and where κ n = nπ/R. These can be summed to give the closed form expressions and These solutions agree with Equations (6) and (7) when the arguments of the trigonometric and hyperbolic functions are small. Again J 0 H 0 , except near r = R, where it goes to zero. The flux at r = R can be written 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 σ. 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 bc which allows the boundary condition J = √ 3H to be satisfied at r = R. The total solution is and H(r, σ) = H 0 (r, σ) + H bc (r, σ).
The additional term J bc must then be a solution of the homogeneous equation with no delta function source term, and it must allow the boundary conditions to be satisfied at the surface. Since J 0 (R, σ) = 0, the surface boundary condition becomes Inserting a frequency dependence J bc ∝ e isσ , 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 where A(s) is the Fourier amplitude. Inserting Equation (17) into Equation (16) leads to the following equation for the Fourier amplitudes, Discretization of Equation (18) for frequency variables σ i and wavenumbers s j leads to a set of coupled linear equations for the A(s j ). We use equally-spaced points δσ = 2σ max /(N − 1) and δs = 2π/(N δσ), where N is the number of points for each grid. The maximum frequency is set as σ max = constant × τ 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 σ max = 60τ 0 were sufficient. Given the Fourier amplitudes A(s), J bc is computed using Equation (17), and the flux is given by The Bessel functions are finite at the center and rise steeply toward the surface when kRs/∆ 1.
2.1. Scaling with Line Center Optical Depth τ 0 We now estimate the scaling of H bc with τ 0 . In the limit J bc H bc , we find that J bc ≈ √ 3H 0 from Equation (16). We estimate H bc from J bc using Equation (19) as where we have used s ∼ 1/σ ∼ 1/τ 0 so that At large τ 0 , it is expected that the correction term will be small, but it will become increasingly important as τ 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 τ 0 . The value of x at which the Doppler and Lorentzian components of the line profile are equal is x cw = 3.3. Setting x cw = x peak and solving for τ 0 gives the value at which the peak of the spectrum falls at the Doppler core boundary, which is τ cp ≈ 10 5 ("core-peak" optical depth). Hence H bc /H 0 ∼(τ cp /τ 0 ) 1/3 is large at τ 0 ≤ τ cp and decreases relatively slowly as τ 0 increases. Additionally, the optical depth at x peak is proportional to (aτ 0 ) 1/3 , so photons here become optically thin when aτ 0 ∼1.

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 τ 0 , being restricted only by the computational demand, which grows proportionally to the number of photons used and τ 0 . For each simulation, a total of ∼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 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 This is normalized so that P (x)dx = 1. Since H(R, x) is per dν, a factor of ∆ gives the expression the correct units.
In Figure 1, the Monte Carlo spectrum is shown along with that of the solutions H d , H 0 , and H 0 + H bc for an optical depth of τ 0 = 10 7 and photons emitted at line center x s = 0. Note that the errorbars shown on the Monte Carlo data points are proportional to √ N , with N being the photon count in each frequency bin, since the photons are all equally weighted. The H 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 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 bc term, 0.00 0.02 0.04 since the Monte Carlo represents the "true" solution, H, and H bc = H − H 0 . It is evident that the divergent solution H 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.
The size of H bc is dependent on τ 0 . H bc is significant even at τ 0 ∼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 τ 0 = 10 5 , 10 6 , and 10 7 . From Equation (21), the size of the term H bc should become smaller with larger optical depths, following a (aτ 0 ) −1/3 scaling. Indeed, agreement between Equation (12) and the Monte Carlo points in Figure 2 improves as τ 0 increases, with H bc providing a fractionally smaller correction to H 0 . One factor of (aτ 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 τ 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 τ 0 = 10 5 and, to a lesser extent, τ 0 = 10 6 . 0.0 0.5 1.0 Next, we show P (x) for x s = 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 |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 τ 0 = 10 7 . The asymmetry of the spectrum is slight for x s = 6 where τ (x s ) = 77, but is larger for x s = 12 outside the line core where τ (x s ) = 19. It is seen here that the difference between the Monte Carlo data and H 0 becomes larger as x s increases. Thus, for large |x s |, inclusion of H 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 τ 0 = 10 6 rather than 10 7 . The difference between the left and right side of the escaping spectrum is now substantial since (aτ 0 ) −1/3 increased by a factor of ∼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 shows how the correction P bc (x), Equation (22) with H → H bc , scales with both τ 0 and x s . In this figure, σ s is shifted by integer multiples of τ 0 (Equation A23) in each panel such that the source falls near the peak of the spectrum for each τ 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 bc relative to H 0 should become smaller with larger optical depths following (aτ 0 ) −1/3 . This factor has been scaled out of the figure such that solutions for different τ 0 and the same σ s should show close agreement in scale on the figure's vertical axis if the relation holds. Indeed, the scaled solutions converge as τ 0 becomes larger, indicating agreement with the (aτ 0 ) −1/3 scaling. The remaining discrepancy present in the vertical axis for fixed σ s results from x peak becoming close to x cw ; at lower τ 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τ 0 ) −1/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Ωdν) on timescales comparable to the light-crossing time t 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τ 0 1.

Derivation of the time-dependent solution
The emissivity for an impulsive source with energy E, source position x s , and frequency ν s is derived in Appendix A. Considering a photon source at x s = 0, we have (Equation A25) The resulting equation for J(r, σ, t) is 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, σ, t) is then Figure 4. The same as Figure 3, but at a lower optical depth τ0 = 10 6 . The shift xs is a much larger fraction of the distance to the spectral peak (aτ0) 1/3 , and thus the asymmetry in the spectrum is much larger. The optical depth at the source frequency is τs = 10 6 , 7.7, and 1.9 for xs = 0, 6, and 12, respectively.
This approximate solution implies a boundary condition at large |σ| where a negative sign is taken for large +σ and a positive sign is taken for large −σ to choose the finite solution as |σ| → ∞. Numerical integrations are performed inward toward σ s over several domains: from large |σ| to σ s , from large |σ| to 0, and from 0 to σ s , depending on whether σ s is positive or negative. If σ s = 0, just two integrations are performed inward from large |σ| to 0. Initial values for integration are obtained either by setting J = 1 and dJ/dσ from Equation (32) at large |σ| or by matching J and dJ/dσ at 0. This gives J and J on either side of σ s , where a prime indicates the derivative ∂/∂σ. By enforcing the matching conditions, Equations (28) and (29), the eigenfunctions J(n, σ, ω) are obtained over the domain of photon frequencies σ. 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, σ, 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 = 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 γ ≡ iω, which is real and positive for damped solutions. At the eigenvalues γ = γ nm , the response J(n, σ, ω) is resonant. We find that near these γ nm poles an approximate expression for the resonant response of the eigenfunctions is where C(γ, σ) varies slowly in γ. If the ω-integral in Equation (25) could be closed at infinity and evaluated using the residue theorem, the result would be Summing over all spatial modes n and over all eigenmodes m for a given n, we obtain This ansatz captures the contributions from n × m simple poles. Taking a derivative with respect to r and evaluating at the surface r = R, we use to obtain the flux, which is Multiplying by 4πR 2 gives the energy per time per frequency emerging from the sphere to be Integrating over time yields a factor 1/γ nm , and by integrating over dν we find This non-trivial "sum rule" provides a check on the values of γ nm and J nm (σ). This expression can also be written as where the contribution of each mode is 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 ∼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).

Numerical calculation
We seek now to calculate the eigenmodes J nm (σ) and eigenfrequencies γ 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, σ, −iγ) over the array σ. We call this response f , and use the index j to represent the value of the response at discrete points γ j over a range of γ. In places where f j > f j−1 and f j > f j+1 , we have bracketed a resonance that occurs in the interval (γ j−1 , γ j+1 ). To refine the value of the eigenfrequency before continuing the sweep in γ, we evaluate f j−1 , f j , and f j+1 at the points (γ j−1 , γ j , γ j+1 ). Assuming the form in Equation (33), a guess at the correct eigenvalue γ nm can be calculated by linear interpolation from where The error of the current guess is |γ guess −γ j |. This error is reduced iteratively by replacing initial points (γ j−1 , γ j , γ j+1 ) with closer estimates while the size of the response grows as the resonance is approached. After iterating an eigenvalue γ nm to convergence, we now find the corresponding eigenfunction J nm (σ). We evaluate Equation 33 at two points γ 1 and γ 2 near the resonance, subtracting them and solving for J nm (σ) to find where C(γ, σ) has cancelled in the difference. The form of a single eigenmode J nm (σ) is oscillatory out to some turning point, σ 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 σ s in Equation (27) and examining the resulting homogeneous differential equation. We obtain 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 Thus, to ensure accuracy in each term of Equation (35), the bounds of σ must be set sufficiently far outside of σ tp such that the function is small at the edges. The scale of an e-folding in J nm (σ) is k/(κ n ∆) = τ 0 /( √ πn), so a grid of σ 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 (σ) to illustrate their relative scales for different m at a fixed spatial eigenmode n. The overall scale of the J nm (σ) are set by the factor E/(kR 3 ) with E arbitrarily set to 1. For H atoms with T = 10 4 K and τ 0 = 10 7 , an eigenfunction has typical size ∼Ea/ R 2 ∆ = 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. The values of the γ nm can be described approximately with Equation (47). Their values depend on m, n, and other physical parameters according to γ nm = 2 −1/3 π 13/6 n 4/3 m − 7 8 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 γ −1 nm by two orders of magnitude. When sweeping through to find resonances, Equation (47) is used to set the scale of the sweep points γ j to ensure no γ nm are missed. The close agreement with the analytic expression shown in Figure 7 indicate that the numerical solutions are accurate.

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 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, γ 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 The first term represents the early-time distribution, which then transitions to an exponential falloff past a point ct diff /R = (aτ 0 ) 1/3 , where t 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 τ 0 . It is shown that the time constant of exponential decay in fitted Monte Carlo escape time distributions converges with γ −1 11 at sufficiently high τ 0 , following a t ∝ (aτ 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τ 0 ) 1/3 . At lower τ 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 τ 0 grows, the effect of core scattering becomes smaller and the approximations hold, agreeing better with the expected (aτ 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 τ 0 , and though these points are not shown at τ 0 = 10 8 and 10 9 due to computational expense, it is expected that the fractional error between the Monte Carlo and the (aτ 0 ) 1/3 scaling would be less than 2% at τ 0 = 10 9 .
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 Integrating over ν then gives unity as required by the sum rule in Equation (40). In Figure 9, the fluence for x s = 0 and τ 0 = 10 7 is shown for a sum up to n = 20 and m = 500, labelled "Timeintegrated", 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.  ). Fluence is the radiation flux integrated over time. Steady-state and timeintegrated spectra for n = 1, ..., 20 and m = 1, ..., 500 are shown with xs = 0 and τ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 τ 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  (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 (xs = 0).
Carlo. The error in this rate constant is a function of τ 0 since the effect from the Doppler core is greatest when it extends into the peak of the spectrum.

Steady-State Source
A primary goal of this work is to present a solution for resonant scattering of photons near the line-center frequency ν 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 bc , which allows the boundary condition J = √ 3H 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 bc , scales as H 0 (aτ 0 ) −1/3 . Thus, for large aτ 0 , only a small correction to H 0 is needed, while larger errors are present in calculations performed at lower aτ 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 φ ≈ a/(πx 2 ∆), our solutions do not accurately model the Doppler core of the Lyα line. Because the peak of the spectral energy distribution of escaping photons is x peak ∼(aτ 0 ) 1/3 , calculations performed at small aτ 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 frequencydependent boundary conditions improves the accuracy of these analytic solutions for aτ 0 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α spectrum emergent from a sphere of uniform optical depth is shown for τ 0 = 10 5 , 10 6 , and 10 7 at a temperature of T = 10 K, corresponding to a = 1.5 × 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 τ 0 = 10 5 and 10 6 . Another example is in Smith et al. (2015), where the peak excess in the Lyα spectrum is particularly noticeable for line center optical depths of τ 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τ 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 τ 0 = 10 7 for T= 10 4 K, corresponding to a = 4.72 × 10 −4 . We note that our correction term H 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.

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 σ is solved to find the flux at the surface of the sphere as a function of t and ν. 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 ν 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 τ 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α 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 τ 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 α 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α transfer. Similar calculations have been done in Weymann (1976) in plane-parallel geometry. However, these solutions are limited to optical depths below 2.5 × 10 3 . For this work, it would be necessary to model line center optical depths of up to 1 million or more.

SUMMARY
We have examined previous solutions to Lyα 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 bc , which explains an excess in flux at the spectral peak and a deficit in the line wing of the calculated spectrum of Lyα radiation as compared with Monte Carlo. The size of H bc /H 0 is of order unity when the spectral peaks are near the Doppler core, and diminishes at larger τ 0 following a (aτ 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 τ 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α transfer at large optical depths with potential applications in radiation hydrodynamic simulations of the atmospheres of exoplanets. The problem is as follows. The radiative intensity I = dE/(dAdtdΩdν) is the energy per perpendicular area dA, per time dt, per solid angle dΩ and per frequency dν (Rybicki & Lightman 1986). The intensity I = I(x, t, n, ν) will be considered a function of position x, time t, photon (unit) direction vector n, and cyclic frequency ν. In the Eddington and two-stream approximations, I(x, ν) J(x, ν) + 3n · H(x, ν), where J = (1/4π) dΩI is the mean intensity and F = 4πH = dΩnI is the flux.
The scattering coefficient, or inverse mean free path to scattering, is The absorption coefficient α 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(n, n ) = (3/16π)(1 + [n · n ] 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 which defines the Case II-b redistribution function R(n, ν; n , ν ) = g(n, n ) 4π (Unno 1952;Hummer 1962). The integral of the redistribution function over outgoing and incoming frequency are dν R(n, ν; n , ν ) = 1 4π g(n, n )φ(ν ) and dν R(n, ν; n , ν ) = 1 4π g(n, n )φ(ν) 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πR(n, ν; n , ν )/φ(ν ) is the normalized distribution for the outgoing n and ν given the incoming n and ν . This probability distribution can be used to define the moments of the frequency shift (Osterbrock 1962) δν n = dν (ν − ν) n R dν R = 1 φ(ν) which are functions of ν, n and n . These integrals can be evaluated in terms of the dimensionless moments of the parallel velocity distribution, defined as u n (x, a) = a/π H(x, a) du u n e −u 2 (x − u ) 2 + a 2 . (A8) The end results for the first and second moments are For small frequency shifts ν − ν, the incoming intensity may be expanded as I(x, n , ν ) I(x, n , ν) + ∂I(x, n , ν) ∂ν (ν − ν) + 1 2 ∂ 2 I(x, n , ν) ∂ν 2 (ν − ν) 2 + ...
To perform the angular integrals, the Eddington approximation for the angular dependence is inserted with the following result j sc = kφJ − kφ∆ u ∂J ∂ν − 6 5 n · ∂H ∂ν + 1 2 ∆ 2 kφ ∂ 2 J ∂ν 2 7 5 u 2 + 3 10 − 12 5 u 2 n · ∂ 2 H ∂ν 2 kφJ − kφ∆ u ∂J ∂ν + 1 2 ∆ 2 kφ 7 5 u 2 + 3 10 The first term in Equation (A13), kφJ, represents re-emission of the photon through de-excitation of the atom. It cancels the −kφJ term in Equation (A1) that corresponds to excitation of the atom. The terms involving frequency derivatives of 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 1 c ∂ ∂t (J + 3n · H) + n · ∇ (J + 3n · H) = −3kφn · H − kφ∆ u ∂J ∂ν + 1 2 ∆ 2 kφ 7 5 u 2 + 3 10 Integrating over angle and frequency then gives where H(x) is the frequency integrated flux, and ∇ · H(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 φ∆ u = − ∂ ∂ν 1 2 φ∆ 2 7 5 u 2 + 3 10 .
As an example, in the damping wing, the line profile can be approximated as φ a πx 2 ∆ (A17) with u 1/x and u 2 1/2, and so this identity is satisfied. The scattering source function can then be rewritten j sc kφJ + 1 2 k∆ 2 ∂ ∂ν φ 7 5 u 2 + 3 10 The following equations will use the approximations for the damping wing. Thus far the transfer equation is 1 c ∂ ∂t (J + 3n · H) + n · ∇ (J + 3n · H) = j em − (kφ + α abs ) (J + 3n · H) + (1 − p) kφJ + 1 2 k∆ 2 ∂ ∂ν φ ∂J ∂ν j em − 3(kφ + α abs )n · H − (pkφ + α abs ) J + 1 2 k∆ 2 ∂ ∂ν φ ∂J ∂ν , where leading order dissipative terms were kept in the second equality. By setting the denominator of this expression equal to zero and solving for the eigenfrequency γ nm contained in σ tp , we find the dispersion relation in Equation (47). While the solution above is specific to the case σ s = 0, we can extend this approach to understand the case where σ s = 0. In the interval σ ∈ (0, σ s ), both the Ai and Bi Airy function solutions must be included, while for σ ∈ (−∞, 0) only Ai is finite. The Bi term causes the asymmetry. As γ → γ nm , however, it is small compared to the Ai term, and hence the eigenfunctions are symmetric.