Best friends from now on: shaped pupils and Lyot coronagraphs

How a shaped pupil solves the "classical" Lyot coronagraph problem of monochromatic field cancellation.
Figure from Zimmerman, et al., 2016: Illustration of how a shaped pupil solves the “classical” Lyot coronagraph problem of monochromatic field cancellation.

An article I authored with my colleagues from the Princeton high-contrast lab has been published by the Journal of Astronomical Telescopes, Instruments, and Systems: Shaped pupil Lyot coronagraphs: High-contrast solutions for restricted focal planes.

For 15+ years, Lyot coronagraphs and shaped pupils have evolved in parallel as competing tools for producing a zone of destructive interference next to a star’s image, with the goal of revealing orbiting exoplanets. After all this time, it turns out they are much better together than apart.

We learned this fact not a moment too soon. The instrument designs for WFIRST are getting locked in place rapidly, and that mission may be our last chance for many years to launch an ambitious space-based exoplanet imager. The new hybrid shaped pupil Lyot coronagraph offers a winning balance of performance, fabrication maturity, and error tolerance.

Generate an interactive gallery of sky survey cutouts from a list of targets

After filtering down some source catalog to a short list of interesting targets, it’s nice to have image cutouts handy, as either a sanity check, or as finder charts. Aladin Lite, hosted by CDS, has a javascript interface that allows you to display image cutouts from some major sky surveys. Here’s a short python script I wrote to generate a gallery of DSS and 2MASS cutouts, with hyperlinks to the full SIMBAD entry. As an input, it takes an ASCII file containing a list of 2MASS catalog IDs.

Here is an example gallery; the python source is below (hosted with gist).

Poster at 227th AAS meeting

AAS meeting, Jan 2016, Kissimmee, Florida
At 227th AAS meeting, Jan 5, 2016, Kissimmee, Florida. Photo by Heather D’Angelo.

Here I am at the American Astronomical Society meeting this week. I presented a poster about the first project I began after arriving at STScI, in collaboration with Marie Ygouf, Laurent Pueyo, Rémi Soummer, and others. I analyzed laboratory data from JPL’s High Contrast Imaging Testbed, where protoypes for the WFIRST coronagraph are being tested. The goal is to learn how well we can subtract residual speckle noise from the coronagraph images, to predict what kinds of planets we can detect.

Imaging the planet around GJ 504

gj504_coadd_medres
Figure 2 from Skemer, Morley, Zimmerman, et al.: De-rotated co-add of all photometric quality GJ 504 images taken in the LNB7 (3.95 μm) filter. The planet, GJ 504 b, is the circled point source to the northwest of the star, at separation ∼2.′′5. Due to the large dynamic range of structure in the image, two regions are displayed with different flux scales: (1) within 2′′ of the star, the Airy rings are displayed with a logarithmic stretch; the outer region (2) that includes the planet is displayed with a linear stretch centered at zero. The corresponding grayscale bars on the right-hand side indicate the flux in units normalized to the peak of the unsaturated star’s PSF.

Andy Skemer, Caroline Morley, myself, and the rest of the LEECH Science Team posted a paper on the arXiv today: Characterization of the Coldest Directly Imaged Exoplanet, GJ 504 b, and Evidence for Super-Stellar Metallicity.

The article, soon to be published in ApJ, describes new observations of the planet candidate discovered to orbit GJ 504 a few years ago (Kuzuhara et al., 2013). By measuring the thermal emission from the companion in three custom filters in the range 3.6–4.1 microns, we’ve made new inferences about the temperature, density, and composition of the atmosphere. The observations were made possible with LMIRCam, an AO-fed infrared camera on the Large Binocular Telescope.

I was responsible for the entirety of the data reduction, taking steps familiar to other high-contrast observers in this wavelength regime:

  • Aligning and combing thousands of raw, short exposures, correcting for background sky noise and detector pattern noise
  • Angular differential imaging to subtract the star (here implemented with the KLIP algorithm)
  • Contrast curve calculation, PSF modeling, photometry, and error analysis

GJ 504 b is on the extreme faint end of all the planetary-mass companions imaged to date (as Figure 4 in the journal article illustrates). Fortunately, the combination of relatively high angular separation (2.5 arcsec from the host star) with the deep sensitivity of LMIRCam–as shown above in the co-add image–was high enough to secure a detection at each filter. The wavefront correction provided by the LBTI AO system was such that Airy rings spanning three orders of magnitude in intensity are evident, out beyond a 1-arcsec radius. Even without PSF subtraction, the planet is evident at a contrast of 1E-5, in the northwest corner of the co-add image.

There is a lot more we could learn from the infrared SED at wavelengths redward of 4 microns. For such a faint object, however, that part of the spectrum is not likely to be accessed from the ground. Therefore, in the future GJ 504 will make an excellent target for JWST NIRCam.

Barnard’s Magic Lantern

Among his many contributions to astronomy, Edward E. Barnard’s photographic studies of the galactic plane were the most transformative and enduring. The images themselves, posthumously compiled in A Photographic Atlas of Selected Regions of the Milky Way, remain stunning today. Aesthetically, at least, they hold their own against contemporary multi-color CCD mosaics of star forming regions. They are all the more impressive if we imagine the impact they had on the scientific community in the early 1900s. At that time, our awareness of the luminous structure of the galactic plane was severely limited by the weak sensitivity of the human eye to these diffuse patterns.
rho oph barnard

Above are a few snapshots I took of the copy of the atlas on the shelves of Princeton’s science library (not an original edition but a reproduction by Cambridge press). I was particularly drawn to the plate centered on rho Oph. Of this field, Barnard wrote,

The region of Rho Ophiuchi is one of the most extraordinary in the sky. The nebula itself is a beautiful object. With its outlying connections and the dark spot in which it is placed and the vacant lanes running to the east from it, it makes a picture almost unequaled in interest in the entire heavens.

Ironically, despite having access to the largest and most advanced observatories of his era, it was a 1.5-inch lens scavenged from a “magic lantern” that led to a major breakthrough in capturing the vast structures strewn along the Milky Way. A magic lantern was an oil-fueled slide projector, a common tool for popular amusement in the decades preceding the cinema. Vistas of emission nebulae and dark lanes, invisible through the 36-inch Lick refractor, illuminated the dry plates at the short focus of the lantern lens. In 1895 he published a description of his apparatus in the second issue of the nascent Astrophysical Journal:

I have elsewhere given some account of the performance of a small lens in photographing large areas of the sky, and for the delineation of very large diffused nebulosities, etc. This lens is a very small one belonging to an ordinary “magic lantern;” it is one and one-half inches in diameter and some four or five inches equivalent focus.

these six photographs with the small lens show us in a most striking manner how the most valuable and important information may be obtained with the simplest means.

In an era when most serious observational work was carried out on long refractors with narrow fields of view, no other astronomer had realized the scientific potential of a low focal ratio objective. Barnard’s engineering efforts in this direction culminated in the construction of the Bruce Telescope, an observatory dedicated to wide-field astrophotography. His wide-field telescopes were true predecessors of many survey-oriented optical observatories in use today, one century on: the Sloan 2.5-meter at Apache Point, the Palomar Transient Factory, Pan-STARRS on Haleakala, VISTA at Paranal, Wide-field Infrared Survey Explorer, etc. In the next decade, LSST and WFIRST will carry that torch.

In 1905, Barnard temporarily relocated the Bruce Telescope from its home at Yerkes Observatory to Mount Wilson in Pasadena. Under the dry Southern California skies he acquired the bulk of the photographs comprising the revered Milky Way atlas. Thanks in large part to these images, the old Herschelian view of the dark nebular features as “vacancies” in the heavens began to be dispelled. Although he was too cautious to entirely dismiss the prevailing theory, he published a classic article in 1919 describing the mounting evidence for the obscuring nature of many of the dark features: http://adsabs.harvard.edu/abs/1919ApJ….49….1B.

The way in which Barnard brought the Milky Way into focus is only one aspect of the man’s life and career I found deeply inspiring to read about. To learn more I recommend William Sheehan’s biography, The Immortal Fire Within.

Sparse aperture mask wavefront sensor

An article I co-authored with Princeton grad student Hari Subedi has been published by the Journal of Astronomical Telescopes, Instruments, and Systems: Coronagraph-Integrated Wavefront Sensing with a Sparse Aperture Mask . A sparse aperture mask, similar to the kind now in common use for high-angular resolution astronomical imaging (a recent example being Hinkley et al. 2015), is also effective at converting phase aberrations into diverse intensity patterns. By intercepting, collimating, masking, and finally refocusing the starlight rejected from a coronagraph focal plane, we therefore have a simple approach to sensing dynamic low-order aberrations of the kind that are problematic for space telescope coronagraphs such as the one envisioned for WFIRST-AFTA.

Fig2_SAM_WFS_diagram_v2_lores

A shaped pupil for a backyard telescope

Dr. Russ Genet at Cal Poly does speckle interferometry observations of binary stars. He came up with the idea to use a shaped pupil to improve the search zone around his primaries, and increase the kinds of targets accessible to his small telescope observing program. Our team agreed it would be an interesting experiment to see what shaped pupils can do for seeing-limited observations of binary stars.
I wrote some 2-D optimization programs (making yet more use out of the approach of Carlotti et al 2013) to come up with a design for Russ’ 11-inch Schmidt-Cassegrain pupil. There was one solution with a fairly simple layout. It improves the open aperture telescope PSF by suppressing the first sidelobe along one axis, as shown in the horizontal contrast cut plot below.C11_SP_12wa60_logc27_bowtie90degC11_contrast_cut_12wa60_logc27_bowtie90degIf Russ and his observing team, led by Cal Poly Master’s student Ed Foley, can freeze out the atmosphere and limit static aberrations, then they have a shot of effectively (almost) doubling the aperture of their Celestron setup.

naive simulated TPF coronagraph image

A student at Princeton is doing his senior thesis on simulated coronagraph data sets, to use as a testing ground for hypothesis testing of somewhat realistic (low SNR) planet signals.

Fomalhaut is the target star for the case study. For its distance of 7.61 pc, 1 AU and 5 AU projected separations correspond to 131.4 mas and 657.0 mas, respectively. For 8-meter telescope observing at 550 nm, \lambda/D = 14.2 mas, so maximum angular separations are 9.27 \lambda/D and 46.33 \lambda/D. Periastron for the inclined, circular 5 AU orbit is \cos(66^{\circ})\times 46.3 \lambda/D = 18.85 \lambda/D.
I handed off a Matlab Fourier propagation model of a concentric ring (spiderweb) coronagraph to my student. This design has a working angle range 6 - 60 \lambda/D to frame the system nicely. Also note that 1 AU projected separation at 1.34 pc distance of alpha Cen corresponds to 1/1.339 / 0.0142 = 52.6 \lambda/D.

concentric_ring_SP_6wa60_c10
This design surpasses 10^{-10} contrast for the ideal wavefront, so of course it’s fun to throw in the canonical Earth-like dot, shown below to the right of the star (whose PSF “saturates” the flux scale).
sun_plus_earth_PSF

Floating point roundoff errors in Fourier optics calculations

Discrete Fourier transforms are essential for simulating light propagation in telescopes and coronagraphs. To cut down the computational burden, and maintain high spatial resolution, we can make use of well-established symmetry properties of Fourier transforms. For example, if we start with a real-valued electric field in the telescope pupil, with even symmetry about the y-axis, we should get an image plane electric field with even symmetry about the y-axis. Or so I thought, until a sanity check failed…
fft2D_real_result

Here’s where my assumptions broke down. The above plot is the real part of the image formed by an ideal shaped telescope pupil. I computed this with a 2-D FFT in Matlab and shifted the result so that the zero-phase pivot is at the array center. Visually, it appears symmetric about the vertical meridian, so that \mathfrak{Re}\left\{ E(x, y) \right\} = \mathfrak{Re}\left\{ E(-x,y) \right\}. And that is indeed what we’d expect, transforming from a flat (zero-phase) pupil plane wavefront with the same even symmetry.

So when I take the right half of the image, and then flip and subtract it from the left half, I should get a wall of zeroes. Instead, disparities appear:

Efield_even_symmetry_check_2Dresid

We have disagreements as high as 2000, corresponding to almost 1% of the peak field amplitude in the image plane. Initially, I was very puzzled by what might cause such large errors. Was it an error in my pupil definition, or some array offset creeping in? Do I not understand the rules of Fourier transforms, after all? After some tests, I figured out that in each image plane pixel I was accumulating a large number of tiny floating-point roundoff errors.

\displaystyle X [ k, l ] = \frac{1}{\sqrt{MN}}\sum_{m=0}^{M-1} \sum_{n=0}^{N-1} x[m, n] e^{-i2\pi(mk/M + nl/N)}

Consider the above textbook formula for the 2-dimensional discrete Fourier transform. There are countless applications, but to an astronomer it describes the electric field diffraction pattern in the image plane of a telescope when the incoming starlight is monochromatic, unpolarized, and the optics are perfect. We can break this into inner and outer sums:

\displaystyle X [ k, l ] = \frac{1}{\sqrt{MN}} \sum_{m=0}^{M-1} \left( \sum_{n=0}^{N-1} x[m, n] e^{-i2\pi nl/N} \right) e^{-i2\pi mk/M} \\ = \frac{1}{\sqrt{MN}} \sum_{m=0}^{M-1} P[m, l] e^{-i2\pi mk/M}

In our case, the telescope pupil x[m,n] has even circular symmetry along the x-axis, so that x[m,n] = x[m,N-n]. When that is the case, and x[m,n] is real, it can be shown that the inner sum, P[m,l], has circular even symmetry: P[m,l] = P[m,N-l]. Furthermore, so does the image: X[k, l] = X[k, N-l].

The problem is, in practice there are several orders of magnitude of separation between the mirrored spatial frequency pair l and N-l. Before the complex exponential can be evaluated, its phase argument must be stored as a floating point value, which entails some roundoff error. To investigate my case, I chose some sample points to test how the inner complex exponential term behaves for my array (with size = 3000 x 3000).

With N = 3000, l = 10, I noticed a particularly high roundoff error at n = 1844. Since I only consider the real part for now, I simply compare cos(2 \pi nl/N) for this mirrored pair of spatial frequencies:

>> wrapToPi(2*pi*n*l/N)
ans =
    0.9215
>> wrapToPi(2*pi*n*(N-l)/N)
ans =
   -0.9215
>> wrapToPi(2*pi*n*l/N) + wrapToPi(2*pi*n*(N-l)/N)
ans =
  -1.1866e-12
>> cos(2*pi*l*n/N) - cos(2*pi*(N-l)*n/N)
ans =
   1.5709e-12

Because the two wrapped angles are symmetric about zero radians, or \left( nl/N~\%~1 \right) + \left( n(N-l)/N~\%~1 \right) = 1, the difference of cosines is mathematically zero. What we actually compute is the best approximation possible with 64-bit double precision floating point. The residual comes out exactly the same in both MATLAB and Python/Numpy. By contrast, if I take the difference of another pair of cosine evaluations, but now with the higher angle wrapping around only one period beyond nl/N, the error is several orders of magnitude smaller:

>> wrapToPi(2*pi*n*l/N) - wrapToPi(2*pi*(n*l/N + 1))
ans =
   7.1054e-15
>> cos(2*pi*n*l/N) - cos(2*pi*(n*l/N + 1))
ans =
  -2.2204e-16

Although that 1e-12 above seems tiny, since we are adding up N phasors in the inner sum, and then adding up M of those, roundoff errors are accumulated 9 million times for each image pixel. Just left of the image center (the l=0 pivot axis of the DFT), N-l is at its highest. Conversely, just right of the image center, the roundoff error fed into the trig function will always be very low.

Smaller array sizes would reduce this problem. But depending on what kind of measurement you’re making, that may not be an option. A better solution is to exploit the pupil’s symmetry and only keep the right hand side of the FFT result (spatial frequencies l = 0, 1, 2,\dotsc, N/2-1). That would maintain a self-consistent symmetry between the pupil and the image. In fact, because the full complex image field of a real pupil is Hermetian, in this particular case the full image image is determined by one quadrant; the imaginary part in the lower half reflection just needs to be negated.

That leads to the question: is their an FFT library for MATLAB or Python/Numpy that either automatically detects, or allows you to specify the symmetry properties of your input array? It seems like an obvious way to eliminate wasteful operations which are particularly prone to quantization noise.

Some useful references on DFTs around the web:
http://mathworld.wolfram.com/DiscreteFourierTransform.html
http://fourier.eng.hmc.edu/e161/lectures/fourier/node11.html
http://web.eecs.umich.edu/~fessler/course/451/l/pdf/c5.pdf