Shadow Imaging of Transiting Objects

The first important degeneracy in the shadow imaging problem results from the reflection symmetry of the star about its horizontal midplane. An opaque shape that transits at an impact parameter b above the midplane produces the same light curve as a "flipped" shape that transits below the midplane.

In planetary transit modeling, this degeneracy can be ignored, because the sign of the impact parameter $b\,=\cos i\left(\tfrac{a}{{R}_{* }}\right)\left(\tfrac{1-{e}^{2}}{1+e\sin \omega }\right)$ is a function of the inclination angle i of the planet's orbital plane and does not describe any inherent property of the transiting planet. However, if we wish to model more general transiting shapes, we must consider the full space of flip-degenerate solutions.

To express the degree of flip degeneracy in a given shadow imaging problem, we consider an image made up of a grid of opaque and transparent pixels, N rows by M columns. (See Figure 1 for an example.) Although there are ${2}^{{NM}}$ unique permutations of opaque and transparent pixels arranged in this grid, each of these permutations does not yield a unique light curve, and in general, a light curve cannot be inverted to produce a unique pixel grid shadow image.

Zoom In Zoom Out Reset image size

We can express the degree of degeneracy by calculating the number of unique light curves, U_LC, possible for this N-row by M-column grid,

$$\begin{eqnarray}\begin{array}{r}{U}_{{\rm{L}}{\rm{C}}}=\left\{\begin{array}{lc}{\left(2\times {3}^{{\textstyle \tfrac{(N-1)}{2}}}\right)}^{M}, & N\,{\rm{o}}{\rm{d}}{\rm{d}}\\ {\left({3}^{{\textstyle \tfrac{N}{2}}}\right)}^{M}, & N\,{\rm{e}}{\rm{v}}{\rm{e}}{\rm{n}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 1 }$$

For intuition, consider first the even-N case. In each of the M columns, there are $\tfrac{N}{2}$ pixels above the midplane; each of these has a counterpart below the midplane with the same impact parameter. There are four possible opacity states for this pair of pixels: 00, 01, 10, or 11. However, because of the flip degeneracy, the 01 and 10 cases produce the same light curve, so only three arrangements are unique—hence, three unique opacity values for each degenerate pixel pair, raised to the power of the total number of pixels above the midplane.

The arithmetic in the odd-N case is the same, except that the middle pixel row has no across-midplane counterpart. Each pixel in that row may only take on an opacity of 0 or 1.

In the case of a square, 3 × 3 pixel grid, there are 2⁹ = 512 unique permutations of opaque and transparent pixels, but only 216 unique light curves.

As a result, the binary-opacity pixel grid solution to any given light curve inversion is not unique, unless the light curve was, in truth, generated by a grid of binary-opacity pixels (τ = 0 or 1) that is symmetrical about its horizontal midplane. Physically, we would only expect such a situation in the case of a perfectly spherical planet, or perhaps a planet–moon system or ringed planet, transiting a star at an impact parameter of 0.

In general, therefore, the inverted pixel grid which generates a light curve is not unique. Figure 1 shows four transiting pixel images which generate identical light curves. Starting with the pixel image at the top and flipping any pixel about the horizontal midplane leaves the light curve unchanged. The pixel image in the bottom panel is the average of the full set of flip-degenerate solutions.

We hope, therefore, to recover shadow images analogous to this bottom panel, which represent a kind of "superposition" of the full set of flip-degenerate solutions to a particular light curve.

There is, however, another degeneracy inherent to the shadow imaging problem by which the set of physically allowable images matching any particular light curve becomes infinitely large. This degeneracy allows a transiting pair of semicircular arcs to generate the same light curve as a single opaque point, and we term it the "arc" degeneracy.

Figure 2 illustrates the geometry of the pair of arcs which generates the same light curve as an infinitesimally small opaque point transiting exactly along the horizontal midplane of the star. Consider this shape to transit from left to right across the star: because the right-hand arc traces the shape of the stellar limb, the entire right-hand arc will ingress upon the star at the same moment, yielding the same vertical ingress feature in the light curve that we would expect from an infinitesimally small transiting planet. (A correspondingly sharp egress feature in the light curve happens when the left-hand arc egresses all at once sometime later.)

Zoom In Zoom Out Reset image size

After the moment of ingress, the top- and bottommost edges of the right-hand arc immediately egress again. However, this egress is balanced by the ingress of the middle of the left-hand arc. If opacity is appropriately distributed along each arc, then the ingress of the left-hand arc and the egress of the right-hand arc may balance exactly. Here, we derive the functional form of the opacity distribution along the arc to allow this exact balance.

Let α (see Figure 2) denote the angle between the horizontal midplane of the star and the point of intersection between the stellar limb and the right-hand arc (which ingresses first). At the moment of ingress, $\alpha =\tfrac{\pi }{2};$ at the moment of egress, α = 0. Let β denote the angle corresponding to the point of intersection on the left-hand arc, and let β range from 0 at ingress to $\tfrac{\pi }{2}$ at egress.

Let θ represent an angle measured from the horizontal midplane of either arc to its outermost point, and let λ(θ) represent the opacity along the arc as a function of this angle. Figure 3 illustrates this setup. Note that λ(θ) cannot be constant, because, for example, during some small time interval dt immediately after the moment of ingress, the length of the right-hand arc which egresses is greater than the length of the left-hand arc which ingresses.

Zoom In Zoom Out Reset image size

Let T be the duration of the transit of the pair of arcs (in other words, the interval between the moment of ingress and the moment of egress). Let the moment of ingress happen at t = 0, and let us define a dimensionless time coordinate $\kappa =\tfrac{t}{T}$ to parameterize the progress of the transit. At ingress then, κ = 0, and at egress, κ = 1.

Following these definitions, we may write

$$\begin{eqnarray}&&\cos \alpha =\kappa ,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\cos \beta =1-\kappa .\end{eqnarray} \tag{ 3 }$$

The total opacity L(κ) transiting the star at a particular moment κ is equal to

$$\begin{eqnarray}&&L(\kappa )={\int }_{0}^{\alpha (\kappa )}\lambda (\theta )d\theta +{\int }_{0}^{\beta (\kappa )}\lambda (\theta )d\theta .\end{eqnarray} \tag{ 4 }$$

For the transit to be flat-bottomed, we require that $L(\kappa )$ be constant, or that $\tfrac{{dL}}{d\kappa }=0$. Differentiating both sides of Equation (4) by κ, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{{dL}}{d\kappa }=\lambda (\alpha )\displaystyle \frac{d\alpha }{d\kappa }+\lambda (\beta )\displaystyle \frac{d\beta }{d\kappa },\end{eqnarray} \tag{ 5 }$$

because λ is time-independent and therefore independent of κ.

Setting this expression equal to 0 and substituting, we obtain

$$\begin{eqnarray}&&0=\lambda (\alpha )\displaystyle \frac{1}{\sqrt{1-{\kappa }^{2}}}+\lambda (\beta )\displaystyle \frac{1}{\sqrt{1-{\left(1-\kappa \right)}^{2}}},\end{eqnarray} \tag{ 6 }$$

or

$$\begin{eqnarray}&&\displaystyle \frac{\lambda (\alpha )}{\lambda (\beta )}=-\displaystyle \frac{\sqrt{1-{\kappa }^{2}}}{\sqrt{1-{\left(1-\kappa \right)}^{2}}}.\end{eqnarray} \tag{ 7 }$$

By the definitions of α and β, we may write

$$\begin{eqnarray}&&\displaystyle \frac{\lambda (\alpha )}{\lambda (\beta )}=-\displaystyle \frac{\sqrt{1-{\cos }^{2}\alpha }}{\sqrt{1-{\cos }^{2}\beta }}=\displaystyle \frac{\sin \alpha }{\sin \beta }.\end{eqnarray} \tag{ 8 }$$

By inspection then,

$$\begin{eqnarray}&&\lambda (\theta )\propto \sin \theta ,\end{eqnarray} \tag{ 9 }$$

whose sign we choose to be positive because physically meaningful opacities are between 0 and 1.

The overall normalization of λ(θ) sets the transit depth of the arcs' light curve. Figure 4 shows a transiting arc pair with $\lambda (\theta )=\sin \theta$.

Zoom In Zoom Out Reset image size

We note that there are two other solutions to λ(θ) that satisfy the condition $\tfrac{{dL}}{d\kappa }=0$. The first is the trivial solution, λ(θ) = 0. The second is a Dirac delta function at θ = 0,

$$\begin{eqnarray}&&\lambda (\theta )\propto \delta (\theta =0),\end{eqnarray} \tag{ 10 }$$

where again the overall normalization sets the transit depth.

For intuition, the two nontrivial solutions to λ(θ) given by Equations (9) and (10) represent two extremes: the least and most compact arrangements of opacity, respectively, that produce the same flat-bottomed, boxlike transit. Any linear combination of these solutions also satisfies $\tfrac{{dL}}{d\kappa }=0$ and generates a boxlike transit.

The above derivation demonstrates that a pair of arcs of variable opacity can match the transit shape of an infinitesimal point of opacity transiting along the horizontal midplane of the star, at impact parameter b = 0. The same logic applies to an infinitesimal point at an arbitrary impact parameter b. Figure 5 illustrates the geometry of this situation.

Zoom In Zoom Out Reset image size

Mathematically, a change in the impact parameter b means that the limits of integration in Equation (4) change,

$$\begin{eqnarray}&&L(\kappa )={\int }_{\arcsin b}^{\alpha (\kappa )}\lambda (\theta )d\theta +{\int }_{\arcsin b}^{\beta (\kappa )}\lambda (\theta )d\theta .\end{eqnarray} \tag{ 11 }$$

Since b is constant, the subsequent steps and resulting solutions for λ(θ) do not change, except that the delta function solution is localized at $\theta =\arcsin \,b$.

We note finally that the arc degeneracy technically only operates for an occulter transiting a uniformly bright star: if the star is limb-darkened, then there is no (unchanging) arc arrangement which can maintain the perfect opacity ingress–egress balance described by Equation (6). However, in practice, the limited time resolution of light curve observations leaves room for significant arc-degenerate behavior in shadow images recovered from real transit data (see Section 5 below).

A third degeneracy inherent to light curve imaging results from the "scale-free" nature of the problem and allows a wide image moving at high velocity to generate the same light curve, within an arbitrarily small measurement uncertainty, as a narrower image moving at lower velocity. We term this degeneracy the "stretch" degeneracy.

The stretch degeneracy is mathematically simpler than the arc and flip degeneracies. Two occulters with the same transit duration T both obey

$$\begin{eqnarray}&&T=\displaystyle \frac{W}{v},\end{eqnarray} \tag{ 12 }$$

where W is the width of the occulter and v is its velocity. The right-hand side of this equation can be multiplied by the same constant in the numerator and denominator without consequence to T. In other words, a "stretched" image traveling fast can generate a transit event of the same duration as a narrow image traveling slowly.

Figure 6 illustrates the stretch degeneracy for a simple, low-resolution circular occulter. Note in particular two features of the "stretched" image: first, that it is semiopaque rather than fully opaque like the unstretched image, and second, that its edges are less opaque than its middle. The semiopacity of the stretched image is necessary in order to match the transit depth of the unstretched image: because the stretched image is wider, it occults more of the stellar surface, so it must let some light through or it will produce a much deeper transit than the unstretched image. Meanwhile, the lightened edges of the stretched image are necessary to better match the ingress and egress shape of the unstretched image's light curve; with arbitrarily high image resolution, it is possible to match the unstretched image's light curve to arbitrary accuracy.

Zoom In Zoom Out Reset image size

In practice, the stretch degeneracy is the least important of the three nontrivial degeneracies we explore in this section, because a fast-transiting, stretch-degenerate image can only match a narrow, slower image's light curve if the image resolution is high enough, as suggested by the example in Figure 6. For real data, image resolution is constrained by the number of observed data points, which causes us to prefer the narrowest, slowest possible image which can match an observed light curve (see Section 4.2 for further discussion).

Finally, we note two trivial degeneracies which do not affect the inference of a shadow image. The first relates to the arbitrary sign of the velocity of the transiter; an image which transits left to right across the star generates the same light curve as the same image, horizontally mirrored, transiting right to left across the star at the same velocity. We choose positive v to indicate that the image transits left to right (see Section 3.2 below).

The second trivial degeneracy relates to a time translation of the entire transit event. As we discuss in Section 3.2, we must choose a "reference time," analogous to a transit midpoint time, along a light curve in order to recover a shadow image; shifting this reference time forward or backward along the light curve results in a shadow image which is shifted right or left, respectively (given our choice of v direction above).

Pixelating, or discretizing, the shadow image is motivated by recognizing that real light curves are themselves discrete time series. A light curve is not infinitely resolved in time, and therefore we should not attempt to recover a shadow image that is infinitely resolved spatially. Similarly, each flux measurement in a light curve has an associated uncertainty; we should not attempt to recover a shadow image with pixel elements too small to be definitively detected within that uncertainty (see Section 4.2 below for further discussion).

Furthermore, discretizing the pixel image enables us to investigate two physical variants of the shadow imaging problem:

• 1.  
  What if the transiting object which generated the light curve is a solid body and therefore our shadow image should only admit of completely transparent (opacity τ = 0) or completely opaque (τ = 1) pixels?
• 2.  
  What if the transiting object is dusty or translucent or is a solid body smaller than the pixel scale and our shadow image can contain pixels of intermediate opacity ($0\leqslant \tau \leqslant 1$)?

These two variations of the shadow imaging problem have different constraints on the pixel opacities and require different mathematical approaches to inversion. In case (1), discretizing the pixel image is necessary to divide it up into opaque and transparent elements. In case (2), discretizing the pixel image enables us to set up the light curve inversion problem as a single-matrix equation (Equation (35) below) and to explore both analytic and numerical approaches to solving this equation. (Similar mathematical formulations exist for both starspot inversion (Vogt et al. 1987) and eclipse mapping, e.g., Berdyugina & Kuhn 2017.)

The same forward model, or procedure for generating a light curve from a pixelated image, can be used in both cases, so we begin there. How do we calculate the light curve of a pixelated image grid transiting a star?

We consider a pixel grid of N rows and M columns transiting a star. We normalize the physical scale of the problem such that the radius of the star is unity.

The grid lives in the X–Y sky-projected plane, with the observer at $Z=+\infty$. The grid moves laterally along the X axis, with ${dX}/{dt}\gt 0$, and does not translate up or down (i.e., ${dY}/{dt}=0$). We illustrate this setup in Figure 7.

Zoom In Zoom Out Reset image size

We treat the grid as moving at a constant lateral velocity $v\equiv {dX}/{dt}$, where dv/dt ≡ 0. This is a reasonable approximation over the timescale of a transit, unless the object resides in a very tight orbit or the object is near the pericenter in a highly eccentric orbit. We define positive v to mean that the grid transits from left to right across the star, such that the rightmost column of pixels ingresses first.

We define the vertical position of the grid such that the top of the highest row of pixels falls at Y = 1 and the bottom of the lowest row of pixels falls at Y = −1. In this way, the grid perfectly overlaps with the star in the vertical dimension.

This definition sets the size of each pixel to have a width w of

$$\begin{eqnarray}&&w=2/N.\end{eqnarray} \tag{ 13 }$$

We emphasize that every pixel has the same square shape with this dimension. For N = 1 then, w = 2 and is thus equal to the diameter of the star.

To refer to individual pixels, we adopt the index notation $i\in [1,N]$ to denote the row and $j\in [1,M]$ to denote the column. To calculate the amount of stellar flux the grid blocks at each discrete time step t_k of the transit observation, we must first calculate the X and Y positions of each grid pixel i, j at each time step.

We may write the Y position of the center of pixel i, j as

$$\begin{eqnarray}&&{Y}_{i,j}=1-(w/2)-(i-1)w,\end{eqnarray} \tag{ 14 }$$

where setting i = 1 recovers ${Y}_{1,j}=1-(w/2)$ and setting i = N recovers ${Y}_{N,j}=-1+(w/2)$. The Y positions of the grid pixels are constant.

For the X positions of the pixels, which evolve in time, we first define a reference X position for each pixel at a reference time $t={t}_{\mathrm{ref}}$ as

$$\begin{eqnarray}&&{X}_{i,j}^{\mathrm{ref}}\equiv {X}_{i,j}[t={t}_{\mathrm{ref}}].\end{eqnarray} \tag{ 15 }$$

We choose the reference time such that the grid is centered on the star at t = t_ref. Therefore:

$$\begin{eqnarray}&&{X}_{i,j}^{\mathrm{ref}}=(j-{j}_{\mathrm{mid}})w,\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&{j}_{\mathrm{mid}}=1+(M-1)/2.\end{eqnarray} \tag{ 17 }$$

We may now write the time-evolving X position of the center of any pixel as

$$\begin{eqnarray}&&{X}_{i,j}({t}_{k})={X}_{i,j}^{\mathrm{ref}}+({t}_{k}-{t}_{\mathrm{ref}})v,\end{eqnarray} \tag{ 18 }$$

where t_k is the kth time index and $k\in [1,K]$. Practically speaking, t_ref is analogous to the transit midtime fitted in conventional transit models.

We may use the above equation for the time-evolving X_{i, j} to solve for the time at which the grid makes first and last contact with the star, t_enter and t_exit. The grid moves from left to right across the star, so at first contact, the Mth column of pixels has an X position equal to −1 − (w/2), and for the last contact the first column of pixels has an X position equal to 1 + (w/2), giving

$$\begin{eqnarray}&&{t}_{\mathrm{enter}}={t}_{\mathrm{ref}}-\displaystyle \frac{1+{Mw}/2}{v},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{exit}}={t}_{\mathrm{ref}}+\displaystyle \frac{1+{Mw}/2}{v}.\end{eqnarray} \tag{ 20 }$$

We assign a time-independent opacity τ_{i, j} to each pixel. τ_{i, j} is a binary value equal to zero or unity—in other words, we construct a grid of perfectly transparent pixels (τ_{i, j} = 0) and opaque pixels (τ_{i, j} = 1).

In total then, our model has ${MN}$ opacity parameters, which are binary-valued (case 1) or real numbers between 0 and 1, inclusive (case 2), and two auxiliary parameters, t_ref and v, which are real-valued.

As pixel i, j transits the star, it occludes a fractional area ${A}_{i,j}({t}_{k})$ of the stellar disk at time t_k; A = 0 for pixels which do not overlap the star, and $A=\tfrac{{w}^{2}}{\pi }$ for pixels which overlap completely (because we choose R = 1, the area of the entire stellar disk is equal to π).

If we assume that the stellar disk is uniformly bright (i.e., there is no limb darkening), we may then compute the light curve ${\boldsymbol{F}}(t)$ of the transiting grid by recognizing that the fractional flux blocked by the grid at each time step t_k is equal to the fractional area of the star occulted by nontransparent pixels (τ_{i, j} > 0), in proportion to their opacity. Therefore, the unocculted flux at time t_k is given by:

$$\begin{eqnarray}&&F({t}_{k})=1-\displaystyle \sum _{i=1}^{N}\displaystyle \sum _{j=1}^{M}{\tau }_{i,j}{A}_{i,j}({t}_{k}).\end{eqnarray} \tag{ 21 }$$

This is the equation for the transit light curve, normalized such that F = 1 out of transit.

We emphasize that while opacities τ < 0 and τ > 1 are mathematically permissible in this equation, they are unphysical: τ < 0 would represent a transiting pixel brighter than the stellar surface it occulted, and τ > 0 would describe a pixel that blocked more than its proper area's worth of starlight.

We compute the area of overlap ${A}_{i,j}({t}_{k})$ of pixel i, j at time step t_k from the (X, Y) position of the pixel's center at t_k, given by Equations (14) and (18), and the pixel's width, given by Equation (13). When a pixel partially overlaps the star, we approximate its overlap area as either a triangle, a trapezoid, or a square missing a triangular corner. We choose the appropriate overlap shape by computing the number of intersection points between the edge of the star and the sides of the pixel and also noting whether the center of the pixel falls inside or outside the star. We then correct this approximation by using the length of the chord between intersection points to calculate the area of the sliver of occluded star yet unaccounted for by this approximation.

For a limb-darkened star, we must also account for the position of each opaque pixel relative to the stellar limb at each time step in order to determine how much flux it occludes. We adopt the small-planet approximation of Mandel & Agol (2002), in which it is assumed that the star's surface brightness is constant across a pixel. In other words, we treat the pixel as occulting a thin, uniform-surface-brightness, annular slice of the stellar disk, where the radius of the annulus is the distance from the center of the stellar disk to the center of the pixel and the annulus is just wide enough to encompass the pixel.

We denote the area of this annulus as A_annulus and its emitted flux as F_annulus. As a rule of thumb, this small-planet approximation is only appropriate for w ≲ 0.2 (i.e., N > 10), which corresponds roughly to an occulter-to-star ratio-of-radii of 0.1, for a circular occulter of the same area as the pixel. The exact ratio-of-radii at which the small-planet approximation becomes inappropriate depends on the impact parameter of the pixel, the size of the pixel, the limb darkening profile of the star, and the bandpass of the observations, so there is no general exact cutoff.

To calculate the light curve in the limb-darkened case, we must renormalize Equation (21): the fractional flux occulted by an opaque pixel is no longer equal to the fractional area of the stellar disk occluded by the pixel, but rather to:

$$\begin{eqnarray}&&{\bar{A}}_{i,j}({t}_{k})=\displaystyle \frac{{A}_{i,j}({t}_{k})}{{A}_{\mathrm{annulus}}}\displaystyle \frac{{F}_{\mathrm{annulus}}}{{F}_{\star }},\end{eqnarray} \tag{ 22 }$$

where ${F}_{\star }$ is equal to the flux of the entire unocculted limb-darkened star relative to the non–limb-darkened star (which must be calculated given a choice of limb darkening coefficients). We note that this equation reduces to ${\bar{A}}_{i,j}({t}_{k})={A}_{i,j}({t}_{k})$ in the case of uniform limb darkening.

The value of the light curve at t_k is then given by:

$$\begin{eqnarray}&&F({t}_{k})=1-\displaystyle \sum _{i=1}^{N}\displaystyle \sum _{j=1}^{M}{\tau }_{i,j}{\bar{A}}_{i,j}({t}_{k}).\end{eqnarray} \tag{ 23 }$$

We provide Python code to calculate the transit light curve of any grid in the case of uniform, linear, quadratic, or four-parameter nonlinear limb darkening in the software package accompanying this paper, EightBitTransit.

To be exact, we note that ${\boldsymbol{F}}$ is a column vector of length K, of which each scalar entry ${F}_{k}\equiv F({t}_{k})$ is given by Equation (23). Let us "unravel" the double sum in Equation (23) by defining a new index l, such that

$$\begin{eqnarray}&&l[i,j]=j+(i-1)M.\end{eqnarray} \tag{ 24 }$$

Given that i ranges from 1 to N, and j from 1 to M, l ranges from 1 to MN.

We may then rewrite Equation (23) as

$$\begin{eqnarray}&&{F}_{k}=1-\displaystyle \sum _{l=1}^{L}{\tau }_{l}{\bar{A}}_{l}({t}_{k}),\end{eqnarray} \tag{ 25 }$$

where we define L ≡ MN. If we further define ${\bar{A}}_{k,l}\equiv {\bar{A}}_{l}({t}_{k})$, then

$$\begin{eqnarray}&&{F}_{k}=1-\displaystyle \sum _{l=1}^{L}{\tau }_{l}{\bar{A}}_{k,l}.\end{eqnarray} \tag{ 26 }$$

Let us now rewrite ${\bar{A}}_{k,l}$ in matrix form:

$$\begin{eqnarray}\bar{{\mathsf{A}}}=\left(\begin{array}{cccc}{a}_{\mathrm{1,1}} & {a}_{\mathrm{1,2}} & \cdots & {a}_{1,L}\\ {a}_{\mathrm{2,1}} & {a}_{\mathrm{2,2}} & \cdots & {a}_{2,L}\\ \vdots & \vdots & \ddots & \vdots \\ {a}_{K,1} & {a}_{K,2} & \cdots & {a}_{K,L}\end{array}\right).\end{eqnarray} \tag{ 27 }$$

$\bar{{\mathsf{A}}}$ is a matrix of shape K by L, where the kth row encodes the state of overlap of the entire pixel grid at time step k and the lth column encodes the overlap state of pixel l across all time steps.

Similarly, we may "unravel" the opacity matrix τ into a column vector ${\boldsymbol{\tau }}$ of length L:

$$\begin{eqnarray}&&{\boldsymbol{\tau }}=\left(\begin{array}{c}{\tau }_{1}\\ {\tau }_{2}\\ \vdots \\ {\tau }_{L}\end{array}\right).\end{eqnarray} \tag{ 28 }$$

We may now re-express Equation (23) in matrix form:

$$\begin{eqnarray}&&{\boldsymbol{F}}=1-\bar{{\mathsf{A}}}{\boldsymbol{\tau }},\end{eqnarray} \tag{ 29 }$$

where

$$\begin{eqnarray}&&{\boldsymbol{F}}=\left(\begin{array}{c}{F}_{1}\\ {F}_{2}\\ \vdots \\ {F}_{k}\end{array}\right)\end{eqnarray} \tag{ 30 }$$

and 1 is a column vector of ones, equal in length to ${\boldsymbol{F}}$.

If we define a vector ${\boldsymbol{R}}=1-{\boldsymbol{F}}$, we may rearrange this equation to read

$$\begin{eqnarray}&&\bar{{\mathsf{A}}}{\boldsymbol{\tau }}={\boldsymbol{R}}.\end{eqnarray} \tag{ 31 }$$

If $\bar{{\mathsf{A}}}$ were invertible, then our work would be done: we could solve Equation (31) directly for the vector of pixel opacities ${\boldsymbol{\tau }}$. However, because of the flip degeneracy, pixel i, j has the same area of overlap at every time step as pixel $(N+1-i),j$, and as a result, $\bar{{\mathsf{A}}}$ always has repeated columns. By the invertible matrix theorem, a matrix with repeated columns is not invertible.

We may proceed by recognizing that $\bar{{\mathsf{A}}}$ and ${\boldsymbol{\tau }}$, given that they describe the entire pixel grid, contain redundant information. We need only solve for the opacities of one half of the pixels (we choose the top half, for convenience). We define a new index

$$\begin{eqnarray}{L}_{\mathrm{half}}\left\{\begin{array}{ll}\tfrac{(N-1)M}{2}+M, & N\,\mathrm{odd}\\ \tfrac{{NM}}{2}, & N\,\mathrm{even}.\end{array}\right.\end{eqnarray} \tag{ 32 }$$

We define a new area-of-overlap matrix ${\bar{{\mathsf{A}}}}_{\mathrm{half}}$, which represents the left half (columns 1 through L_half, inclusive) of $\bar{{\mathsf{A}}}$, and a new opacity vector ${{\boldsymbol{\tau }}}_{\mathrm{half}}$, which represents the corresponding top half of ${\boldsymbol{\tau }}$. We may then write:

$$\begin{eqnarray}&&{\bar{{\mathsf{A}}}}_{\mathrm{half}}{{\boldsymbol{\tau }}}_{\mathrm{half}}={\boldsymbol{R}}.\end{eqnarray} \tag{ 33 }$$

Given that, in general, $K\ne {L}_{\mathrm{half}}$, we may multiply both sides of this equation by ${\bar{{\mathsf{A}}}}_{\mathrm{half}}^{{\rm{T}}}$ to yield

$$\begin{eqnarray}&&{\bar{{\mathsf{A}}}}_{\mathrm{half}}^{{\rm{T}}}{\bar{{\mathsf{A}}}}_{\mathrm{half}}{{\boldsymbol{\tau }}}_{\mathrm{half}}={\bar{{\mathsf{A}}}}_{\mathrm{half}}^{{\rm{T}}}{\boldsymbol{R}}\end{eqnarray} \tag{ 34 }$$

so that both sides of the equation are column vectors of length L_half and ${\bar{{\mathsf{A}}}}_{\mathrm{half}}^{{\rm{T}}}{\bar{{\mathsf{A}}}}_{\mathrm{half}}$ is a square matrix.

For notational simplicity, let ${\mathsf{B}}\equiv {\bar{{\mathsf{A}}}}_{\mathrm{half}}^{{\rm{T}}}{\bar{{\mathsf{A}}}}_{\mathrm{half}}$, and let ${\boldsymbol{C}}\,\equiv {\bar{{\mathsf{A}}}}_{\mathrm{half}}^{{\rm{T}}}{\boldsymbol{R}}$, such that

$$\begin{eqnarray}&&{\mathsf{B}}{{\boldsymbol{\tau }}}_{\mathrm{half}}={\boldsymbol{C}}.\end{eqnarray} \tag{ 35 }$$

We have therefore reduced our shadow imaging problem to the problem of solving a system of linear equations for the entries of the column vector ${{\boldsymbol{\tau }}}_{\mathrm{half}}$. These entries, reshaped into the matrix τ, correspond to the opacities of the pixels making up the top half of the grid, which define the image.

In the sections below, we elaborate on the two steps of shadow imaging: first, selecting the grid parameters and second, solving Equation (35) for the pixel opacities subject to our chosen physical constraints.

In general, the auxiliary parameters t_ref and v can be set to reasonable approximations of their "true" values, and the pixel image will slightly shift and stretch, respectively, relative to the "truth," without disturbance to its principal morphology. This means we can proceed by fixing these terms and optimizing the opacities ${\boldsymbol{\tau }}$ only. We may then, depending on the success of the solution ${\boldsymbol{\tau }}$, perform further iterations, varying the grid parameters each time, to reach an optimal grid with optimal auxiliary parameters. We discuss here some constraints of the grid parameters which allow us to estimate their values initially.

The first constraint we consider is that the number of pixel elements should not exceed the number of data points obtained during the transit event of interest. For regularly sampled data, such as those of Kepler, we may write the sampling constraint as

$$\begin{eqnarray}&&{NM}\leqslant \displaystyle \frac{{t}_{\mathrm{event}}}{{t}_{\mathrm{cadence}}},\end{eqnarray} \tag{ 36 }$$

where t_event is the timescale of the event we wish to image and t_cadence is the cadence of the time series, i.e., the interval between successive observations.

The second constraint we consider is that a pixel should not be too small to detect individually. In other words, the transit depth of a single opaque pixel should not be smaller than the uncertainties on the flux measurements. In principle, smaller pixels could be resolved over repeated transit observations, but this approximation again aids in selecting a unique initial grid size from which to begin optimizing the grid opacities.

Mathematically, we can express the precision constraint as:

$$\begin{eqnarray}&&\displaystyle \frac{{w}^{2}}{\pi }\gtrsim \sigma \end{eqnarray} \tag{ 37 }$$

where σ is the typical photometric uncertainty. Combining Equation (13) with this constraint gives

$$\begin{eqnarray}&&N\lesssim \sqrt{\displaystyle \frac{4}{\pi \sigma }}.\end{eqnarray} \tag{ 38 }$$

For reference, using a 60 ppm uncertainty, this yields N ≲ 146. (In practice, we are usually limited to much smaller values of N by the number of data points in the observed light curve.)

The third constraint we consider is the size of M. Our grid must be wide enough to create a total duration sufficient to explain the event timescale t_event. We require that ${t}_{\mathrm{exit}}\,-{t}_{\mathrm{enter}}\geqslant {t}_{\mathrm{event}}$, or

$$\begin{eqnarray}&&\displaystyle \frac{2+2(M/N)}{v}\geqslant {t}_{\mathrm{event}}.\end{eqnarray} \tag{ 39 }$$

Similarly, we consider that a single pixel needs to be able to traverse the entire disk of the star within the event timescale. The actual duration of a single pixel's transit will depend on the pixel's latitude Y, but to simplify things, we consider an equatorial pixel of infinitesimal size and use an approximate symbol for the inequality, to give

$$\begin{eqnarray}&&v\gtrsim 2/{t}_{\mathrm{event}}.\end{eqnarray} \tag{ 40 }$$

Together, these expressions constrain the velocity to the range

$$\begin{eqnarray}&&2/{t}_{\mathrm{event}}\lesssim v\leqslant 4/{t}_{\mathrm{event}}.\end{eqnarray} \tag{ 41 }$$

As a general strategy then, we choose a grid velocity v equal to 2/t_event, and t_ref to correspond to the minimum of the observed light curve. To choose N and M, we recognize that, for a chosen N, we may solve for M such that the grid continuously overlaps the star, by rearrangement of Equation (20). We can then adjust N to accommodate the constraint that NM should be less than the number of observed data points. Once the grid dimensions have been chosen, we re-execute the inversion for different velocities until the fit ceases to improve.

Because of the resolution constraint, we prefer the slowest grid velocity v which returns a reasonable fit to the observed light curve, because this slow velocity corresponds to the highest image resolution N. This is, in a sense, an image prior which prefers narrow, slow images to their fast, stretch-degenerate counterparts.

Once we have reasonable first estimates for ${t}_{\mathrm{ref}}$, v, N, and M and have chosen a limb darkening law to describe the stellar disk, we may use Equation (22) to solve for ${\bar{A}}_{i,j}({t}_{k})$ for each grid pixel at each light-curve time step. At this stage of shadow imaging, it is helpful to think of the grid pixels as containers for as-yet-to-be-determined opacity: each transits the star in a definite way according to the grid parameters, so ${\bar{A}}_{i,j}({t}_{k})$ and hence $\bar{{\mathsf{A}}}$ are well defined, but its opacity is not yet known.

For a chosen t_ref, v, N, and M, we restrict our attention to the observed light-curve data points that satisfy t_enter < t < t_exit. In other words, we consider only the points in time during which the grid partially overlaps the star, because the transiting grid cannot influence points outside this range.

To determine the opacities, we must solve Equation (35) for the entries of the opacity vector ${\boldsymbol{\tau }}$. Given that this matrix equation is linear, in principle it can be directly, analytically solved.

However, direct solution of Equation (35) cannot accommodate constraints on the pixel opacities. Namely, there is no way to restrict the entries of ${\boldsymbol{\tau }}$ to the physically meaningful range [0,1] (case (1)), let alone to the binary values (0 or 1; case (2)). Mathematically, introducing these constraints transforms the problem into a nonlinear optimization problem, which is not susceptible to solution by a linear matrix equation. We furthermore find that transforming the opacity variables through a logistic function, which maps the real numbers to the range [0,1], results in numerical instabilities in our attempts to solve Equation (35) both directly and iteratively (e.g., with the Simultaneous Algebraic Reconstruction Technique [SART]; see Section 4.3.2 below).

Furthermore, we find that choosing grid parameters t_ref, v, N, and M that deviate even slightly from the true values leads to a completely nonsensical recovered ${\boldsymbol{\tau }}$. Direct analytic solution is therefore not robust enough to apply to a light curve of unknown origin, where our initial guesses for the grid parameters are unlikely to be so accurate.

We therefore explore less exact but significantly more robust algorithmic approaches to solving for ${\boldsymbol{\tau }}$. Below, we discuss each of these algorithms in turn. The first two address case (1), where pixels may take on intermediate opacities, and the latter three address the more restrictive case (2), where pixels are constrained to be binary-valued.

In Figures 8 and 9, we compare their performance in recovering a number of known test grids from noiseless light curves. In these recovery tests, the parameters N, M, v, and t_ref are assumed to be known. Eight of the test grids are binary-valued, and three (the low-resolution planet–moon, planet-ring, and comet test grids) include intermediate-opacity pixels.

Zoom In Zoom Out Reset image size

In Figures 8 and 9, we have chosen to generate our test light curves with a uniformly bright star, i.e., without limb darkening. We make this choice because non–limb-darkened light curves are sharper and less rounded than limb-darkened light curves, and the inversions result in correspondingly sharper image grids, among which the differences between the images generated by our four recovery algorithms stand out most clearly.

Zoom In Zoom Out Reset image size

We find that introducing realistic limb darkening results in recovered images very similar to those shown in Figures 8 and 9, with two notable qualitative differences: First, for the limb-darkened case, opacity tends to be pushed farther out toward the top and bottom edges of the recovered image. This effect is most obvious in the arc combinatorics images. Second, the recovered images appear blurrier, which makes intuitive sense given the more rounded features of a limb-darkened transit event compared to those of a non–limb-darkened transit.

4.3.1. Arc Averaging

4.3.2. SART

The next algorithm we test is called SART (Andersen & Kak 1984). SART was originally developed for medical computed tomographic imaging. Specifically, SART reconstructs a 2D image from projections of X-rays through the body—this is directly analogous to our shadow imaging problem, where "projections" of the pixelated image against the stellar disk are the individual data points in the light curve.

SART operates iteratively upon an initial guess for the opacity vector ${{\boldsymbol{\tau }}}_{\mathrm{half}}$, which encodes the opacities of the pixels of the top half of the image grid. Beginning from this initial guess, it computes subsequent corrective updates to the individual entries of ${{\boldsymbol{\tau }}}_{\mathrm{half}}$.

The $\left(q+1\right)\mathrm{th}$ iteration of τ_l, the opacity of pixel l, is given by

$$\begin{eqnarray}&&{\tau }_{l}^{q+1}={\tau }_{l}^{q}+\displaystyle \frac{{\displaystyle \sum }_{k=1}^{{L}_{\mathrm{half}}}\left({B}_{{kl}}\tfrac{{C}_{k}-{\displaystyle \sum }_{\lambda =1}^{{L}_{\mathrm{half}}}({B}_{k\lambda }\,{\tau }_{\lambda }^{q})}{{\displaystyle \sum }_{\lambda =1}^{{L}_{\mathrm{half}}}{B}_{k\lambda }}\right)}{{\displaystyle \sum }_{k=1}^{{L}_{\mathrm{half}}}{B}_{{kl}}}.\end{eqnarray} \tag{ 42 }$$

The scalar ${\tau }_{l}^{q}$ is the lth entry of ${{\boldsymbol{\tau }}}_{\mathrm{half}}$ at iteration q, representing the opacity of pixel l; the scalar B_kl is the kth-row, lth-column entry of ${\mathsf{B}};$ and the scalar C_k is the $k\mathrm{th}$ entry of ${\boldsymbol{C}}$. ${\mathsf{B}}$ and ${\boldsymbol{C}}$ are defined in Equation (35).

For intuition, the update term in Equation (42) is equal to the average correction to pixel opacity τ_l over all rows and all columns of ${\mathsf{B}}$. (Hence, the sum in the denominator is over all rows of column vector ${{\boldsymbol{B}}}_{l}$, and the sum in the numerator term's denominator is over all columns of row vector ${{\boldsymbol{B}}}_{k}$.) The numerator, specifically, is the average value over all pixels in the grid of a sort of "residual" between the observed light curve and the model. This residual is equal to C_k minus the scalar projection of ${{\boldsymbol{\tau }}}_{\mathrm{half}}^{q}$ along ${{\boldsymbol{B}}}_{k}$. In effect, these two averages allow for a correction to the opacities, which is averaged over all time steps of the light curve and all pixels in the grid.

By running the SART algorithm for a large number of iterations (usually ∼10⁴ for a 16 by 16 pixel grid), we achieve good convergence to the observed light curve for a number of test cases. The rms error between the light curve of the input image and the light curve of the SART solution declines monotonically over the SART iterations, indicating that SART achieves a progressively better fit to the light curve as it proceeds.

We find that starting from an initial guess of all τ_l = 0.5 works well, because the step-by-step updates to ${\boldsymbol{\tau }}$ are generally small, so the algorithm does not wander far into unphysical parameter space (i.e., τ_l < 0 or τ_l > 1). In the event that the resulting SART solution does have slightly unphysical opacities, we redistribute the excess positive or negative opacity uniformly among the pixels whose centers fall within a distance of w/2 of the arc pair that intersects at the unphysical pixel. This redistribution renders the SART solution fully physical without drastically changing its light curve. Because SART exploits information in the light curve, not just its derivative, it is not necessary to renormalize the SART solution's pixel opacities.

SART solutions exhibit horizontal symmetry as a result of the flip degeneracy and semicircular arc–like features as a result of the arc degeneracy. Like the arc averaging algorithm, SART tends to smear out sharp features in the true input image along arcs, resulting in pixel grid solutions which are smoother, with more dispersed opacity than the true image. (SART solutions are even smoother than the corresponding arc-averaged solutions.) As as a result, SART fails, for example, to match the sharply flat-bottomed transits of the 16 by 16 pixel circle and square test grid light curves (Figure 9), producing slightly rounded light-curve shapes instead. On the other hand, because SART allows the pixel opacities to take any continuous value between 0 and 1, it can accurately reproduce the light curves of nonbinary test grids, like the planet–moon, planet-ring, and comet test grids.

4.3.3. Brute-force Search

4.3.4. Parsimonious Opacity Assignment

4.3.5. Arc Combinatorics

We begin with the TRAPPIST-1c/e/f triple transit, for which the expected shadow image is known. We hope to recover an image of three transiting planets, analogous to the diagram presented in Gillon et al. (2017, Extended Data Figure 1).

In attempting to invert this light curve, we have useful prior information beyond the expected image. First, because of the repeated transit observations and N-body dynamical simulations presented in Gillon et al. (2017), the periods and eccentricities, respectively, of planets c, e, and f are well constrained. This enables us to calculate the Keplerian orbital velocities of c, e, and f, which we can use as the v of our transiting grid. (We note that because these three orbital velocities are different, the pixel image we are attempting to recover changes during the transit, so we will only be able to recover an approximate image for any single choice of v.)

Second, because the physical behavior of this system is so well understood and the other properties of these planets (${R}_{{\rm{p}}}/{R}_{\ast }$, b, $a/{R}_{* }$) are so well constrained by transit modeling, we can generate an extremely finely time-sampled model light curve, based on a BATMAN model (Kreidberg 2015), of this triple transit, which matches the observed light curve. We can use this high-resolution light curve to test the effects of grid resolution on the success of shadow imaging: when the light curve is finely sampled, we can recover a much higher-resolution grid than when the light curve is sparsely sampled. Finally, we can adopt the same approach to determining the quadratic limb darkening coefficients for TRAPPIST-1 as Gillon et al. (2017) did in their analysis, interpolating TRAPPIST-1's stellar properties from the tables of Claret et al. (2012).

In Figure 10, we present three inversions of the BATMAN-modeled high-resolution TRAPPIST-1c/e/f triple transit light curve, conducted with the grid v equal to the Keplerian velocity of planets c, e, and f, respectively. We choose N = 16 because it is a high enough resolution that pixel width $w\lessapprox {R}_{{\rm{p}}}/{R}_{\ast }$ for planet e, the smallest planet (${\left({R}_{{\rm{p}}}/{R}_{\ast }\right)}^{2}\,=0.52$, according to the transit modeling of Gillon et al. 2017). We show the results of the arc averaging algorithm here, because it produces the cleanest and most interpretable shadow images, although results from the other three algorithms are qualitatively similar.

Zoom In Zoom Out Reset image size

In the shadow images, which transit the star moving left to right (i.e., the pixels at the right-hand edge of the image transit first), clear ingress and egress arcs for each planet are visible, in the expected order: first, planets c and f ingress together; then, e ingresses; c egresses; f egresses; and e egresses.

The three planets move at three different velocities to produce the light curve, but the grid moves as a unit, so none of the three inversions perfectly matches the light-curve model. When the velocity is correct for a particular planet, that planet's image is a pair of arcs whose points of intersection fall at the planet's impact parameter as measured by Gillon et al. (2017), demonstrating that shadow imaging of that planet is successful within the constraints of the arc and flip degeneracies.

When the grid v is slower than the planet's velocity, the planet's ingress and egress arcs are spaced too closely together; this effect is most visible for planet c in the top panel, where the grid moved at planet f's velocity. In the light curve, the overlapping arcs manifest themselves in a too early dip, caused by c's egress arc entering too quickly, and in a too deep transit depth between the egresses of c and f, caused by c's ingress arc remaining in front of the star for too long.

Conversely, when the grid velocity is faster than the planet's, the planet's arcs are too widely separated; this effect is most visible for planet f in the bottom panel. This time, the light curve is too shallow between the egresses of planets c and f, because f's ingress arc egresses too soon.

We next investigate what happens if we attempt to invert the observed Gillon et al. (2017) light curve of this triple transit, which is noisy and much more coarsely time-sampled, rather than a high-resolution BATMAN-modeled light curve. Additionally, we ask what happens if we attempt to recover a shadow image without knowing the true velocity of the transiting object: what happens if we use the guidelines presented in Section 4.2 instead?

We invert the observed TRAPPIST-1 triple transit light curve at a range of velocities: the slowest is 31.9 day⁻¹, corresponding to 2 divided by the entire triple-transit event duration (in accordance with the guidelines presented in Section 4.2), and the fastest is 135.9 day⁻¹, corresponding to 4 divided by the duration of planet c's transit by itself. At each velocity, we choose the maximum grid resolution N that, when combined with v to solve for M, allows the transiting grid to partially overlap the star at all time steps of the light curve while still maintaining NM at less than the number of observed data points. Accordingly, the resolution N decreases as v increases, because M increases with v to maintain full light-curve coverage.

In Figure 11, we present the results of these inversions. There are a number of interesting features about these results. We note, first of all, that SART is consistently the most successful inversion algorithm—this is true across the range of tested grid velocities. Furthermore, the SART shadow image consistently resembles the expected shadow image illustrated in Figure 10, even at low image resolutions. Arc combinatorics is somewhat successful at matching the observed light curve at the slowest tested velocity (corresponding to the highest grid resolution) but fails otherwise.

Zoom In Zoom Out Reset image size

The other algorithms fail consistently across the range of tested velocities. For arc parsimony and arc combinatorics, this results because these algorithms assign binary opacities (0 or 1) to individual pixels, rather than assigning continuous opacities. (We note that the shadow images presented in Figure 11 do not have binary opacities because the final step of both the arc parsimony and arc combinatorics algorithms is to average the binary shadow images produced at each time step dt and renormalize the average to match the observed transit depth.)

When the pixel resolution of the grid is too low, a single pixel's transit depth can exceed the transit depth of a shallow event like the TRAPPIST-1c/e/f triple transit (maximum transit depth ∼2%). As a result, the smallest unit of opacity that the arc parsimony or arc combinatorics algorithm can assign is too deep, and these algorithms cannot reproduce the observed light curve. Instead, they tend to assign opacity to pixels along the top and bottom of the image grid, which have the smallest impact on the light curve. This is especially visible in the high-v arc combinatorics panels in Figure 11.

Meanwhile, the arc averaging algorithm also fails to match the observed light curve, regardless of velocity. This is because the arc averaging algorithm, unlike SART, is not robust to noise in the light curve; noise is tantamount to light-curve fluctuations at much higher frequency than can be accommodated by the grid velocity. While SART is able to average out high-frequency noise over many corrective iterations, arc averaging calculates only one arc arrangement per time step dt; if these arrangements are wildly different for neighboring time steps, as they will be for noisy light curves, arc averaging fails.

From these investigations, we conclude that SART is the most robust light-curve imaging algorithm. In particular, light curves with large measurement uncertainties and/or shallow transit depths should only be inverted with SART.

Next, we proceed to a light curve with an unknown generative shadow image: that of KIC 8462852, Boyajian's Star (Boyajian et al. 2016). This star exhibits aperiodic, deep transit events of unknown origin; hypotheses to explain these events include a family of transiting comets (Bodman & Quillen 2016; Boyajian et al. 2016); circumstellar rings (Katz 2017); an intervening occulter not orbiting Boyajian's Star directly, such as a structure in the interstellar medium or an object with a dusty disk (Wright & Sigurdsson 2016); circumstellar debris following the star's earlier engulfment of a planet (Metzger et al. 2017); and alien megastructures orbiting the star (Wright et al. 2016).

We focus on the two deepest dimming events observed in the Boyajian's Star light curve during the Kepler mission, which Boyajian et al. (2016) named Dip 5 and Dip 8. Because these events are aperiodic, the appropriate grid velocity v is not obvious; furthermore, in both dips, the light curve smoothly tapers to a sharp point, so the "beginning" and "ending" points of the event are not obvious. Correspondingly, we start from a velocity $v=2/{t}_{\mathrm{event},\max }$, where ${t}_{\mathrm{event},\max }$ is a wide window around the deepest part of the transit, outside of which the flux of the star has essentially returned to 1. (These time ranges are plotted in Figures 12 and 13.) We then test several other velocities doubled from this starting point. We interpolate quadratic limb darkening coefficients for Boyajian's Star from Sing (2010).

Zoom In Zoom Out Reset image size

The inverted images for Dips 5 and 8 are presented in Figures 12 and 13, respectively. There are a number of interesting features in these images.

Zoom In Zoom Out Reset image size

First, there is a circular ring–like feature, of the same radius as the star, that appears generally in images from all four algorithms when the grid velocity v is too slow. (See Dip 5, $v\leqslant 1.6$ days, and Dip 8, $v\leqslant 1.2$ days, for examples.) This happens because when the grid velocity is too slow, the grid struggles to produce narrow features in the light curve; the rate of change of the state of the grid overlap is simply too slow. Under this constraint, the circular ring is the grid pattern that matches a narrow light-curve feature best, in the sense that it generates the narrowest possible V-shaped light-curve feature.

For intuition, consider a copy of the circular ring with the addition of some opaque interior pixels: at some time steps, these interior pixels will be entirely contained within the stellar disk, and their effect on the light curve during these time steps will be constant. In other words, their transit will be flat-bottomed. This is not the case for the ring, whose overlap state changes at every time step of the transit; the ring is the "opposite" of a flat-bottomed semicircular arc pair in this sense.

Velocities which produce a ringed shadow image are therefore too slow. This is also obvious from the light curve of the shadow image, which is wider than the observed transit event.

When v is fast enough, we find that all four algorithms produce qualitatively similar shadow images for both Dip 5 and Dip 8 and furthermore that the shadow images of the two dips are similar to each other. As in the case of the TRAPPIST-1 triple transit, the arc parsimony and arc combinatorics algorithms generate "noisy" shadow images where the light curve is shallow, because they cannot assign opacity in units smaller than one fully opaque pixel. For Dip 5, SART is clearly the best match to the shadow image; arc averaging produces a light curve which is too narrow, likely because the near-transparent pixels farther from the center which would have produced the "wings" of the transit event have been averaged away in the last combination-and-normalization step of the algorithm. Meanwhile, for Dip 8, the shadow images from both arc averaging and SART match the light curve well.

We strongly caution that there is no straightforward way to interpret these images, for two reasons. First, these images are subject to both the flip and arc degeneracies; second, the grid resolution is low (N = 5 for Dip 5; N = 6 for Dip 8) because we are limited by the 30 minute cadence of Kepler observations, and technically this resolution is so low that the small-planet limb darkening approximation used to calculate the light curves of these shadow images is inappropriate. Nevertheless, these limitations should only affect the distribution of opacity among the semiopaque pixels in the shadow images of Figures 12 and 13; pixels which are fully transparent in the shadow images should remain so, even if we were to obtain a much higher-resolution time series of these events.

We therefore note that the "gaps" of near-zero opacity (i.e., nearly transparent regions) which symmetrically frame the opaque transiting blob at the center of the shadow images in Figures 12 and 13 suggest that structured occulters are responsible for Dips 5 and 8 of the light curve of Boyajian's Star. The gap structure appears to be necessary to produce a shadow image which matches the ingress and egress shape of both Dip 5 and Dip 8; we note, for example, that this gap structure is missing in the arc-averaged shadow image of Dip 5 at $v=3.2{{\rm{d}}{\rm{a}}{\rm{y}}}^{-1}$, and the ingress and egress features of the corresponding light curve are too sharp to match the observed light curve.