Translator Disclaimer
Open Access Paper
26 July 2016 Point spread function determination for Keck adaptive optics
Author Affiliations +
One of the primary scientific limitations of adaptive optics (AO) has been the incomplete knowledge of the point spread function (PSF), which has made it difficult to use AO for accurate photometry and astrometry in both crowded and sparse fields, for extracting intrinsic morphologies and spatially resolved kinematics, and for detecting faint sources in the presence of brighter sources. To address this limitation, we initiated a program to determine and demonstrate PSF reconstruction for science observations obtained with Keck AO. This paper aims to give a broad view of the progress achieved in implementing a PSF reconstruction capability for Keck AO science observations.

This paper describes the implementation of the algorithms, and the design and development of the prototype operational tools for automated PSF reconstruction. On-sky performance is discussed by comparing the reconstructed PSFs to the measured PSF’s on the NIRC2 science camera. The importance of knowing the control loop performance, accurate mapping of the telescope pupil to the deformable mirror and the science instrument pupil, and the telescope segment piston error are highlighted. We close by discussing lessons learned and near-term future plans.



W. M. Keck Observatory (WMKO) was the first to implement both natural guide star (NGS) and laser guide star (LGS) AO systems on a large telescope in order to achieve angular resolutions in the near-infrared that match the angular resolution of the Hubble Space Telescope at visible wavelengths. A total of 633 refereed science papers have been published through 2015 using the Keck AO systems. WMKO has endeavored to continually improve the capabilities of these systems.

One of the challenges for AO science instruments is the limitations due to incomplete knowledge of the PSF making them difficult to use AO for accurate photometry and astrometry in both crowded and sparse fields, for extracting intrinsic morphologies and spatially resolved kinematics, and for detecting faint sources in the presence of brighter sources. The characteristics of the PSF vary with observing condition, instrumental parameters, and guide star properties and hence the need for PSF estimation for each science exposures. Moreover, the PSF depends on the position in the science field due to field dependent instrument distortions and anisoplanatism due to the difference in turbulence between the guide star and the science target.

The goal of the PSF determination program at WMKO is to provide a PSF estimate for every point in the science field through post processing of the AO telemetry data taken in parallel with the science observations. The resulting capability of producing a high-fidelity PSF estimate will provide dramatic science gains and stands as one of the next major challenges for achieving further breakthroughs in high angular resolution science.

The PSF determination efforts at WMKO can be broadly classified into two categories: on-axis and off-axis cases. In the on-axis case, the guide star is co-aligned with the science object. In the off-axis case the natural guide star (NGS) is located within an ~30″ radius field around the science target for NGS AO and an ~60″ radius field for laser guide star (LGS) AO. Both the on-axis and the off-axis cases could use a NGS or a LGS. The off-axis LGS could be further divided into (1) tip/tilt NGS off-axis and LGS on-axis, and (2) tip-tilt NGS and LGS off-axis (the off-axis angle for the tip/tilt NGS and the LGS need not be the same). The on-axis PSF determination is led by WMKO with a sub-contract to the University of Applied Sciences Western Switzerland (Hes-so) and is funded by the National Science Foundation (NSF). The off-axis PSF determination effort is led by UCLA, including a sub-contract with the Optical Sciences Company, and is funded by the W. M. Keck Foundation. The NSF funding also supports designing an operational tool for PSF determination that would make use of both PSF determination efforts. More recently, we have initiated technical collaborations with Marcos van Dam at Flat Wavefronts, Luc Gilles and Lianqi Wang at the TMT, and Carlos Correia and Olivier Martin at LAM in validating the PSF reconstruction algorithm. This paper gives a broad view of the progress achieved in implementing this capability for Keck AO science observations, with a special emphasis to the on-axis NGS case and the work performed since our last update at the AO4ELT4 conference in October 2015 [1][2].

The outline of the paper is as follows: The concept of the PSF reconstruction is briefly described in Section 2, core on-axis algorithm development in Section 3, operational tool/facility development in Section 4, recent AO performance optimization/characterization relevant to PSF reconstruction in Section 5, highlights of the on-axis results in Section 6, off-axis algorithm development in Section 7, science verification plan in Section 8, and a brief summary and future plans in Section 9.



The PSF is the response of an imaging system to a point source. The image of a complex astronomical object can be seen as a convolution of the true object and the instrumental PSF. There are several terms involved in defining the complex astronomical optical system and the uncorrected atmospheric turbulence.

The approach taken at WMKO involves (1) computing the guide star PSF from wave front sensor (WFS) measurements using the technique introduced by Véran et al. [3] for a curvature WFS and modal correction AO system (PUEO on the CFHT), which is applicable to a Shack-Hartmann WFS and zonal correction system such as the Keck AO systems [1][2][4], and (2) computing the science PSF by applying anisoplanatic corrections to the guide star PSF using a modified version of the Arroyo software package written by Matthew Britton [5][6] and off-axis instrument characterization by the UCLA team [7].

Assuming that the residual phase is stationary over the telescope pupil, the long exposure Optical Transfer Function (OTF) can be written as the product of the instrumental OTF (telescope and instrument optics) and an OTF associated with the residual post-AO wavefront error, i.e.


The long exposure PSF is estimated from the reconstructed OTF through an inverse Fourier transform. The OTFAO is related to the residual phase structure function, DAO (νλ) as follows:


The wavefront phase errors are assumed to be composed of complementary and orthogonal components, namely, controlled and uncontrolled modes. The residual errors of the controlled modes (tip/tilt and higher order or deformable mirror, DM, errors) are estimated from the AO control loop data, and the errors of the uncontrolled modes (fitting error and aliasing) are estimated through modeling of astronomical seeing knowing the instrument characteristics. i.e.


Where λ is the wavelength of science observation, ν is the angular frequency vector on the sky, <εiεj> is the covariance of the WFS measurements, Ui,j is the spatial correlations of the DM modes (the influence functions), and DTT (νλ), DFE (νλ) and DAL (νλ) are the tip/tilt, fitting error and WFS aliasing structure functions, respectively.

Aberrations not seen by the wavefront sensor can be estimated through on-sky phase diversity measurements. In the case of LGS AO, focal anisoplanatism must be included and for off-axis cases at least tip/tilt anisplanatism is required. For the cases where the LGS is off-axis, angular anisoplanatism is also needed. The anisoplanatism terms are estimated through Cn2 profiles from an on-site MASS/DIMM atmospheric profiler [1][2][4].

A schematic diagram illustrating various components of the PSF reconstruction process is shown in Figure 1. The top, middle, and the bottom levels show various components developed at Hes-so (Haute école spécialisée de Suisse occidentale), WMKO, and UCLA, respectively.

Figure 1:

A schematic diagram of the PSF reconstruction process showing various components developed at different locations.




The core on-axis software consists of three major components, namely, (1) the atmospheric seeing estimation tool, (2) static low order aberration estimation tool, and (3) the PSF reconstruction tool. The first tool is essential for the second and the third, and the second one is essential for the third.


Atmospheric Seeing Estimation Tool

The seeing tool estimates the Fried parameter and the outer scale through statistical analysis of the commands to the DM of the AO system during closed loop operations. The principle of the method was proposed by Véran et al. [3], and demonstrated on a modal system (PUEO, a curvature sensor based system). We extend this method to retrieve the outer scale of turbulence as well. We find that measuring the Fried parameter and outer scale during science observations along the direction of the science target is crucial to achieve the necessary accuracy on the properties of the reconstructed PSFs as these parameters are non-stationary, both in time and space. Figure 2 shows that the Fried parameter estimated from the closed loop AO telemetry compares well with the values estimated from long exposure NIRC2 images. The implementation of the seeing tool for our PSF reconstruction program is presented in two journal papers [8], [9].

Figure 2:

Fried parameter estimated from closed loop AO telemetry is compared to the value estimated from long exposure NIRC2 images.



Static aberration estimation

The achieved Strehl ratio of the Keck AO systems under excellent seeing condition are limited by static phase aberrations introduced by the optical elements of the telescope and the instrument, calibration errors, and the inability of Shack Hartmann wavefront sensors to measure segment piston errors. The phase aberrations could be measured from the science images as was demonstrated for the Hubble space telescope (HST) [10, 11, 12] although measuring the phase errors in the presence of the residual atmospheric turbulence is much harder.

The phase-diversity technique [13] was chosen to measure the low order static aberrations in the system as this has the potential to be used with non-point sources, as this is the case for on-sky objects due to residual turbulence, and no additional hardware is involved in making defocused images. The defocused images are taken by moving the wavefront sensor focus stage typically by 5 mm (0.6 waves) for measurements taken at 1.6455 μm. This simplistic approach of not simultaneously taking in-focus and out-of-focus images poses some challenges in practice – especially under varying seeing conditions. We take about 20 long exposure images each for the in-focus and defocused observations and form in- focus and out-of-focus pairs taken under similar seeing condition. A journal article on the implementation of the phase diversity scheme at WMKO is being prepared.

On-sky in-focus and out-of-focus NIRC2 data, taken on 2013 Feb. 3 UT, is presented below to illustrate the performance of the algorithm. The measurements were taken using the Fe II narrow band filter (effective wavelength = 1.6455 μm) at three defocus values: -2 mm, -4 mm & -6 mm, corresponding to an rms wavefront error of 0.24, 0.48 & 0.72 waves, respectively. The retrieved phase maps are shown in Figure 3. The rms WF error of the three phase maps (left to right) are 97, 105 & 102 nm, respectively. The consistency of the phase maps at different defocus value is apparent. Also, the model and the zonal reconstructions give similar results. The reconstruction code produces large and unphysical excursion of the phase on the very edge of the segmented pupil, and we are working on a solution to mitigate this issue. Figure 4 compares the retrieved average phase map from 2013 Feb. 3 UT with that of 2013 Aug. 01 UT. The similarity between the two phase maps suggests that there is a component of the static aberration that is stable over at least six month period.

Figure 3:

The retrieved phase map for the three epochs using the phase diversity algorithm. Left to right: phase diversity measurements taken with defocus values -2 mm, -4 mm & -6 mm (0.24, 0.48 & 0.72 waves at 1.6455 μm), respectively.


Figure 4:

The retrieved phase map for the two engineering runs (2013 Feb. 03 & 2013 Aug. 01) using the phase diversity algorithm.



PSF reconstruction

The PSF reconstruction tool, developed in IDL for NIRC2 science observations, uses AO telemetry data saved as IDL save files in addition to instrument configuration information stored in the fits header of the science exposures. The PSF reconstruction program does not look at the science data as the objective is to predict the science PSF. The MASS profiler data is required to measure focal anisoplanatism for the LGS cases, and tilt and angular anisoplanatisms for the off-axis cases. The MASS data is currently manually downloaded from the Mauna Kea Weather Center website but efforts are underway to automatically acquire them. The PSF reconstruction tool gets the Fried parameter by calling the seeing tool (Section 3.1) and the instrument OTF from the phase-map generated using phase diversity technique (Section 3.2). The tool also requires the Uij matrix (the spatial correlations of the deformable mirror (DM) modes - the influence functions) already computed for the Keck WFS configuration and stored as an IDL save file.

The main task of the tool is to generate AO OTF for a given science exposure by computing various structure functions, namely the (1) tip/tilt (2) DM, (3) WFS aliasing, (4) fitting error, (5) focal anisoplanatism, (6) tilt anisoplanatism, and (7) angular anisoplanatism structure functions. The tip/tilt, focal, and angular anisoplanatisms have also been implemented at the WMKO for testing purposes. We will be comparing these terms with that developed at UCLA and have the option to use either of the implementation in the final operational tool.

The DM structure function is computed from the covariance of the WFS measurements using the Uij matrix. The fitting error is computed using a Monte-Carlo algorithm where instantaneous turbulent phase screens are generated, and projected on to the influence function basis made from an empiral model of the DM influence function. The projection is subtracted from the initial turbulent phase, and the instantaneous residual phase power spectral density (PSD) is computed. The procedure is run until a convergence is reached on the long exposure PSD. The structure function is then computed from the 2D correlation of the residual phase computed from the Fourier transform of the PSD. The final OTF is computed using Eq. 1 through 3.



In addition to the core algorithm software, several tools have been developed at WMKO with the goal of transitioning the PSF demonstration project into a facility class operational tool. An initial version of the prototype operational tool was developed in the IDL environment through multiple scripts. In order to provide more stability, a database and a task scheduler were implemented subsequently as a part of the prototype. The effort involves (1) modifications to the AO telemetry archival process, (2) data management and workflow development (including a task scheduler, mySQL database, and an interactive web-based graphical tool), and (3) computational facility development (including storage disks/tapes and servers) [2].


AO Telemetry Archival

The telemetry recording system (TRS) records telemetry data on AO nights. Data is sampled at 10 MHz and saved to a PostgreSQL database. The amount of data stored in the database is large enough such that the database disk has to be cleared on a rolling seven day cycle. To preserve the telemetry data for a longer period of time and to reduce the disk space requirements, the data is selectively retrieved and transferred to disk storage as described below.

A script retrieves and saves TRS data for the integration period of a given science FITS file. The script is started via a cronjob and runs daily. The script uses IDL routines to find the appropriate FITS files, verify that they are on-sky images and that the AO loops are closed, and then retrieve the TRS data. The retrieval of the TRS data involves querying the database through a C-interface to the PostgreSQL database. The processing tasks are currently running on a storage server.

The cronjob starts at 8 am daily and the data extraction typically takes a few hours depending on the amount of science data collected the previous night.


Data Management and workflow

To facilitate the processing of archived and new data, a prototype software infrastructure has been developed to collect, ingest, monitor and manage the required science and telemetry data. PSF reconstruction (PSF-R) tasks are scheduled for processing immediately after ingestion into the system. A workflow control task dispatches the scheduled tasks to multiple hosts; depending on available resources up to 48 tasks can be executed in parallel. The task’s progress and status can be monitored and displayed via a web-based graphical tool, which can also be used to cancel tasks or select tasks for re-processing. The tools use MySQL, PHP and JavaScript database/language. These tools enable easy monitoring of the status of telemetry extraction and PSF-R processes, and enable easy rerun of telemetry extraction and PSF-R processes, if necessary. The concept of scheduling multiple tasks (one task per each science exposure) on multiple computers across the Keck network guarantees the completion of PSF-R of the science data from the previous night before sunset.


Computational Facility and Data Storage

The computational facility includes (1) the servers used for computation, and (2) the storage system to accommodate at least ~3 years of data readily accessible for PSF-R.

The PSF reconstructions are performed on two computer systems in parallel across the internal network: (1) a virtual computer system in Linux environment, and (2) a dedicated storage server, a Linux host, that holds the AO telemetry data and some Keck Observatory Archive (KOA; [14]) data.

The AO telemetry data collected during science exposures are archived to the PSF-R data storage system. Re-processing of the existing/old data may be required if the PSF-R application is modified and hence we would like to have the data readily available online on the network. The storage requirement is ~ 10 TB/year/telescope to hold about 3 years of data, therefore a total of 60 TB total storage for Keck 1 and Keck 2. As the storage cost tends to go down with time, we started off with half this capacity.

The nightly AO telemetry data are transferred to the storage disk before noon (HST) the following day. This provides the necessary time in the afternoon for lev0 data processing. We opted for an ASL Sovereign 4898SRT storage server with ~ 33 TB of useable space. The 33 TB storage space is carved up into smaller volumes so that we can freeze the data and back it up to LTO-6 tapes for the offline storage requirement. A new LTO-6 tape drive has been procured and physically connected to the Sovereign 4898SRT server for data backup for disaster recovery purposes. The LTO-6 tape drive is a 2.4TB non-compressed data storage device; therefore we have configured Network File System (NFS) storage volumes to be approximately 2 TB each. The goal is to be able to store a NFS volume on a single tape.



We discuss some of the AO performance improvements being made to improve PSF reconstructions in addition to improving the image quality. A toolkit has been developed at WMKO to interface the core algorithm to the Keck software platform and to improve the accuracy of PSF-R. The software includes routines to (1) account for control loop delays and tip/tilt mirror dynamics (Section 5.2), (2) apply centroid gain corrections (Section 5.3), and (3) potentially measure high order static phase aberrations (Section 5.4).


Telescope pupil registration to the DM and NIRC2 pupil masks

The Keck telescope pupil was decentered with respect to the NIRC2 pupil masks by ~ 34 mm in x and ~ 478 mm in y on NIRC2 in telescope primary mirror space during the period most of the engineering data was taken. Additionally the telescope pupil nutated by ~ 159 mm peak-to-peak in the telescope primary space as the AO rotator rotates. A NIRC2 pupil mask was therefore not used for PSF-R observations. Subsequently the K2 AO bench was realigned to accurately map the telescope pupil to the DM and the NIRC2 pupil masks. This is crucial to apply the algorithm to science observations typically taken using one of the NIRC2 pupil masks.


Control loop delay verification and tip/tilt mirror dynamics update

One of the uncertainties in the PSF-R has been incomplete knowledge of the mirror dynamics and loop delays. Our earlier attempts to verify and/or improve the control loop model parameters both from on-sky data and internal data have been inconclusive primarily due to daytime beam-train vibrations and residual turbulence in the case of on-sky data.

More recently, the tip/tilt loop delays for different wave front sensor camera clocks have been validated through an internal test with carefully chosen intensity levels to have adequate measurement noise in the data; the data was also taken over the weekend to minimize vibrations from activities at the summit. We took telemetry data using the calibration source (SFP) and running the wavefront controller with different frame rates, programs and loop gains. These measurements were then compared with the modeled control law [15]. An excellent agreement is attained between the measured and modeled rejection transfer functions by modifying the model for the tip-tilt mirror dynamics. The updated model for the tip/tilt mirror dynamics on Keck II, obtained by trial and error, is as follows:


where s is the complex number frequency. A power spectrum of the tip/tilt residuals taken at 438 Hz with a very high loop gain of 0.7 using the artificial light source is shown in Figure 5. The measurements were taken at a low light level in order to be dominated by the measurement noise. The blue and red curves represent the rejection transfer function using the old and new tip-tilt mirror dynamics respectively. It can be seen that the new curve does a much better job of matching the measured power spectrum.

Figure 5:

Measured power spectrum when running at 438 Hz with a loop gain of 0.7. The modeled rejection transfer function using the old and new tip-tilt mirror dynamics are over-plotted in blue and red respectively.


A similar test was performed for the Keck I AO system and the updated model for the mirror dynamics on Keck I is as follows:


The loop delays for the STRAP tip/tilt sensor (used for the LGS case) are expected to be of the order of a few microseconds and as the sensor integration time is typically a few milliseconds, the latency for this case is negligible.

The DM loop delays have also been verified through DM loop noise modelling. The measured power spectra are compared with models obtained by filtering noise through the theoretical rejection transfer function. The results for a sample case taken at 1054 Hz is shown in Figure 6. The match of the control loop model to the data is not entirely satisfactory. The measured power spectrum shows sharp peaks and a broad hump centering around 375 Hz. We will be investigating the source of the strong peaks measured on Keck II and in understanding why the power spectra do not look as expected.

Figure 6:

Measured power spectrum of DM loop at 1054 Hz with a loop gain of 0.4. The modeled rejection transfer function is over- plotted in red.



Automated centroid gain update

As the tip/tilt loop delays for the wavefront sensor is now verified through laboratory experiments and the tip/tilt mirror dynamics has been updated, the control law modeling of on-sky tip/tilt residuals provides more reasonable estimation of centroid gain from on-sky data for the NGS case – especially for the cases where measurement error dominates over the bandwidth error. The control loop modeling serves the purpose of validating the automated centroid gain update tool based on the seeing estimation and daytime spot size calibration.

We have implemented an operational tool that addresses the source of the problem (i.e. automatically updates the centroid gain in real time based on the seeing). The updated centroid gain (CG) is estimated by quadratically adding contributions from the atmospheric seeing (ω0) at 640 nm and the static (daytime) spot size (0.646" – primarily due to charge diffusion) as follows:


Where FF is the fudge factor to be set through simulation and is currently set to unity. Figure 7 (left) shows the performance of the tool during Keck II AO testing on 2016 May 31 UT. The default centroid gain is 0.527 that is correct only when the seeing is 0.75". The operator could change the seeing settings to 0.5", 1.0, 1.5, or 2.0" but this is seldom optimized in practice. We intentionally left the settings at the default value (dashed green lines in Figure 7) to test the performance of the new tool. When we acquire a target on the WFS, the target acquisition tool sets this parameter to the default value and the new automated centroid gain update tool estimates the seeing and updates the WFS centroid gain (blue filled squares) during the science exposure. Also plotted in this figure is the centroid gain estimated through post- processing of AO telemetry data (red filled circles). As the tool has an integrator, high frequency seeing charges are rejected by design. Figure 7 (right) shows an example for the current default state where the AO system uses a wrong centroid gain.

Figure 7:

Left: Automatically updated centroid gain is compared to the default and the ones estimated through post-processing of the AO telemetry during these observations. Right: Same as the left figure but for a night when the centroid gain was not updated.


The automated centroid gain update tool uses the Keck seeing monitor (identified in Figure 1) which is similar to the one presented in Section 3.1 and is an extension of the work of Schöck et al. [16] using the open loop measurements. The tool has been validated using simulations carried out in yao ( over a range of seeing, guide star magnitudes and frame rates (with the appropriate WFS camera program in use). The results show that the seeing estimates are excellent (better than 10% accuracy) and consistent.

Simulations were run in yao to find the correct value of the centroid gain as a function of seeing empirically. The correct value was gauged in two different ways.

First, 300 nm of 90-degree astigmatism was added to the wavefront sensing path. The centroids induced by this aberration were then used to define the centroid origins. We then run simulations with different centroid gains, and save the residual wavefront. The correct centroid gain should lead to a wavefront with no average astigmatism.

Second, the power spectrum of the tip-tilt residual was compared with the modeled power spectrum to see what loop gain it corresponds to. First, the simulations were run using an ideal WFS, which looks at the mean slope of the wavefront over the subapertures (i.e., it does not compute slopes). Then the simulations were run using Fourier optics to simulate the WFS, and the centroid gain that produces the same rejection transfer function is deemed to be correct.

The results from the simulation are presented in Table 1. It can be seen that over the scientifically useful operating range of the AO system (r0 ≥ 7.5 cm), the centroid gain varies by a factor of two.

Table 1:

Centroid gain estimated from the NCPA method and the power spectrum method as a function of Fried’s parameter at 500 nm.

r0 (cm)No turb.25201512.5107.55
From NCPA0.3440.3670.3850.4270.4660.5330.6570.934
From Power spectrum0.3440.3670.3850.4270.4660.5080.6110.865

These results were compared with the expected centroid gain due to the seeing (Eq. (5)) after scaling the seeing for 640 nm from the estimated value at 500 nm and optimizing the fudge factor (FF). The results are plotted in Figure 8. It can be seen that the centroid gain values calculated using the power spectrum and NCPA methods agree very well with each other, especially in good seeing. Where there is a discrepancy, the values found using the NCPA method should be used, since it is more important to reduce the effect of uncorrected non-common-path aberrations than to have the loop gains set correctly to a few percent. The centroid gain using the expression in Eq. (5) with a fudge factor of 0.8 fits the simulated data. This is overdrawn in Figure 8 as a continuous curve.

Figure 8:

Estimated centroid gains from simulations and seeing-based approach (i.e. using Eqn. (5) with a fudge factor of 0.8.)



Telescope segment piston error

There are two major components of static aberrations in the case of the NIRC2 AO instrument. The first one is the astigmatism (~ 250nm rms wavefront error) introduced by the infrared dichroic that transmits the science beam and the aberrations of the NIRC2 camera itself. The second one comes from the optical elements that are prior to the AO bench and hence not probed during daytime AO calibration.

The first component of static aberration is carefully accounted for in the daytime AO calibration procedure by measuring and applying corrections to the WFS centroid offsets. The challenges to this approach are: (1) spot size dependency, (2) WFS operating off-null, and (3) saturation of the DM actuators under poor seeing conditions. As the seeing changes the spot size also changes, resulting in significant leakage of the instrument static aberrations for on-sky observations. The automated centroid gain update tool demonstrated in Section 5.2 addresses this challenge. The centroid gain update is also beneficial to the overall performance of the AO control loop. There is no easy solution to the other two challenges in the current Keck AO system. In spite of these careful AO calibrations, we do see some static component in the science path on the sky that is not seen by the WFS. The second component of static aberration comes from the telescope itself – primarily due to the segment piston error.

Our approach in the past has been to take on-sky long exposure phase diversity measurements to characterize the residual low order static aberrations (Section 2.2). More recently we started questioning this approach. What if segment piston error is the major source of the residual static aberration? We have thought about this during the early part of the project but discorded as the telescope is typically phased to an rms wavefront error of ~ 30 nm. However, there are some serious questions about the stability of telescope segment phasing and the compensation applied to account for the terrace mode.

Here we also attempt measuring high order static aberrations through short exposure NIRC2 image with some diversity using a phase retrieval method. The ability to measure segment piston errors using the science instrument would be beneficial to improving the telescope phasing process as well. Segment piston error could potentially be a significant problem for the segmented extremely large telescopes especially if frequent phasing of the segments is needed.

We have implemented a Gerchberg-Saxton phase retrieval algorithm to validate the performance of the phase diversity algorithm and to explore the possibility of measuring high order aberrations. In-focus images suffer from measurement noise and hence (unlike in the case of phase diversity) we opted to use a single exposure NIRC2 image with some focus diversity for this phase retrieval approach. The approach has its own weakness of difficulties to use with non-point sources for long exposures and speckle noise for short exposures. Nevertheless, we anticipate extracting complementary information using this approach. The test results using the phase retrieval algorithm are presented below for two cases: (1) simulated data with added piston error and (2) on-sky data.


Single exposure phase retrieval: simulations

Simulated NIRC2 images at 1.6455 μm with segment piston errors are shown in Figure 9. The dataset includes one in- focus image and two defocused images (±4 mm), corresponding to an rms wavefront error of 0.48 waves.

Figure 9:

Simulated NIRC2 images with added piston error (150 nm). Left to right: in-focus, 4 mm defocus, and -4 mm defocus.


We applied piston errors to seven of the 36 Keck segments to form the letter “P” with an rms wave front error of 150 nm. The amplitude of the individual piston errors are 0.229, 0.343, -0.343, -0.229, -0.458, 0.229 and 0.458 μm for the segments 5, 8, 15, 18, 22, 32 and 34. In addition, 0.2 radian rms random phase has been added to the telescope pupil. The poked segments can be easily seen from the applied phase maps shown in Figure 10.

Figure 10:

Left: Input phase map. Middle: Retrieved phase map. Right: Projected segment piston error.


The algorithm converges very well and retrieves the piston errors from single exposure defocused image with an image correlation of 99.95%. Figure 10 shows the applied phase map with +4 mm defocus (0.48 waves at 1.6455 μm; left) and the retrieved phase map (middle). Also shown is the projected piston error (right) on to the telescope segments. The results for the -4 mm defocus case are very similar to the case of +4 mm defocus. Even the single exposure in-focus image using the algorithm starts to recover the applied phase, but some sort of diversity is certainly needed, at least for the NIRC2 pixel sampling, to retrieve higher order aberrations such as segment piston error.


On-sky short exposure NIRC2 data

On-sky short-exposure NIRC2 data were taken with ±5 mm defocus offsets (0.6 waves at 1.6455 μm) on 2016 May 31 UT for the purpose of phase retrieval. The retrieved phase map for the +5mm defocus offset is shown in Figure 11 (left). Also shown in this figure is the projected telescope segment piston error (right). The reconstructed phase map shows noticeable segment piston error.

Figure 11:

Left: The retrieved phase map using a modified GS algorithm. Right: projected segment piston error.


It came to light very recently that the compensation for the terrace mode of the Keck II primary mirror is being wrongly applied to the primary mirror control system resulting in a relatively large segment piston error (peak to peak values of ~ as much as 1 μm across the full telescope pupil) for low elevation observations. However, the primary component of this wavefront error is tip/tilt which is not relevant here. While the tip/tilt removed terrace mode has a well-defined saw-tooth pattern, the corresponding rms wavefront error of ~ 70 nm at lower elevations is comparable to that of the typical random segment piston error. Effects are underway to model this error term to improve PSF reconstruction of the

existing engineering data further. The ability to measure the residual segment piston error using AO science instruments is beneficial not only for PSF reconstruction but in providing valuable input for the overall performance of the telescope for science operations.



The on-sky observations for the on-axis NGS case were taken in 2013 with limited supplementary data taken recently (on 2016 May 31 UT), and that for the LGS case were taken in 2014 and 2015. All these measurements were taken with the narrow camera of the NIRC2 imager (pixel scale ~10 mas) using the Fe II narrow band filter (effective wavelength = 1.6455 μm). The measurements include NIRC2 NGS or LGS images taken with a range of guide star magnitudes and NIRC2 NGS in-focus and out-of-focus images for static aberration estimation using phase diversity.

The on-axis results from analysis of the on-sky NGS data taken on 2013 Feb. 03 UT and 2013 Aug. 01 UT, and the LGS data taken on 2015 July 7 UT are presented in this section. The NGS data taken on 2013 Sep. 14 UT is excluded as the static measurements during this engineering run is not satisfactory. Also, the measurements taken at 1054 Hz are excluded as the internal test show significant noise in the AO telemetry at framerates (Section 5.2.)

Figure 12 compares the Strehl ratio and the FWHM of the reconstructed and the NIRC2 PSFs. The contribution from the diffraction pattern is quadratically subtracted from the reconstructed and sky FWHM values. Overall there is a good match of the reconstructed and the NIRC2 PSFs in terms of Strehl ratio and FWHM. The correlation coefficient for the Strehl ratio and FWHM comparisons is 0.97 and 0.99 respectively. The rms value of the percentage difference between the reconstructed and the sky FWHM is ~10%. Overall, the reconstructed FWHM is ~ 15% larger than that of the sky value. When the fixed offset (~ 3 mas) is removed, the percentage difference goes down to ~ 5%. In the case of Strehl ratio, the rms value of the percentage difference between the reconstructed and the sky values is ~12%. The high Strehl ratio measurements, taken at relatively high frame rates, show some systematic bias that is being investigated.

Figure 12:

Left: The reconstructed Strehl ratio is compared with the NIRC2 Strehl ratio for the 2013 Feb. 03 UT and 2013 Aug. 01 UT NIRC2 data. Right: Same as the left plot but for the FWHM after quadratically subtracting the contribution from the diffraction pattern.


Figure 13 & Figure 14 compare reconstructed PSF with the Sky PSF for a bright high Strehl ratio (Figure 13) and a faint low Strehl ratio (Figure 14). The normalized reconstructed PSFs compare very well with the sky PSFs for the two cases.

Figure 13:

The reconstructed NGS PSF (middle) is compared with the sky PSF (left) for a bright high Strehl ratio case. The images were taken at 1.655 μm. Right: a scan across the PSF through the peak is shown in logarithmic scale.


Figure 14:

Left: The reconstructed NGS PSFs are compared with the sky PSFs for a faint low Strehl ratio case. The images were taken at 1.655 μm. Right: a scan across the PSF through the peak is shown in logarithmic scale.


The reconstructed LGS A PSFs taken on 2015 July 7 UT are compared with the sky PSFs in Figure 15: (left) and the 1D intensity profiles is shown in Figure 15 (right). As in the case of the NGS observations, the sky PSFs were taken with the NIRC2 camera with Fe II filter (λeff = 1.6455 μm; δλ= 0.0256 μm). The tip-tilt star was on-axis for these observations. Again, the reconstructed PSFs match fairly well with the sky PSFs. Also shown in this figure (right; dotted line) are the reconstructed profiles ignoring the focal anisoplanatism.

Figure 15:

Left: The reconstructed LGS PSFs are compared with the sky PSFs for a typical Strehl ratio and a low Strehl ratio case. The images were taken at 1.655 μm. Right: a scan across the PSF through the peak is shown in logarithmic scale. The cross symbol, the solid line, and the dotted line represent the sky profile, the model profile and the model profile ignoring focal anisoplanatism.




Off-axis PSF prediction tools are developed at UCLA. The Galactic Center Group @ UCLA is about to finish a project to optimize the extraction of AO astrometry and photometry by incorporating measurements of the Earth’s atmospheric turbulence profile and instrumental aberrations to build a model of the PSF spatial variations in AO image obtained with the NIRC2 instrument on the Keck II telescope.

The UCLA team has developed algorithms to predict how the PSF varies and incorporated them into a new software package: AIROPA (Anisoplanatic and Instrumental Reconstruction of Off-axis PSFs for AO). This software package makes use of independently measured atmospheric data like MASS and DIMM profiles and predicts the differential OTFs between different field positions. Based on ARROYO (M. Britton 2006), a set of C++ class libraries that aim to support simulations of electromagnetic wave propagation through turbulence and through optical systems, AIROPA includes C++ tools for predicting both natural guide star and laser guide star PSFs, as well as application program interfaces (APIs) for IDL. Furthermore, it provides extensive IDL modules for Fourier-based PSF manipulation and for PSF modeling from any externally generated aberration maps. These IDL routines and APIs are integrated with a modified version of StarFinder [17].

A flow chart of AIROPA’s software structure is shown in Figure 16. Based on additional data describing the atmospheric turbulence at the time of observation, a set of spatially variable PSFs is modeled. AIROPA is a semi- empirical approach that extracts information on the on-axis PSF from the science image itself (PSF Extraction). In the future we plan to replace this part with the on-axis model generated from AO telemetry data.

Figure 16:

A flow chart of the off-axis software showing various components.


The AIROPA software, a modified version of StarFinder (StarFinder_AIROPA), and related IDL interface routines have been successfully installed at WMKO to be used for other off-axis science programs as well. With these steps we have entered the phase of on-sky tests of our software package, and we expect to have a validated off-axis prediction tool by end of 2016.



Astronomers measure fundamental properties like mass, density, luminosity and structure via quantitative measurements of positions, brightness, morphology and kinematics to understand the physical properties of the universe. The accuracy of the quantitative science measurements is limited by the accuracy to which the PSF is known and by its variability. A highly accurate and detailed 2-dimensional (2D) PSF provided with each science observation would be ideal. In practice this will be very challenging and in many cases may not be needed. It is therefore important to understand what features of the PSF have the largest impact on different types of science measurements and results.

The type of PSF knowledge that would provide the most benefit to a particular science case is highly specific. Overall a careful characterization of the PSF would benefit both observing efficiency and measured uncertainties. Table 2 shows the science cases identified for the initial science verification.

Table 2:

Science Cases.

Research AreaObservablePSF parametersObserving ModeScience collaborators
Dynamical mass measurements for Brown Dwarfs and Low-Mass StarsAstrometryFWHM(near) On-axis NGS/LGSM. Liu & T. Dupuy
Planet searchAstrometryFWHM & Strehl ratioOn-axis NGSS. Ragland & B. Bowler
Galaxy Morphology and KinematicsMorphology & KinematicsFWHM & Strehl ratioOn-axis/off- axis LGSS. Wright
Stellar population in the Galactic centerAstrometry, Photometry, Morphology & Kinematics2D profileOff-axis LGSA. Ghez
Gravitational lensingAstrometry, Photometry, Morphology & Kinematics2D profileoff-axis LGST. Treu

We have collected initial science data with and without a coronographic mask for a planet search program (PI: S. Ragland) and more recently the target list for this case has been extended through a collaboration with B. Bowler and M. Liu. We plan on testing the algorithm on a set of bright binary stars in the NGS configuration first before testing the other NGS and LGS on-axis and off-axis cases.



The on-axis PSF reconstruction algorithm development is complete and the code is being tweaked to improve the performance. The reconstructed PSFs are compared with the sky PSFs in terms of FWHM and Strehl ratio. High level prototype operational tools to reconstruct on-sky NIRC2 PSFs have been developed and some initial tests were carried out.

We find that accurate estimation of (1) atmospheric seeing (< 10% accuracy), (2) non-common path aberration, and (3) telescope segment piston error are essential for PSF-R. A well-optimized AO system with known control loop characteristics is also crucial for PSF-R. The ability to measure residual telescope piston error using the AO science instrument is beneficial for PSF-R as well as for the overall performance of the telescope for science operations.

The next steps are:

  • (1) Document the PSF reconstruction process in a series of journal articles.

  • (2) Integrate off-axis components developed at UCLA and validate the algorithm with on-sky engineering data.

  • (3) Conduct an operational tool design review.

  • (4) Carryout science verification.


The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

This material is based upon work supported by the National Science Foundation under grant AST-1207631, the W. M. Keck Foundation, and the National Aeronautics and Space Administration (NASA) for KOA support.

The authors wish to acknowledge our technical/scientific collaborators at UCLA (Max Service and Tommaso Treu), tOSC (Matthew Britton), UCSD (Shelley Wright), and IFA (Michael Liu), UT Austin (Trent Dupuy and Brendan Bowler), TMT (Luc Gilles and Lianqi Wang) and LAM (Carlos Correia and Olivier Martin). Special thanks to Kuochou Tai for telescope phasing discussions, and Hien Tran and Bob Goodrich for their support for data storage.

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.



Jolissaint, L., Ragland, S., Wizinowich, P., “Adaptive Optics Point Spread Function Reconstruction at W.M. Keck Observatory in Laser & Natural Guide Star Modes: Final Developments,” in Proc. AO4ELT4, (2016). Google Scholar


Ragland, S., Jolissaint, L., Wizinowich, P., Witzel, G., Chock, J., Mader, J., Morris, M. R., and Kwok, S., “Point spread function determination for Keck adaptive optics: overview,” in Proc. AO4ELT4, (2016). Google Scholar


J.-P. Véran, F. Rigaut, H. Maître, and D. Rouan, “Estimation of the Adaptive Optics Long-exposure Point Spread Function using Control Loop Data,” Journal of the Optical Society of America A, 14 3057 –3069 (1997). Google Scholar


Jolissaint, L., Ragland, S., Wizinowich, P. & Bouxin, A., “Laser guide star adaptive optics point spread function reconstruction project at W. M. Keck Observatory: preliminary on-sky results,” in Proc. SPIE, 4SJ (2014). Google Scholar


Britton, M., “Arroyo,” in Proc. SPIE, 290 –300 (2004). Google Scholar


Fitzgerald, M.P., Witzel, G., Britton, M.C., Ghez, A.M., Meyer, L., “Modeling anisoplanatism in the Keck II laser guide star AO system,” in Proc. SPIE, 24 (2012). Google Scholar


Sitarski, B., “Modeling instrumental field-dependent aberrations in the NIRC2 instrument on the Keck II telescope,” in Proc. SPIE, 9148 –265 (2014). Google Scholar


Jolissaint, L. and Blanc, Ph., “Zernike polynomials, annular pupils and optical turbulence,” submitted to JOSA A, (2016). Google Scholar


Jolissaint, L., Christou, J., Ragland, S., and Wizinowich, P., “Optical turbulence characterization using adaptive optics deformable mirror telemetry,” submitted to JOSA A, (2016). Google Scholar


Grey, L.D. and Lyon, R.G., “Correction of misalignment dependent aberrations of the Hubble Space Telescope via phase retrieval,” 201 (1989). Google Scholar


Fienup, J. R., Marron, J. C., Schulz, T. J., and Seldin, J. H., “Hubble Space Telescope characterized by using phase- retrieval algorithms,” ApOpt, 32 1747 (1993). Google Scholar


Krist, John E., and Burrows, Christopher J., “Phase-retrieval analysis of pre-and post-repair Hubble Space Telescope images,” ApOpt, 34 4951 (1995). Google Scholar


Mugnier, L. M.; Sauvage, J.-F.; Fusco, T.; Cornia, A.; Dandy, S., “On-Line Long-Exposure Phase Diversity: a Powerful Tool for Sensing Quasi-Static Aberrations of Extreme Adaptive Optics Imaging Systems,” Optics Express, 16 18406 (2008). Google Scholar


Berriman, G. B., Gelino, C.R., Goodrich, R.W., Holt, J. Kong, Mihseh; Laity, A.C.; Mader, J.A.; Swain, M.; Tran, H.D., “The design and operation of the Keck Observatory archive,” in SPIE Proc., 0AB (2014). Google Scholar


van Dam, M. A., Le Mignant, D., and Macintosh, B.A., “Performance of the Keck Observatory adaptive-optics system,” Applied Optics, 43 5458 (2004). Google Scholar


Schöck, M., “Atmospheric turbulence characterization with the Keck adaptive optics systems. I. Open-loop data,” Applied optics, 42 3705 –3720 (2003). Google Scholar


Witzel, G., “The AIROPA software package: milestones for testing general relativity in the strong gravity regime with adaptive optics,” 9909 –62 (20116). Google Scholar
© (2016) COPYRIGHT Society of Photo-Optical Instrumentation Engineers (SPIE). Downloading of the abstract is permitted for personal use only.
S. Ragland, L. Jolissaint, P. Wizinowich, M. A. van Dam, L. Mugnier, A. Bouxin, J. Chock, S. Kwok, J. Mader, G. Witzel, Tuan Do, M. Fitzgerald, A. Ghez, J. Lu, G. Martinez, M. R. Morris, and B. Sitarski "Point spread function determination for Keck adaptive optics", Proc. SPIE 9909, Adaptive Optics Systems V, 99091P (26 July 2016);

Back to Top