Study of atmospheres in the solar system, from stellar occultation or planetary transit

Stellar occultations and transits occur when a planetary body passes in front of a star (including our Sun). For objects with an atmosphere, refraction plays an essential role to explain the drops of flux and the aureoles observed during these events. This can be used to derived key parameters of the atmospheres, such as their density, pressure and temperature profiles, as well as the presence of atmospheric gravity waves and zonal winds. Here we derive from basic principles the equations that rule the ray propagation in planetary atmospheres, and we show how they can be used to derive the physical parameters of these atmospheres.


Introduction
As planetary bodies move in space, they may pass in front of another object, as seen from an observer on Earth or from an instrument on board a spacecraft. Different terminologies are used to describe these phenomena. An occultation occurs when a body blocks the light from a background object. A typical example is a stellar occultation, where a planetary body passes in front of a star. In these cases, the physical disk of the star usually appears as much smaller than the size of the occulting body itself. For instance, the angular diameter of a star projected at the typical distances of an asteroid amounts to a kilometer at most, while the asteroid itself may have a diameter of tens of kilometers.
Conversely, transits occur when the foreground object is angularly much smaller than the background object. Famous examples are transits of Mercury or Venus in front of the Sun as seen from Earth. In the last two decades, transits of exoplanets in front of their stars have been a very powerful tool to discover exoplanets, assess their sizes and orbital periods, measure perturbations from other exoplanets around the same star, or detect chemical species in their atmospheres.
Finally, the term eclipse refers to a body casting its shadow on another body. A well known example is given by lunar eclipses 1 . While eclipses of the Galilean satellites have been observed since their discovery by Galileo, these phenomena among satellites of Jupiter and Saturn have been widely observed In the last decades to pin down their orbital elements and assess secular trends caused for instance by tidal effects.
Here we focus on stellar occultations and transits involving bodies with atmospheres. During a stellar occultation by an opaque object, the star abruptly disappears when reaching the limb of that object. More precisely, the sharpness of the disappearance and reappearance of the star are only limited by the stellar diameter and diffraction effects, and typically last for a fraction of a second only.
In contrast, occultations by objects with an atmosphere are gradual and may last for several minutes. In fact, even if the occultation is diametric, the occulted star may remain faintly visible during the entire event, due to the refraction of the stellar rays by the atmosphere. Contrarily to what is often thought, the gradual character of atmospheric occultations is usually not caused by absorption (due for instance to hazes). In fact, even a completely transparent atmosphere can cause the gradual disappearance of the star (or on the contrary, its brightening), through refraction effects.
In fact, the atmosphere acts as a lens that may focus or defocus the stellar flux. As we will see, the observed phenomena bear some ressemblance with gravitational lenses, where the ray bending stems from gravity instead of gas refraction.

The equations of refraction from basic principles
We recall here a few basic principles that provide the equations of propagation of a luminous ray in an atmosphere. The first principle was stated by Fermat, who noted that during its propagation between two points A and B , a luminous ray minimizes the time of travel between A and B . Since the velocity of light is v = c/n, where c is the speed of light in the vacuum and n is the index of refraction, Fermat's principle may be expressed as the fact that the optical path l , where is stationary (and usually minimal) during the ray propagation. Fermat's principle can be derived from the undulatory nature of light (Huygens' principle), and more generally, from the more modern principle of least action, widely used in Quantum Mechanics 2 . In practice, the problem of minimizing the integral l above is solved by the classical Euler and Lagrange equation, which provides ( [2] p. 1.20-24) with τ = nû, whereû is the unit vector tangent to the ray at r and ∇n is the gradient of n( r ).
In the case where the ray propagates in a plane, the equation above can be written in another way. During an elementary displacement d s, the ray is deflected by an elementary angle d ω, so thatû suffers a deviation d ω given by wherev is the unit vector perpendicular toû. Using Eqs. 2 and 3, we obtain after elementary calculations which will be used in the rest of this chapter. Note that Eq. 4 is equivalent to Eq. 2 only in the case of a planar propagation of the ray. This is not true anymore for a 3D propagation, where the torsion of the ray must also be accounted for.

Refraction by planetary atmospheres
The effects of refraction by a planetary atmosphere during Venus transits have been studied quite long ago, with articles dating back to the eighteenth and nineteenth centuries, see Section 7.
Applications to stellar occultations, on the other hand, started to be discussed one century ago or so, in particular by Anton Pannekoek [3] in 1903 and by Charles Fabry [4] in 1929 .
Observations of stellar occultations were difficult at that time, though, because they required a fast and sensitive photometric recording device. Only on 20 November 1952 was an occultation of the star σ Arietis by Jupiter recorded [5]. Another occultation was monitored on 7 July 1959, involving Venus that passed that time in front of the bright star α Leonis (alias Regulus) [6]. More than one decade was necessary to have another event recorded, on 13 May 1971, when Jupiter occulted the star β Scorpii [7][8][9].
The formal equations describing the effects of atmospheres during occultations are found in various works from the beginning of the 1970's [8][9][10]. They consider the simplest possible case: a planet with a spherically symmetric and transparent atmosphere. Its refractive index n(r ) then depends only on the distance r to the planet's center. A ray coming from infinity from the left with impact parameter p 0 is refracted by the atmosphere and returns to infinity at the right with a total deviation ω(p 0 ), after passing at a closest distance r 0 (p 0 ) from the planet's center (Fig. 1). From the symmetry of the problem, the ray propagates in the plane that passes through the center of symmetry of the atmosphere, as depicted in Fig. 1.

Figure 1.
A ray coming from infinity from the left with impact parameter p 0 is refracted by the atmosphere of a planet. It reaches its closest approach at r 0 and goes to infinity again with a total deviation of ω(p 0 ) that depends only on the incoming impact parameter. Note that by convention here, the angle ω is negative while i is positive.

Bouguer's rule
Let us denote i the angle between the position vector r of the current point along the ray path and the tangent to the path, see Fig. 1. We have whereû is the unit vector in the direction of the ray propagation,r is the unit radial vectorr = r /r , andẑ is the unit vector perpendicular to the plane of Fig. 1. Using τ = nû, we define B as As the ray progresses by a displacement d s along the path, B varies at the rate This directly results from the fact that the atmosphere is spherically symmetric, so that ∇n is parallel to r , and from the fact that τ is parallel toû by definition. We then obtain the Bouguer's rule, d B /d s = 0, i.e. nr sin(i ) = constant along the path.
Very far away from the planet (at left in Fig. 1), there is no atmosphere (n = 1) and r sin(i ) = p 0 . At closest approach, i = π/2 and B = n 0 r 0 , where n 0 = n(r 0 ). Hence which relates the impact parameter p 0 , the closest approach distance r 0 of the ray and the index of refraction n 0 at that point. By integrating d ω/d x from +∞ (corresponding to a ray coming from infinity at the left of Fig. 1) to r 0 (the closest approach to the planet), we obtain half of the total deviation, so that the total deviation is

Total deviation of the ray
where we recall that p 0 = r 0 n 0 . Suppose that we know the refractive structure of the atmosphere, i.e. the profile n(r ). Then, for each impact parameter p 0 , we can determine the closest approach distance r 0 by solving the equation r 0 n(r 0 ) = p 0 . This entirely determines ω(p 0 ) using a numerical scheme to calculate the integral above. In other words, to each impact parameter p 0 , we can associate a deviation angle ω(p 0 ). Figure 3 summarizes the principle of differential refraction in an atmosphere, and defines the various geometrical quantities used in the text. For the commodity of plotting, the deviation is sketched as an abrupt change of propagation, which it is actually gradual (Fig. 1).

Abel inversion
Stellar occultations provide the deviation ω of the ray and the corresponding impact parameter p at the various instants of the event, thus yielding p(ω). The problem is now to derive the value of n(r 0 ) using all these values p(ω). This is done through an Abel inversion of Eq. 6, which is detailed for instance in [8][9][10] and provides which is a relation between n 0 and p 0 = r 0 n 0 , from which n 0 corresponding to the distance r 0 , i.e. n(r 0), is retrieved. The geometry of a refractive occultation. The stellar rays come from infinity at left with impact parameter p. They are differentially refracted as they probe deeper atmospheric layers, and reach the observer located at z after traveling the distance D.
Outside the atmosphere (e.g. at z out ), the deviation ω is zero. Right: A sketch illustrating the relation between the angular diameter θ of a star and the angular diameter of its refracted image, θ . See text for details.
During an occultation by a transparent planetary atmosphere, the flux of the star gradually dims due to the differential deviation of the stellar rays. Thus in this case, the dimming of the flux is not caused by absorption or scattering (due for instance to hazes) or but by refraction.
Consider in Fig. 3 a planet with an atmosphere that deviates a stellar ray with impact parameter p by an angle ω (which is negative by convention, see Figs. 1 and 3), and reaches the observer at z. We have which provides If the atmosphere is transparent, the luminous flux contained in the beam of width d p is retrieved in the beam of width d z. Consequently, the stellar flux at z is "diluted" and yields the irradiance 3 : taking a stellar irradiance outside the occultation normalized to unity. Thus we obtain ω(z) by In practice, it is enough to start the integration just outside the occultation, at some level z out where we have φ = 1 and ω = 0, see Fig. 10. This is reached rapidly, as planetary atmospheres decay exponentially in density, and become undetectable by the observer above a certain level.
Once ω(z) is known, p(z) is given by Eq. 8, which provides ω(p). This is finally introduced in Eq. 7 to obtain the refractivity profile n 0 at radius r 0 .
Here we drop for sake of simplicity the index 0 and we use r instead of r 0 . The refractivity of the gas at a given radius is defined as ν(r ) = n(r )−1. It is related to the molecular density of the gas n g by ν = K n g , where K is the molecular refractivity of the gas under consideration. Consequently, the Abel inversion eventually provides the density profile of the atmosphere through On the other hand, the hydrostatic equation provides the pressure p by integration of the equation d p d r = −µn g (r )g (r ), (12) where µ is the molecular mass and g (r ) the acceleration of gravity at radius r . Finally, the ideal gas equation Both Eqs. 12 and 13 are first order differential equations. As such, they require a boundary condition, i.e. the value of p and T at some prescribed radius r , respectively. This is not too much of a problem for the pressure, as the atmosphere decays exponentially with radius. Thus, we can safely take p = 0 at z out , as we did for ω.
This approximation cannot be used for the temperature, as it is usually not known at z out , and it is certainly not zero. In fact, integrating Eq. 13 requires an independent knowledge of T at some given radius. Otherwise, an infinity of profiles T (r ) can explain the same observable (here, the occultation light curve). This ambiguity can be resolved for instance by using other ground-based observations or spacecraft measurements that have access to the level probed by the occultation. Another approach is to propose physical arguments (such as a radiative transfer model) that restrict the range of plausible values of T at some level.

Conservation of energy, primary and secondary stellar images
We now consider the problem of stellar images during a refractive occultation. For this, we have to reverse the diagram displayed in the left panel of Fig. 3, as shown in the right panel of this same figure. Let us consider an observer at coordinate z (point A) who watches though the atmosphere a star at infinity with angular diameter θ, subtended by the blue and red rays in the figure. These two rays are deflected by the atmosphere and reach the observer at A, where they subtend an angle θ . This angle defines the angular diameter of the stellar image after refraction. We now rotate the red ray around B by an angle θ, A and C will superimpose onto A and C , respectively. In that case, BC is parallel to the outgoing blue ray.
For large distances D, θ = d z/D and θ = d r /D. Using Fig. 3 and p ∼ r , we have from the conservation of energy in a transparent atmosphere This means that the stellar image is compressed by a factor of φ perpendicularly to the limb of the planet. Consequently, the brightness of the stellar image through the atmosphere (i.e. the flux received per unit surface and unit angle at the observer) is the same as the brightness of the stellar image outside the atmosphere. The equation above actually states the conservation of the specific intensity (also called radiance or brightness) of the ray as it propagates through the transparent atmosphere, a theorem due to Clausius (see [11] and the discussion in [2] p. 1.25). Another, equivalent way to state this conservation law is to say that the observed irradiance is proportional to the angular dimension of the stellar image seen through the atmosphere.
For the moment, we have considered that the refraction acts in the plane of the figure. In reality, the limb curvature also causes a deviation of the rays, perpendicular to that limb. It is then easy to extend the result obtained above to a full 2D image, replacing the angle θ (resp. θ ) by the solid angle Ω (resp. Ω ) subtended by the star (resp. its image). Then Eq. 14 can be re-written as which states again that the received flux during an occultation by a transparent atmosphere is directly proportional to the apparent size of the refracted stellar image.

The straight line approximation
For ground-based stellar occultations, the deflection angle ω (Eq. 10) is very small because the distance D is very large. This angle is actually of the order of the angular diameter of the observed body, typically a few arc seconds, i.e. less than 10 −5 radian. From Eq. 5, this implies that r 0 is very close to p 0 , to within p 0 ω, so that n 0 is very close to unity. Consequently, in Eq. 6, we can write n ∼ 1 and p 0 ∼ r 0 . Moreover, d n/d r = d ν/d r , so that One can change the variable r to l = r 2 − r 2 0 (Fig. 4), which leads to where r = r 2 0 + l 2 . Note that the equation above can be used for ray tracing purposes, once a density profile n g (r ) -and thus a refractivity profiles ν(r ), see Eq. 11) -has been prescribed.

The small scale height approximation
In most of the cases, planets have atmospheres with a roughly constant density scale height H , defined as H = − n g d n g /d r .
In usual cases, T varies much more slowly than n g , so that H also varies slowly with r . Moreover, it is usually much smaller than the planet typical radius, i.e. Then where ν 0 = ν(r 0 ). Furthermore, we will see that most of the refractive deviation occurs over a distance ∆l that is significantly smaller than r 0 . Thus, r = r 2 0 + l 2 = r 0 1 + l 2 /l 2 ∼ r 0 (1+l 2 /2r 2 0 ), so that r − r 0 ∼ l 2 /2r 0 and r 0 /r ∼ 1. Using those approximations, introducing in Eq. 16 the expression of d ν/d r obtained above and using This shows that the deviation angle mainly comes from a interval along the ray path of characteristic length For Jupiter, r 0 ∼ 70, 000 km and H ∼ 30 km, so that ∆l ∼ 3, 500 km, which is significantly smaller that r 0 , as announced. For Pluto or Triton, r 0 ∼ 1500 km and H ∼ 20 − 50 km we obtain ∆l ∼ 400 − 700 km. The approximation ∆l r 0 is not so good in those cases, but Eq. 18 still captures the correct orders of magnitude.

The Baum and Code equation
Baum and Code [5] derived a simple equation that describes how the stellar flux decrease when observed from Earth when observing a stellar occultation by a planetary atmosphere. Using Fig. 3, Eq. 9 and p ∼ r (the straight line approximation), we obtain We note that in Eq. 18, the rapidly varying factor in an exponential atmosphere is ν 0 , not r 0 . More precisely, we have d ω/d r ∼ −ω/H , so that, dropping from now on the index 0 in Eq. 18 This shows that the stellar flux has dropped by a factor of two (the "half light level", denoted here by a subscript 1/2) for This means that the ray corresponding to the half light level has been deviated by one scale height H when it arrives at the observer. This occurs for corresponding to a molecular density n g,1/2 = ν 1/2 /K . From Eqs. 8 and 21 and from the definition of the half light level, we have where ∆z = z − z 1/2 . Using the fact that ω/ω 1/2 ∼ ν/ν 1/2 = exp[−(r − r 1/2 )/H ], and from Eq. 22, we have and thus from Eq. 21 This permits to express r − r 1/2 as a function of φ in Eq. 24, and finally get known as the Baum and Code equation [5]. Classical numerical schemes can be used to invert this equation so that to provide φ as a function of ∆z, the distance traveled by the observer in the shadow plane (Fig. 5). Once this is done, Eq. 24 provides the radius of closest approach r probed by the ray as a function of z: Asymptotic expressions of φ can be obtained for ∆z → +∞, i.e. when the star is observed far away from the planet. Then φ approaches unity, so that Eq. 27 yields Consequently, the light curve approaches very rapidly (in fact, exponentially) the unocculted stellar flux unity as the star probes a few scale heights only above the half light radius. To take an example, suppose that the photometric quality of the occultation light curve is such that a drop of flux by at least 1% is necessary to be detected, a common situation in practical cases. This means that we must have exp(−∆z/H ) > 0.01 for detecting the stellar drop. In other words, we can probe only levels below At the other extreme, in the case ∆z → −∞ (i.e. when the star probes deep layers of the atmosphere), Eq. 27 provides Thus, the stellar flux goes to zero rather mildly (∝ 1/|∆z|) when compared to its exponential behavior near φ ∼ 1. The asymptotic behaviors of φ are summarized in the right panel of Fig. 5.

Applications to planetary atmospheres
Eq. 23 shows that the larger the scale height H , the denser the half light level probed in the atmosphere. Conversely, the larger the distance D, the smaller n g,1/2 . This explains why groundbased observations, for which D is very large (up to billions of kilometers), can probe very tenuous pressure levels and still cause significant stellar drops. Using the ideal gas equation and the classical expression H ∼ k B T /µg (r 1/2 ), Eq. 23 provides the expression of the pressure P 1/2 probed by the half light rays: where M is the mass of the body and G is the constant of gravitation. Using this equation and the parameters listed in Table 1, we obtain order-of-magnitude estimations of P 1/2 . Table 1. Estimation of the half light pressure level P 1/2 probed during ground-based stellar occultations. In Table 1, we assume that the objects are observed from Earth near opposition, providing a heliocentric distance D that is roughly the orbital radius of the object minus 1 au. In spite of very large ranges of values for the masses, radii, molecular masses and distances, we see that the combination of these parameters eventually provides a rather narrow range for P 1/2 , typically a fraction of a Pascal, corresponding a few µbar, using 1 Pa = 10 µbar.
The analysis presented in this Section assumes an atmosphere with constant scale height H and adopts approximations that permit the use of the Baum and Code equation. In view of the increasing quality of the occultation light curves and the significant departure of certain density profiles from having a constant H , it is now customary to use "brute force" ray tracing that numerically integrate Eqs. 6 and 13, and then Eqs. 8-10. This allows one to generate synthetic light curves with any atmosphere profiles, and compare them with observations, even in the case of strong local variations of H . Moreover, ray tracing has the advantage to account for the limb curvature, which is not necessarily circular, or not even smooth, as developed in the next Section.
In any case, the stellar occultation method has been very unique and productive for studying planetary atmospheres. Besides the historical cases evoked in Section 3, and without being exhaustive, we may cite the structure and extinction of the Martian upper atmosphere [12], the study of waves in Uranus' stratosphere [13], the discovery of Pluto's atmosphere in the 1980's [14][15][16], its seasonal three-fold pressure increase between 1988 and 2020 [17][18][19], the wave forcing by solar-induced sublimation at Pluto's surface [20,21], the structure and evolution of Neptune's stratosphere [22], the structure, zonal wind regime and haze properties of Titan's stratosphere [23,24] and Triton's atmosphere [25].
An important point is the complementarity between those ground-based observations and space exploration. They may access different regions of the studied atmospheres and thus, provide a synoptic description of these atmospheres. Stellar, solar or radio occultations have been performed by various spacecraft. As they are much closer to the body than terrestrial observers, they probe much deeper layers as the quantity D in Eq. 31 is much smaller.
For instance, while ground-based Titan occultations typically probe a few µbar to some 100 µbar pressure levels, the solar occultations observed by the Cassini spacecraft could reach layers with pressure of more than 10 mbar [26].

Primary and secondary images
From now on, we will consider the curvature of the planetary limb, which creates a "central flash" effect. First, as illustrated in Fig. 6, we see that a spherical planetary atmosphere generally produces two images, a primary image (sometimes called the near-limb image) and a secondary (or far-limb) image. More complex situations where several images are produced by non-spherical atmospheres will be considered later in this chapter. primary image secondary image Figure 6. Stellar rays coming from infinity at left can be refracted to the observed at z following two paths. One produces the primary (or near-limb) image and corresponds to the less refracted ray, plotted here in red. The other ray (in blue) produces a secondary (or far-limb) image that comes from the other side of the planetary disk, and thus suffers a stronger deviation. If z is negative, the primary and secondary characters of the images are swapped.
The figure 7 displays the compression and stretching suffered by the two stellar images in the spherical case, which eventually explains the variation of flux observed during an occultation (Eq. 15). The compression of the image perpendicular to the limb (which decreases the received flux) is due to the differential refraction that defocuses the stellar rays (Fig. 3) by a factor 1 + D(d ω/d r ) (Eq. 20). On the other hand, the stretching of the image along the limb (which increases the received flux) is caused by the limb curvature that focuses the stellar rays toward the shadow center by a factor f (Fig. 7). Thus, the normalized irradiance received by the observer from any of the two (or more) images produced by the limb is  Figure 7. The stellar images caused by refraction during a stellar occultation by an atmosphere. Here, we consider a planet with radius ∼1200 km and an atmosphere with a scale height H ∼ 50 km. The gray circles with a star symbol at their centers delineate the apparent stellar disk projected at the planet distance, here with a radius r = 180 km. The plus symbol marks the planet center, assumed here to be spherical. The occultation proceeds as the star moves from right to left relative to the planet. The value of r has been greatly exaggerated here in order to illustrate the mechanism at play. In real cases, r is at most a few kilometers when considering occultations by solar system objects. At any moment (for instance at points 1, 2 and 3), the stellar disk has two images. One is the primary (or nearlimb, in red) image caused by the refraction due to the nearest point of the limb. The other is the secondary (or far-limb, in blue) image caused by the opposite point of the limb, i.e. rays that pass by the other side of the planet before reaching the observer. The stellar images are compressed perpendicular to the limb due to differential refraction (Eq. 14), and they are stretched parallel to the limb due to the focusing caused by the limb curvature (Eq. 32). The stellar flux caused by a particular image is then proportional to the area encircled in the stellar image. This latter can be calculated using the dashed lines. They show that the stretching greatly increases as the star approaches the planet center, as projected in the sky plane. This leads to the detection of a central flash, see Fig. 8.
The first equation above assumes that the stellar radius r projected at the planet distance 4 is small compared to the atmospheric scale height H , so that the factor d p/d z in Eq. 9 can be considered constant across the star diameter. This is usually the case for planetary occultations, where r is of the order of a few kilometers and H is of the order of a few tens of kilometers. If r is comparable to or larger than H (as it is the case in Fig. 7), then from Clausius' theorem, φ is given by the surface area of the image normalized to πr 2 . No analytical expression of this surface area is available in this case, but it can be easily obtained numerically.
The second equation 32 is restricted to the case of a spherical atmosphere, for which f = r /|z| from the examination of Fig. 7. We have ignored this term so far because it is very close to unity near the half light level. From Eq. 22, z 1/2 = r 1/2 − H . Thus at that level f = r 1/2 /|z 1/2 | = r 1/2 /(r 1/2 − H ) ∼ 1 because usually H r 1/2 . In other words, the focusing effect due to the limb curvature is in general noticeable only near the shadow center. Morever, the second equation 32 is valid only if |z| is larger than a few times r . We will see in the next subsection that the finite stellar size actually prevents the singularity that occurs at |z| = 0 in Eq. 32.
We have also added in Eq. 32 a term exp(−τ) which accounts for the presence of hazes that scatter and absorb the light, where τ denotes the optical depth of the atmosphere along the line of sight. In several cases examined here, the atmosphere is transparent (τ = 0). However, we will give examples on non-zero τ, which leads to a decrease of flux not only because of refraction, but also because of scattering and absorption.
Once the irradiances of each image produced by the limb have been calculated using Eq. 32, they are summed up to derive the synthetic light curve to be compared with the observations. Note that in some cases, the planet is angularly large enough to distinguish the various stellar images moving along the limb. Then, Eq. 32 can be used to produce synthetic light curves for each of these images.

Central flashes: the spherical case
The equation 32 predicts that φ diverges to infinity at z = 0, causing a "central flash". This is true in the limiting case of a point-like source and geometric optics. In actual cases, the star has a finite angular size, so that f (and φ) actually remains finite at the shadow center 5 .
The entire atmosphere can be seen as a lens that focuses the stellar rays toward the shadow center, that can be seen as the focal point of that lens. More precisely, there is one layer (called the flash layer hereafter) in the atmosphere that has the right focal length D, i.e. that causes the right deviation of the ray so that the observer can see the central flash. Thus, the flash can be observed for any value of D, the flash layer being located deeper and deeper as D decreases. Note that the flash ceases to be observed if the flash layer reaches the planet surface, in which case the stellar image vanishes behind the limb.
The figure 8 summarizes the flash process. As long as the apparent stellar disk (delineated in gray in Fig. 8) does not overlap with the planet center (the plus symbol), the primary and secondary images are disconnected. As soon as the stellar disk covers the center, these two images merge into a luminous ring surrounding the planet (panel (f )). This ring actually reveals the flash layer, i.e. the layer that focuses the stellar rays towards the shadow center. This phenomenon is akin to the "Einstein-Chwolson ring" caused by the gravitational lensing of rays coming from a remote galaxy or star by an intervening massive object. So, although the causes of the bending are different (refraction vs. gravitation), atmospheric occultations and gravitational lenses share the same basic process 6 .
For a transparent atmosphere, the brightness of that ring is the same as the brightness of the unocculted star (see the Clausius theorem discussed after Eq. 14). so the flux received at perfect alignement star-planet-observer remains finite, contrarily to what is expected from Eq. 32. From Eq. 14 the width of the Einstein-Chwolson ring is where φ c is the stellar flux at the shadow center without the focusing term f in Eq. 32. Thus, the total area of the ring is (2πr cf )w EC , where r cf is the radius of the flash layer. As the surface area of the star projected at the planet is πr 2 , the maximum height of the flash at shadow center is, normalized to the unocculted stellar flux and from Clausius theorem, 5 However, even for a point-like source, the flux does not diverges at z = 0 due to diffraction effects that are not considered here. 6 For massive planets like Jupiter, the ray bending caused by the gravity of the planet is not negligible compared to the bending caused by refraction. However, its derivative d ω/d r is negligible and thus not affect significantly the stellar flux in Eq. 32. Figure 8. The sequence leading to a central flash. This is the same as Fig. 7, but for a star that goes just behind the planet center (cross). In steps (a)-(c), the primary image is first compressed parallel to the limb, and then stretched parallel to it due to, as illustrated by the dashed lines. As the star aligns with the planet center, the primary and secondary images continue to stretch, until they connect (step (e)) when the stellar limb intersects the planet center (i.e. the cross). At this point, the two images merge into one and form a luminous "Einstein-Chwolson ring" around the body, as illustrated in step (f ).

Einstein-Chwolson ring
where the second equations stems from Eq. 33. For order of magnitude considerations, we can use φ c ∼ H /r cf (Eq. 30), so that The height of the flash decreases as r increases because the flash gets more and more convolved by the stellar disk. Typical values of H are ∼20-50 km depending on the planet, while r is typically of the order of one kilometer. Thus, φ cf may reach values as large as one hundred or more at the very center of the shadow of a spherical and transparent atmosphere. This is indeed the case for Pluto and Triton's atmospheres, as seen later in this section.

Central flashes: the non spherical case
If the atmosphere is not spherical, the equation 32 is still valid, but the factor |z| must be understood as the distance of the observer to the center of curvature of the flash layer. A simple case is when the flash layer assumes a spheroid shape, i.e. an ellipsoid with equatorial and polar radii a and b, respectively (Fig. 9). due to the flattening of the solid planet itself (as it is the case for Mars [12]), or may be maintained by zonal winds, i.e. an atmospheric flow parallel to the equator. These winds create a centrifugal acceleration in a reference frame rotating with the planet. It results in a flattening of the atmosphere under the combined effect of gravity. Elliptical shapes have been used to describe central flashes observed during occultations by Mars [12] or Neptune [27].  Each stellar ray is then deflected perpendicular to the limb of the planet, and converge towards the centers of curvature of the limb (called the evolute) to which the rays are tangent. This creates a caustic where the stellar flux suffers a sudden increase. In the example of Fig. 9, the flash layer appears with an elliptical shape whose equation is The evolute of the ellipse has then the following equation (see e.g. [28]), which is shown in red in Fig. 9. This said, there is no reason why the simple elliptical shape applies in all circumstances. For instance, at its solstice, Titan has weak zonal wind regime in its summer hemisphere (which is then essentially spherical), and a strong jet in the winter hemisphere around the latitude 60 degrees [23,24]. The resulting shape of the atmosphere is delineated in red in the left panel of Fig. 10, with an expansion factor of twenty applied for a better viewing. The resulting intensity map (right panel) is then quite different from the elliptical case shown in Fig. 9.
An observer who is far away from the shadow center receives the flux from the two classical stellar images (primary and secondary) moving in opposite directions (the blue arrows in panel (a) of Fig. 11). The crossing of the caustic causes the sudden appearance of two bright stellar images that moves in opposite directions (red arrows in panel (b)), so that four images are now seen. As the observer proceeds towards the other side of the caustic, the two top-most images approaches each other (blue and red arrows in panel (c)). They coalesce into a bright image at the crossing of the caustic, before disappearing suddenly. As the observer recedes away from the caustic, only the two classical primary and secondary images remain (panel (d)). The right panel of Fig. 11) illustrates how the positions of the four images can be determined at any moment from the shape of the caustic in the case (b). Figure 11. Left: the motion of the stellar images during an occultation by Titan, observed from Gifberg (Republic of South Africa) on 14 November 2008 [24]. These stellar images are reconstructed from the observations of the flash. They could not be seen individually in the data, as Titan was too small (about one arcsec) to be resolved by the instruments used during this campaign. As long as the observer is outside the region delimited by the caustic (see right panel), only the two classical primary and secondary stellar images are detected. When the observer is inside this region, four images contribute to the total flux. They move rapidly along the limb, eventually leading to the coalescence and disappearance of two of them as the observer leaves the caustic domain. Right: A close in view of the shadow center. It shows the position of the observer relative to the caustic, at each steps (a), (b), (c) and (d) illustrated at left. The straight arrows point to the stellar images seen by the observer at step (b). Those straight lines are the four solutions that pass through the point b, while being tangent to the caustic.

Einstein-Chwolson ring and Einstein cross
In the case of stellar occultations by bodies such as Titan, Pluto or Triton, the angular resolution of classical imaging is usually not sufficient to resolve the disk of these objects (that are at the level of one arcsec or less) and thus see the stellar images moving along the limb.  Fig. 11 (left panel). Right: the same where the lens is now a double galaxy (the diffuse objects at the center). Again an Einstein cross is visible, with four images of the quasar 2M1310-1714, but also an Einstein-Chwolson ring, which is the image by the lens of the extended galaxy which hosts the quasar. Credit: ESA/Hubble and NASA.
As mentioned earlier, gravitational lenses act on objects like galaxies that are much more extended angularly than the planetary bodies mentioned above. It is then possible to resolve the images and obtain a direct illustration of the Einstein-Chwolson ring. If the source is a point-like object (like a quasar), it is even possible to see the various images provided by the foreground lens, for instance the four images seen in panels (b) and (c) of Fig. 11, often dumbed as the "Einstein cross". An example of Einstein cross is provided in the left panel of Fig. 12. The right panel displays an image where both the Einstein cross and the Einstein-Chwolson ring are seen.

Effect of atmospheric waves
The left panel of Fig. 13 displays the intensity map of the Titan's flash already shown in Fig. 10. It stems from a flash layer which has a smooth profile. Observations of various occultations by giant planets or Titan, however, reveal irregular structures of the flash. More precisely, rapid  Fig. 10). In the case of Titan, these corrugations are caused by fluctuations induced by atmospheric gravity waves. They cause many streaks in the flash region that result in flux fluctuations, or spikes, in the light curves (Fig. 14).  (Fig. 13). Note that the flash observed in the visible light is much weaker than its counterpart observed in the infrared. This difference stems from absorption by hazes that are more opaque at 0.89 µm than at 2.2 µm [24]. Lower panel: the occultation by Triton observed on 5 October 2017 at Constância in Portugal. This station passed at a mere 6-km distance to Triton's shadow center. The red line is a fit to the data assuming a spherical and transparent atmosphere. The height of the flash represents more than three times the flux of the unocculted star, a current record for this kind of observations. The flux fluctuations seen in the light curve are caused by the Earth atmosphere, not by Triton's atmospheric waves that are much weaker than for Titan. This observation shows that, contrarily to Titan, Triton's atmosphere is essentially spherical and transparent [25].
fluctuations of the stellar flux (or "spikes") are superimposed to the general smooth increase of signal observe during the flash episode, see an example in Fig. 14. These fluctuations are caused by internal gravity waves that propagate in Titan's upper atmosphere. They create small "corrugations" of the central flash layer that break down the smoothing varying centers of curvature of the limb into many centers of curvature. This results into a blurring of the central flash intensity map. In the case of Titan, these corrugations amount to some hundreds meters and cause the blurring illustrated in Fig. 13 [24].

Opacity
In Eq. 32, we mentioned the existence of the factor exp(−τ) which stems from the possible presence of an absorbing material in the atmosphere. In fact, depending on the body, this term An example of haze absorption is given in the upper panel of Fig. 14. A clear difference between the two flashes is observed, due to the differential extinction between the I (visible) and the K (near infrared) bands 7 . More precisely, the chromatic dependence of τ vs. wavelength is such that the atmosphere is essentially transparent in the infrared, while being quite absorbant in the visible.
Since the flash strongly increases the flux when the star is deeply immersed in the atmosphere, it is a useful tool to probe haze properties, and in particular its chromatic dependence. As the zero stellar flux is usually ill-defined (due to the contribution of the occulting body), it is difficult to assess the haze optical depth outside the flash region, where the residual stellar flux is small.
The other example of Fig. 14 is a flash observed during an occultation by Triton. Contrarily to Titan, Triton's flash is completely explained (to within the noise in the data) by a spherical and transparent atmosphere: in that sense, Triton's atmosphere appears as a "perfect lens".

Principle
We now turn to the case where the occulted background star is angularly much larger than the foreground occulting body. As mentioned in the Introduction, this situation is described as a transit (instead of an occultation). This occurs for instance when a exo-planet passes in front of its star, or when the planet Venus is seen transiting in front of the solar disk.
As an example, the geometry of a Venus transit is sketched in Fig. 15, where D and D are the distance of Venus to the Earth and to the Sun, respectively. Without refraction by Venus' atmosphere, an observer would receive a ray from a point S on the Sun that intersects the plane perpendicular to the line of sight passing through Venus' center at ordinate y i . Due to refraction, however, the ray is deflected and appears to come from another ordinate y. In other words, the point S that should be seen at y i has an image that is seen at y.
Taking by convention ω negative (as in Fig. 1) and α and β positive, we have ω = −(α + β). In the limit of small angles, we have α = (y − y i )/D and β = (y − y i )/D, thus On the other hand, we can express ω as a function of y (Eq. 25), so that the equation above can be re-written which defines r 1/2 and where g is the dimensionless geometric factor 8 Eq. 36 provides implicitly y as a function of y i , that, is the position of the image of S as seen in Venus' atmosphere by the observer. This equation can be solved numerically. Alternatively, we can introduce the quantity φ used before (but now as in auxiliary variable) by defining it as (1/φ) − 1 = exp[−(y − r 1/2 )/H ] in analogy to (Eq. 26). From y − y i = y − r 1/2 + r 1/2 − y i , we finally obtain This is similar to the Baum and Code equation 27, except from the appearance of the geometrical factor g and from the fact that the term (1/φ − 2) has been replaced by (1/φ − 1). We will refer to this equation as the "modified Baum and Code equation". For a given y i , the inversion of the modified Baum and Code equation provides φ, which in turn yields y through In this equation, we can choose without loss of generality y > 0. In this case, y i > 0 (resp. y i < 0) corresponds to the primary (resp. secondary) image of S. More generally, this equation can be used to relate the vector position r i of S projected at Venus to the vector position r of its image (Fig. 16). Denoting r i = || r i ||, we can re-write Eq. 38 in a vectorial form. We can encapsulate in the same equations the cases of the primary images and secondary images. This is done by defining a parameter = +1 (resp. = −1) for primary (resp. secondary) images. Then we pose so that the equations relating r and r i are In numerical schemes, one can get rid of the factor in Eqs. 39 and 40 by adopting r i > 0 for primary images and r i < 0 for secondary images. This trick must be used with care, however, as r i classically denotes the modulus of r i and thus, is in principle always positive.
For a given r i , the first equation of the system 40 (i.e. the modified Baum and Code equation) provides u. Once this is done, the second equation provides the vector position r of the image of the point S located at r i . For D → +∞, we have g = 1, and the second equation of the system above reduces to Eq. 28, as expected.
It is instructive to consider the asymptotic behavior of r for ∆r /H approaching −∞, corresponding to φ approaching zero (or equivalently, u approaching +∞), i.e. images of S that are deeply immersed into the atmosphere, referred to as the "deep image regime" hereafter.
The modified Baum and Code equation provides a first estimation u ∼ −∆r /g H by neglecting ln(u) with respect to u. Introducing this expression of u back into the modified Baum and Code equation, we obtain g Hu ∼ −∆r −H ln(−∆r /g H ). Finally, using this approximation of g Hu in the second equation of the system 40, we see that in the deep image regime, the image is located at distance from the planet center.
Because of the weak logarithmic dependence, even for large negative values of ∆r , the deep image of S probes atmospheric layers that are only a few scale heights below the half light radius r 1/2 , At this point, other effects can dominate the aspect of the image. For instance, if the atmosphere is tenuous enough, the image may hit at some point the surface of the planet and merely disappears. In some other cases, the deep atmosphere may be hazy or cloudy, and the image may enter opaque regions, causing also its disappearance. This point is now discussed further in the particular case of Venus transits.  Figure 16. Left: the point S of the solar limb has two images after being refracted by the planet's atmosphere, assumed here to be spherical. The primary (resp. secondary) image is located at r p (resp. r s ). Right: the case shown here is the transit of Venus in front of the Sun, using here a radius of 6130 km for the opaque atmosphere [29]. The quantity h denotes the position of Venus' center (cross) relative the solar limb. By convention, h is positive (resp. negative) if Venus' center is above (resp. below) the solar limb. Two edgings appears along Venus' limb, an upper bright one (the "aureole") and a lower dark one (the "anti-aureole"). They are obtained by using Eqs. 40 to calculate the primary and secondary images of the solar limb. If Venus' center were below the projected solar limb (as it is in the left panel), then the primary and secondary nature of the images would be swapped. Assuming a transparent atmosphere, the conservation of radiance implies that the aureole has the same brightness as the Sun. Likewise, the anti-aureole has the same brightness as the background sky, in black color here. For a better viewing of the aureole and the antiaureole, the scale height H of Venus' atmosphere used here has been largely exagerated (100 km) compared to the actual value of about 4 km [29]. Also, Venus' disk is plotted here in gray for better visibility, but it is by no way luminous in actual observations. It has actually the same brightness as the background sky.

The Lomonossov effect
The discovery of Venus' atmosphere is traditionally attributed to Mikhail Lomonossov, who observed the Venus transit of 6 June 1761 from St Petersburg Observatory [30]. The credit of this original discovery by Lomonossov is debated, though [29,[31][32][33]. However, Lomonossov's basic interpretation was correct, that is, the luminous ring appearing along Venus' limb as the planet emerges from the apparent solar disk is caused by atmospheric refraction. This effect should not be confused with the extension of Venus' crescent near inferior conjunction, caused by the quasi forward-scattering of light by hazes, another evidence that Venus does possess an atmosphere. This point is discussed in the next subsection, and we will see that forward-scattering plays a negligible role during Venus transits, when compared to refraction.
The Lomonossov effect is visualized by calculating the primary and secondary images of the point S as this point is moved along the solar limb (Fig. 16). The time series of Fig. 17 shows in more details the evolution of the two "edgings" resulting from the Lomonossov effect. The bright edging (called the "aureole" for short hereafter) is always outside the solar limb, while the dark edging (called the "anti-aureole") is always inside the solar limb. In practice, the anti-aureole goes unnoticed in the case of Venus, as it is confounded with the planet dark side. However, it must be accounted for when simulating transit light curves involving exo-planets. For instance, the panel (c) of Fig. 17 shows that when the planet center projects itself on the solar limb, the contributions of the aureole and the anti-aureole cancel out. Then, the flux taken away by the planet just corresponds to the flux blocked by its own disk.
The visibilities of the aureole and anti-aureole depend on the value of r opa , the radius of the opaque atmosphere. The image of a given point S on the solar limb is located at a distance r We can consider two cases, depending on the sign of h, the height of Venus' center above the solar limb, see Fig. 18. For h < 0 (left panel of Fig. 18), the aureole is the primary image of the limb, so that = +1, and the condition of visibility of the aureole is where we define the cutoff value r cut . Geometrical considerations based on the examination of Fig. 18 (left panel) show that the aureole is visible if its angle with the vertical is larger than the cutoff angle θ cut = arccos |h| r cut , As r opa decreases, r cut decreases as well, until it reaches the value of |h|. At this point, and for all values of r cut between 0 and h, we have θ cut = 0 and the aureole is uninterrupted along the upper limb of Venus, as illustrated for instance in panel (b) of Fig. 17.
For r cut = 0, the aureole is complete when h = 0 and extends over π radians, a situation illustrated in panel (c) of Fig. 17. A new regime sets in for h > 0, as the aureole now corresponds to secondary images of the solar limb ( = −1), so that Eq. 42 reads This situation is depicted in the right panel of Fig. 18. The aureole now extends over an angle 2θ cut , where again we have θ cut = arccos(|h|/r cut ).
As r opa decreases, r cut increases until it is larger than the planet diameter. Then, the aureole remains visible even if the planet disk is completely detached from the solar disk. This is indeed the case for Venus. During the 8 June 2004 transit, estimates of the various quantities entering the expression of r cut were derived [29]: H ∼ 4 km, r 1/2 ∼ 6170 km and r opa ∼ 6130 km. Using the geometrical factor g = 0.716 relevant to that event, we obtain r cut ∼ 57, 000 km. This means that the aureole could be observed even in Venus' disk is completely detached from the solar disk, a situation illustrated in panel (f ) of Fig. 17.
Beyond the distance of 57,000 km, the aureole disappears behind Venus' opaque atmosphere. As Venus is at ∼ 0.3 au from Earth under these circumstances, this corresponds to an angular separation of a mere 4 arcmin between Venus and the solar limb.
The equation θ cut = arccos(|h|/r cut ) shows that for r cut very large, θ cut approaches π/2, but it cannot goes beyond this value. In other words, the aureole can be observed only along the Venus limb opposite to the Sun, and cannot exceed an angular extension of π radians. This is expected from the fact that solar rays passing near the lower limb of Venus (panel (f ) of Fig. 17) cannot be refracted back to Earth, as the curvature of the ray is always pointing to the direction of increasing refractive index (Eq. 2).
Another aspect, however, limits the observation of the aureole. Even though the brightness of the luminous arc is very high (it is actually the brightness of the solar surface), it becomes very narrow for large values of r i . The maximum width w of the aureole is reached at the top of Venus' limb ( Fig. 18), where r i = h. Thus, the maximum width reached by the aureole is w = r − r opa , where r is given by Eq. 41, with ∆r = −h − r 1/2 , so that where we recall that here h > 0. As an example, we take h = r opa , corresponding to the situation where Venus' disk is tangent to the solar limb (called the first and fourth contacts in the terminology of Venus' transits). Adopting the numerical values mentioned before, we obtain w ∼ 6.5 km. This width cannot be resolved using any imaging technique from Earth. Under usual seeing conditions, the Point Spread Function (PSF) due to our atmosphere of typically one arcsec. This corresponds to about 210 km projected at the planet, which was at 0.288 au from Earth on 8 June 2004. In the example taken above, the PSF dilutes the apparent brightness of the aureole by a factor of 210/6.5 ∼ 30. As Venus' disk recesses away from the solar disk, the width w of te aureole decreases and the dilution factor increases. At some point, the aureole, although intrisincally very bright, becomes too narrow, and thus too faint, to be detected.
In the particular case of Venus' transits observed from Earth, the particular values of r 1/2 , r opa and H are such that w in Eq. 45 is (coincidentally) quite small. Thus, a small variation of r opa induces a large relative variation of w. This creates a patchy aspect of the aureole (Fig. 19) and is a way to map the altitude of the cloud deck along Venus' limb.
The same kind of behaviors occur symmetrically for the anti-aureole. In particular, if Venus' disk is deep inside the solar disk, it disappears behind the opaque layer. Then, it is very difficult by direct imaging to know that the planet has an atmosphere, since the brightness of the Sun, even when observed through the atmosphere, remains unchanged. In this case, the only way to detect the atmosphere is to observe the passage of Venus in front of a sunspot. The image of the spot will be distorted in the same as a stellar disk is distorted when observed through the atmosphere (Fig. 7). This would be an alternative to study Venus' atmosphere. However, such an event has a low probability to occur. Note that another way to detect Venus' atmosphere would be to perform spectroscopic transit observations to detect the gaseous CO 2 .
As a final remark, we note that the general formalism developed here can be applied to exoplanets transiting in front of their stars. Some detailed calculations and applications to exoplanetary atmospheres are given in [34].

Refraction vs. forward-scattering
Near inferior conjunction, Venus' crescent extends beyond 180 degrees, the value expected for an opaque, airless planet. This was first reported by Johann Schröter in 1790, three decades after Lomonossov's observations [35], and then by various observers during the nineteenth century, see the review by Henry Norris Russell in 1899 [36] 9 and the example displayed in Fig. 20.
Russell considered refraction and haze scattering as a possible causes for the extension of Venus' crescent. He then inferred (rightly) that such extension was likely caused by haze scattering rather than refraction. This is in line with the results obtained in the previous subsection: as soon as Venus projects itself at more than a few arcmin from the solar limb, the "refractive aureole" disappears due to the presence of an opaque layers in the atmosphere, leaving only a "haze aureole".
However, during a transit the aureole is largely dominated by refraction, not by haze scattering. Let us consider for this the solar disk with radius R S and apparent surface area S = πR 2 S , emitting rays with radiances (or intensities) I S towards an observer at distance r S , who receives the solar flux over a surface area σ (Fig. 20). Ignoring factors of order unity that account for the angles of emission or reception, the luminous power received from the Sun at σ is Φ S = I S Sσ/r 2 s , so the solar flux f S (also called the solar constant) is f S = Φ S /σ = I S S/r 2 s . Let us denote f V the solar constant at Venus, at distance r V from the Sun. Then, Venus' lit hemisphere re-emits towards σ rays with radiance I V = pΦ(α) f V /π. Here, p is the geometric albedo and Φ(α) is the phase function of the atmosphere, where α is the angle Sun-Venus-observer. Consequently, the ratio I V /I S is As R S ∼ 700, 000 km and r V ∼ 0.7 au ∼ 10 8 km, and because pΦ(α) is of order unity, we obtain I V /I S < 10 −2 1. As mentioned before, I S is also the brightness of the refractive aureole, while I V is the brightness of the lit side of Venus, as seen from the surface element σ. The radiance (or brightness) on the rays forward-scattered by hazes to Earth is itself much less than I V , so that the brightness of the haze aureole is several orders of magnitude fainter than the brightness of the refractive aureole. This makes impossible the detection of the haze aureole, for instance along the limb that is nearest to the Sun when Venus' disk is completely outside the solar disk (panel (f ) of Fig. 17).
More discussion on the respective roles of refraction and scattering in the case of exoplanet transits is provided in [34]. In particular, we note that in exoplanetary cases, the solar radius R S and Venus heliocentric distance r V in Eq. 46 are replaced by the stellar radius R and the planet distance r P , respectively. As R and r P may be comparable, the contribution of forwardscattering relative to refraction may be not overwhelmingly small, and thus should be considered when generating synthetic transit light curves.

Conclusion
The various aspects of refraction during stellar occultations and transits described in this chapter illustrate the long-standing interest of the astronomical community in refraction phenomena. Transits permitted to discover Venus' atmosphere in the eighteenth century, and occultations revealed Pluto's atmosphere in the years 1980's, thus paving the way to the NASA New Horizons flyby of the dwarf planet in 2015.
Meanwhile, occultations and transits still continue to raise great interest. They are used to monitor long-term seasonal evolutions of the atmospheres of Titan, Triton and Pluto, among others. They also probe subtle dynamical effects such as gravity waves that are just impossible to track using any other Earth-based methods.
Occultations and transits have a bright future in store for us. The discovery of tenuous atmospheres around remote Trans-Neptunian Objects will only be possible by observing stellar occultations, to a sensitivity as good as a few nanobars. Also, time series such the one displayed in Fig. 17 can be applied to the study of atmospheres of exoplanets. By comparing transit photometric light-curves with models, one can constrain key parameters such as the scale height, cloud altitude and density profiles of the atmospheres of these remote worlds.