Refraction of Line and Continuum Light in Exoplanet Atmospheres

Exoplanet transit spectra are calculated including the effect of atmospheric opacity and refractive lightbending. While previous studies considered the case of continuum light, here the effect of an atomic resonance line is included. The model assumes a clear atmosphere and includes H and He, which contribute static polarizability and Rayleigh scattering, as well as the Na D doublet, which contributes dynamic polarizability and a resonant cross section. The image locations and magnifications are found using the lens equation. The model including lightbending is compared to the standard model in which the light travels on straight lines. It is found that near the line center, where the polarizability is large, bending angles are nevertheless small since the optical depth τ = 1 trajectory is at such a high altitude where the particle density is low. Moving away from the line center, the Na D resonance dominates the opacity over ∼400 Å, and over most of this wavelength range the polarizability is dominated by hydrogen and helium and is nearly wavelength-independent. However, the density of the τ = 1 trajectory is wavelength-dependent, and hence the bending angle increases strongly away from the line center. The wavelength-dependent flux deviation between the straight-line and lightbending models occurs at the level of ΔF λ /F 0 ∼ 10−5 for a planet at orbital separation a = 10 R ⊙.


Introduction
When an extrasolar planet transits the parent star, the diminution of stellar flux may be used to infer the presence of atoms and molecules in the planet's atmosphere through the wavelength-dependent planetary radius (e.g., Burrows et al. 2007).Most previous studies have constructed transmission spectra assuming the stellar light travels in a straight line through the planet's atmosphere, ignoring the effect of refractive lightbending.While lightbending has been shown to be a small effect (e.g., Hubbard et al. 2001), it may nevertheless be detectable at sufficiently high spectral and time resolution (Hui & Seager 2002) and may give additional information about the planet's atmosphere, such as the ratio of atmospheric scale height to radius.
Previous studies of refractive lensing have illuminated the possibility of strong lensing and large deflection angles for trajectories that travel deep in the planet's atmosphere.This may give rise to images of the stellar flux that curves around the limb of the planet (Hui & Seager 2002;Sidis & Sari 2010).The observability of these strongly bent trajectories is suppressed when absorption is included, as the optical depth τ  1 trajectories experience more modest bending angles due to the lower particle densities.Alp & Demory (2018) have made a recent comparison of these models to Kepler data, finding no detectable signatures of lensing when averaged over the broad Kepler band.
While previous studies have focused on continuum light, for which the wavelength-dependence of the absorption cross section and refractive index are slowly varying, here the role of resonance lines is studied by including the strong increase near resonance for both quantities.The objective is to understand if refraction may create features in the spectrum or lightcurve, which are confined to a small wavelength interval and do not require comparisons of the flux across the entire optical band.
The plan of the paper is as follows.Section 2 contains the prescription for the index of refraction and extinction cross section.Sections 3, 4, and 5 contain a review of the bending angle, lens equation, and image locations for the isothermal atmosphere.Results for the transmission spectrum and lightcurves are given in Section 6. Conclusions are given in Section 7.

Cross Section and Refractive Index
The dispersion relation of electromagnetic radiation in a dielectric medium is N = ck/ω, where N is the index of refraction, k = 2π/λ is the wavenumber, and ω = 2πν = 2πc/λ is the angular frequency.The index of refraction is related to the polarizability of the particles as (Jackson 1975) where n i and α i are the number density and polarizability of each species and it will be assumed that |N − 1| = 1 in this work.The total number density of the gas is n, and the mixing ratios are x and ξ He ; 0.16.Over the optical range of wavelengths, dynamic effects are a 1%-10% effect (Born & Wolf 1999) and will be ignored here.
Refraction causes light to bend away from regions of higher phase speed and toward regions of lower phase speed.In terms of the density n, the phase speed , and hence, in the case of "normal dispersion" Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
with α > 0, light bends toward the planet where density is higher.For the case of "anomalous dispersion" with α < 0, light bends away from the planet toward lower densities.
Ignoring the frequency dependence, the static polarizability of Na is α Na,static = 24.1 × 10 −24 cm 3 , which is larger than that of H and He.The large polarizability is offset though by the low (solar) abundance of Na, ξ Na ; 3.4 × 10 −6 , and the nonresonant contribution of Na to the polarizability is smaller by a factor of ∼10 5 .However, the resonance allows an increase in the response, larger by a factor ∼c/v th ∼ 10 5 near the line center, where v th is the thermal velocity of the Na atoms.Hence, refraction by H and He will dominate away from each line center, and refraction by Na will come to dominate near the center of each line, assuming solar abundances.
For the dynamic polarizability of Na, the oscillator response (Jackson 1975) is first approximated for frequencies near the resonance, and then thermally averaged to account for the Doppler shifts due to the thermal motion of the atoms, with the result ( ) where the second expression is valid on the damping wing.
Here ω 0 = (ω 1 + 2ω 2 )/3 is a weighted center of the Na D doublet, x i = (ν − ν i )/Δ is the frequency in Doppler units relative to line i, and Δ = ν 0 v th /c is the Doppler width, which is the same for each line.The thermal velocity of Na is = and m Na is the mass of a Na atom.The oscillator strengths of the transitions are denoted f 1 and f 2 .The damping parameter is a = Γ/(4πΔ), where Γ is the downward transition rate.The function is the Faddeeva function (Abramowitz & Stegun 1965) 1. Numerical results for N − 1 are shown in Figure 1 for T = 2000 K and n = 1 mbar/k B T = 3.6 × 10 15 cm −3 .There, λ 2 is the D 2 line of the doublet (see Table 1).There is a wavelength interval of width ∼0.1 Å near each line center in which the index of refraction varies significantly in wavelength and can be much larger than the constant H/He value.Since N − 1 can take on both positive and negative values there, refraction due to the Na resonance can cause bending either toward or away from the planet.
Since N − 1 ∝ n, refractive effects are most important at high gas density.However, the trajectories that travel through such high-density regions are subject to large absorption.The absorption coefficient for the gas of H 2 , He, and Na is and a mixture-averaged cross section has been defined as ´-[ The mass of the Na atom (dominant isotope) is m Na = 23 m p , where m p is the mass of the proton.
Figure 2 shows the mixture-averaged cross section over a broad range of wavelengths near the Na D doublet.The Na D opacity dominates Rayleigh scattering over a wavelength interval ∼400 Å wide around the Na D doublet.

Bending Angle and Optical Depth
As the atmospheric refraction will only cause small changes in the photon direction, its effect can be estimated as the kick imparted near the densest point of the trajectory evaluated along a straight line orbit x(t) = b + nct, where b is the impact parameter relative to the center of the planet and n is the photon direction, with n • b = 0.
Wave propagation for a dielectric medium has been discussed by Weinberg (1962).The dispersion relation connecting ω to k can be written as For |N − 1| = 1 and assuming weakly damped waves, the photon trajectories are then and a kick vector, perpendicular to n, of The integration is over the unperturbed trajectory through the atmosphere.
Table 1 Na D Line Properties

Note.
Data taken from Morton (1991).The wavelengths listed are in a vacuum.
For simplicity, an isothermal atmosphere model with constant scale height and mixing ratios will be assumed: where the pressure scale height is H = k B T/μm p g, the gravity is g GM r p 0 2 = , and the density has a value n 0 at the reference radius r 0 .A M p = M J and r 0 = R J planet is assumed.We choose r 0 for the radius at which the pressure P = 1 bar.For an isothermal atmosphere, the polarizability and absorption The fiducial bending angle at the reference radius is and all the position dependence is now in the exponential.For the case of normal dispersion, α > 0, the path is bent toward the planet.However, near the center of each line is a wavelength range over which anomalous dispersion occurs, and the trajectory is bent away from the planet.The images formed in this case differ from the case of normal dispersion.
Although it has been assumed that the waves are weakly damped, i.e., the polarizabilityʼs real component is much larger than the imaginary component,  Re Im a a, this is not the case in the Doppler core of each line.However, the absorption coefficient is so large there that the τ = 1 surface is at low density where refractive effects are small.Hence, we will use the weakly damped assumption even in the line core where it is only a rough approximation.
The optical depth can be similarly evaluated along the straight-line trajectory, with the result where the optical depth at the reference radius is and the mixture averaged cross section is used.The pressure level of the τ = 1 trajectory versus wavelength has been plotted in Figure 3.
The intensity is decreased by a factor e − τ along each trajectory, and hence little flux is received from τ ? 1 trajectories.For a τ(b, λ) = 1 trajectory, the bending angle Θ 1 (λ) can be rewritten as Figures 4 and 5 show the bending angle over a broad and narrow wavelength range, respectively, for the τ = 1 trajectory.The behavior of Θ 1 can be understood using Figure 1 as well as the behavior of the Voigt function.In the Doppler core, e x 2 s µ -rises rapidly, while α ∝ x goes to zero, and so the refraction effect is small.On the damping wing, while Na dominates both α and σ, α ∝ x −1 and σ ∝ x −2 , so that the bending angle increases away from the line center as Θ 1 ∝ x.Further from the line center, α will be dominated by H and He and become constant, while σ is still dominated by Na, and the bending angle will increase away from the line center as Θ 1 ∝ x 2 .Eventually, both α and τ will be dominated by H and He and will vary slowly with wavelength.p s at which τ = 1 along the slant trajectory.The lines show opacity contributions from Rayleigh scattering by H and He, the Na D resonance line, and the sum of the two.The wavelength spacing is 1 Å.
While the above discussion is for the Na D line assuming solar abundances, the conclusion is more general.The bending angle Θ 1 of the τ = 1 trajectory will always decrease toward the line center, and hence the bending angle in the lines is always smaller than the neighboring continuum.For other species with different abundances and lines with different frequencies and oscillator strengths, these parameters only affect the size of the wavelength region over which the bending angle is smaller than the continuum.Even in the case of an abundant species with large oscillator strength, which would dominate the polarizability, Θ 1 would still decrease toward the line center.Hence, the role of the spectral lines is to provide narrow regions of wavelength where the bending angle changes sharply as compared to the neighboring continuum.The isothermal assumption, and fixed mixing ratios with altitude, were chosen for simplicity in order to make analytic progress.A more detailed model would instead compute the atmospheric composition, including chemical reactions as well as photoionization (see, e.g., Huang et al. 2017).In this case, the atmospheric temperature will vary from near the equilibrium temperature near the 1 bar level with a large rise above μbar pressure levels due to thermospheric heating.The photoionization and heating will lead to the dissociation of H 2 into H atoms and eventually, protons, and heavy atoms such as Na will also become progressively more ionized.The effect on the model will be to increase the scale height, which will tend to slow the upward decrease of density, but at the same time, ionization will decrease the density, especially for neutral Na.While the index of refraction of H is comparable to that of H 2 , the ionization of Na will give smaller opacity and push the τ = 1 trajectory lower in the atmosphere.The dependence of the refraction effect on atmospheric properties and star-planet distance has been discussed in detail in Sidis & Sari (2010).As this paper is more concerned with exploring a specific physical effect, wavelength-dependent refraction, these more detailed considerations are beyond the scope of the current work.

Image Flux
Consider the coordinate system shown in Figure 6 in which the observer sits at the origin, the center of the source is at De z , and the position in the source plane relative to the center of the source is x = e x x + e y y.The lens plane is at de z relative to the observer, the center of the lens is displaced from the center of the source by the vector ℓ = e x ℓ, and the impact parameter relative to the lens is b = e x b x + e y b y .Hence, ℓ is the distance of the planet from the center of transit in the sky.The planet-star separation along the z-axis is D − d = D and so the star-planet distance is While the small angle approximation is good for the observer-lens part of the trajectory, it is less so for the source-lens angle.
The unit vector from lens to observer is and that from the source to the lens is Hence, the lens equation, Θ = n o − n s can be written in the form where the approximation d/D ; 1 has been used, the length t 0 ≡ (D − d)Θ 0 and the vector s = ℓ − x have been defined.The solutions of Equation (21) will be considered in Section 5.In the absence of bending, straight line motion from source to observer has impact parameter b = −s so that ℓ + b = ℓ − (ℓ − x) = x.Figure 7 shows an example of the two solutions for the impact parameter as a function of s.
For stellar intensity I, the flux solely due to the star, as measured by the observer, is F 0 ; ∫dΩI, where the integration is over the solid angle dΩ = dxdy/D 2 of the star as seen by the observer.The goal is to understand how this flux is changed by absorption and refraction in the planet.Since radiative intensity is constant in the absence of absorption, the flux in an image is just related to the solid angle on the sky as seen by the observer.Performing the integral over the lens plane, this solid angle involves x y e 1 , , . 2 3 x y The ratio of the lensed and absorbed flux to the background stellar flux is then x y e 1 , , . 2 4 x y with β = 0.63.This "Eddington" limb-darkening law, derived for a gray atmosphere, is a simple and approximate formula used for solar-type stars in the optical (Mihalas 1978).This gives Straightforward, but tedious differentiation gives the Jacobian factor to be and the lower sign is for b s b = -.The integral over the star is computed using cylindrical coordinates, using the Python function dblequad, with the error parameter rtol = atol = 10 −12 .Given r and f position on the star, the trajectory can be solved to find b x (x, y) and b y (x, y).The magnification and optical depth τ can then be computed for that (b x , b y ) so that the integrand is known.
The effect of lensing can be compared to the straight-line case by setting the Jacobian to unity and using the straight-line trajectory b = −s.

Location of the Images
This section discusses the detailed solutions of the lens Equation (21).There are four cases to consider: bending toward (t 0 > 0) or away (t 0 < 0) from the planet, and the image near the unperturbed trajectory, b s b = -, or the image on the opposite side, b s b = .Since the flux in the anomalous dispersion case (t 0 < 0) is heavily absorbed, only t 0 > 0 will be considered here for brevity.
The left-hand side is a monotonically decreasing function of b for which there is a solution for either ±s.
First, consider the solution for the source position directly behind the center of the planet, at s = 0.There is a symmetric ring solution around the planet at impact parameter b = b ring .Equation (26) becomes ( ) ) These trajectories will need fairly large bending angles ∼r 0 /(D − d) ∼ 10 −2 in order to be bent toward the observer, and hence b will have to be smaller than r 0 by many scale heights.An approximate solution is then found by replacing b ring ; r 0 on the right-hand side to find . 28 ring 0 0 0 -These trajectories will have large optical depth and the intensity emerging from the planet will be greatly attenuated by the e − τ factor.Hence, the planet acts as an opaque disk for these trajectories with s  r 0 . For The flux in each image depends on both the Jacobian and also the amount of absorption.Solutions with b  r 0 have both small Jacobian and also heavier damping than solutions at larger b, and hence will have much smaller flux.

Numerical Results
Results are presented for a planet with mass M p = M J , 1 bar radius r 0 = R J , and atmospheric temperature T = 2000 K, orbiting a star of radius R s = R e at distance D − d = 10R e .Constant mixing ratio with solar abundances is used for the atmosphere.
The fluxes F 1 and F 2 correspond to the two roots discussed in Section 4. They are compared to the flux F nb assuming straightline trajectories.The fluxes are normalized by F 0 , the flux of the star alone.The fluxes depend on wavelength λ and the distance of the planet from the center of transit, ℓ, here presented in units of the stellar radius R s .The results are presented for a wide wavelength range 4000-8000 Å with a spacing of 100 Å.
As the focus of this paper is the effect of refractive lightbending, the intrinsic stellar spectrum I is taken to be constant in wavelength, so that deviations due to absorption and refraction may be more easily perceived In order to investigate the effect refraction, the deviation of F 1 + F 2 and F nb is shown in Figure 8.The focus is on the largest deviations, and so only ingress and egress are shown.
The main wavelength dependence is that the deviation becomes very small near the center of the line.This is due to the increase in the Na D cross section and the τ = 1 trajectory moving higher in the atmosphere where the density and bending angle are smaller.Even though the refractive index N − 1 becomes larger near the line center, this effect is overshadowed by the much stronger increase in absorption cross section.
Starting just before the planet begins to transit the star (ℓ/R s = 1.11), the deviation is small.It rapidly increases in size at ℓ/R s = 1.1, and is negative, meaning that lightbending causes a decrease in flux.As more of the planet enters the disk of the star (ℓ/R s = 1.09, 1.08), the deviation becomes smaller and eventually changes sign so that the lensed flux is larger than the straight-line flux when more than half of the planet is inside.The deviation grows for ℓ/R s = 0.93 − 0.91 and reaches a positive maximum for ℓ/R s = 0.9.Beyond that point, where the planet's disk is completely inside that of the star, the deviation is much smaller.As discussed in Hui & Seager (2002), if the planet was entirely contained within the stellar disk and the stellar disk was perfectly  uniform, there would be no lensing signal.It is the deviation of the stellar limb darkening around the disk of the planet that gives rise to a nonzero signal, and it is suppressed due to the planet's small size as compared to the star.
Figure 9 shows the lightcurve, flux deviation versus ℓ, for a range of wavelengths.The peaks at ℓ/R s = 0.9 and 1.1 are apparent, as discussed for Figure 8.The most resonant case, λ = 5900 Å, has the smallest deviation as the τ = 1 trajectory encounters such small densities.The deviation is seen to grow further away from the line center, with λ = 5800 and 6000 having similar sizes, and 5700 and 6100 even larger.Hence, there is a large difference in the deviation over small wavelength ranges.

Conclusions
This paper discusses the effects of refraction and absorption for the Na D resonance line and the neighboring continuum.The index of refraction (N) and absorption cross section are reviewed including contributions that are slowly varying in frequency from the dominant species H 2 and He, and the more rapidly varying resonant contribution from the trace species Na near the Na D doublet.
In the normal dispersion case (N > 1), light is refracted toward the planet, and this is found to be the more important case with the largest bending angles and effect on the transmission spectrum.Anomalous dispersion (N < 1) does occur over part of the Doppler core of the line, and over this wavelength range, the light is bent away from the planet.However, the absorption cross section is so large over this wavelength range that the τ = 1 trajectory passes high in the atmosphere where particle densities are low and the bending angle is small.This is due to the much faster rise of the cross section toward the line center as compared to the index of refraction.As a result, even though the opacity may be dominated by the Na D line, the index of refraction is dominated by the slowly varying static polarizability from H and He.
The thin lens equation is used along with the "kick" to the trajectory, which is applied near the point of closest approach of the photon trajectory to the planet.The location of the images and the flux in each image are found, including the refractive lightbending effect as well as absorption along each trajectory.These results are compared to the standard calculation where photons travel along straight lines, including the effect of absorption.
It is found that the effect of lightbending is very small near the center of the line, and increases rapidly away from the line center so that the effect is confined to where the Na D opacity dominates that of the neighboring continuum.Further, the effects are largest during ingress and egress and specifically near the the start and end of each.
For the specific parameters chosen here for the star, planet, and orbit, the effect of lightbending is small, with fractional size ( ) -R R 1 10ppm p s 2 1 ~Q ~for the flux.Detection of lightbending is challenging as it requires measuring small changes in flux at high spectral resolution during ingress and egress.As discussed in depth by Hui & Seager (2002), the bending angle formula depends on the atmospheric scale height H (see Equation (18) for the τ = 1 trajectory), and may allow H to be constrained by comparison of the models to data.More detailed models than presented in this paper would be required to accurately constrain atmospheric parameters.In comparison to continuum light, for which lightbending is more important for the smaller Rayleigh extinction to the red, the size of lightbending effects is diminished over the narrow wavelength regions of spectral lines such as the Na doublet.

Figure 1 .
Figure 1.Dielectric constant N − 1 of the mixture as a function of wavelength near λ 2 .A temperature T = 2000 K and a pressure P = 1 mbar is used.The dotted line is the frequency-independent contribution from H and He.The dashed line is the resonant contribution from the λ 2 line of the Na D doublet.The thick solid line is the sum of these two contributions.The thin solid line shows N − 1 = 0.

Figure 2 .
Figure2.Optical depth of the mixture as a function of wavelength for gas with temperature T = 2000 K.The dotted line is the slowly varying Rayleigh scattering cross section from H and He.The dashed line is the resonant contribution from both lines of the Na D doublet.The solid line is the sum of these two contributions.The spacing of the points is 1 Å, which greatly limits the cross section near the resonances.

Figure
Figure 3. Pressure ( ) k T r H 2 B 0

Figure 4 .
Figure 4. Bending angle for the τ(b, λ) = 1 trajectory as a function of λ, over a broad range of wavelength.The temperature is T = 2000 K.

Figure 5 .
Figure 5.The same as Figure 4 over a narrow range of wavelength around one line of the Na D doublet.
Approximating d ; D, allowing a limb-darkened intensity I, and including the e − τ factor due to absorption by the lens, gives the formula for the flux to be

Figure 6 .
Figure 6.Geometry of the star-planet-observer system near transit.Top panel: the relationship between x, ℓ, b, and s is depicted for the case b s b = .Bottom panel: the lensing geometry in the z-direction.See the main text for further description.

Figure 7 .
Figure7.The two solutions for the impact parameter as a function of source position s.The lines labeled b 1 and b 2 are the impact parameters of the two images.The line b ring is the s = 0 solution, to which b 1 and b 2 converge as s → 0. Note that the image on the opposite side of the planet as the source asymptotes to b ; s, while the solution on the same side of the planet moves slowly inward of r 0 as s increases.The parameters used were r 0 = R J , H = 0.0039r 0 , and t 0 /H = 16.

Figure 8 .
Figure8.Effect of lightbending on the transit spectrum.The abscissa is wavelength λ in angstroms.The ordinate is the deviation of the flux including lightbending, F 1 + F 2 , and the flux ignoring lightbending, F nb .Each line is labeled by the value of ℓ/R s , the impact parameter of the planet from the center of the stellar disk.

Figure 9 .
Figure9.Effect of lightbending on the lightcurve.The abscissa ℓ/R s is the distance of the planet's center from the center of the star in units of the stellar radius.The ordinate is the deviation of the flux including lightbending, F 1 + F 2 , and the flux ignoring lightbending, F nb .Each line is labeled by the wavelength λ in angstroms.
The mixture-averaged polarizability is α = ∑ i ξ i α i .The static polarizabilities are used for H 2 and He,