博士論文

Multi time-step tomographic
wave-front reconstruction for wide
field Adaptive-Optics in astronomy(天文学における広視野補償光学のための複タイムス
テップトモグラフィー波面再構成)

大野 良人
平成 27 年

Abstract
In this thesis, we develop and evaluate a new tomographic wave-front reconstruction method for wide-field adaptive optics (WFAO) systems, called multi time-stepwave-front reconstruction, to solve the lack of information in tomographic reconstructions and expand the field of regard (FoR) of WFAO systems. Also, weevaluate a method to estimate wind speeds and directions at multiple altitudes
using the measurements from multiple wave-front sensors. This wind estimation
method is required to implement the multi time-step reconstruction. The basic
idea of the multi time-step reconstruction is to increase the information using themeasurements at both the current and previous time-steps, which requires windspeeds and directions of wind and the frozen flow assumption.
First, we demonstrate the influence of the lack of information in the tomographic reconstruction. Then, we present that the multi time-step tomographic
reconstruction method is eﬀective to reduce the tomographic error due to the lackof information.
Second, we show the result of the numerical simulation of the tomographic
reconstruction on a multi-object adaptive optics (MOAO) system on a future
extreme large telescope. The numerical simulation shows that the multi time-stepreconstruction increases Strehl ratio (SR) over a field of regard of 10 arcminutes
diameter by a factor of 1.5–1.8 if we know the wind speeds and directions.
Third, we evaluate the multi time-step reconstruction and the wind estimation
method in the laboratory experiment on RAVEN, which is a MOAO demonstrator.
We can successfully measure the wind speeds and directions in the laboratory
experiment. Also, the multi time-step reconstructor can increase an ensquaredenergy (EE) in a 140 mas box by 0.03–0.05.
Finally, we present the on-sky performance of the multi time-step reconstruction and the wind estimation with RAVEN on the Subaru telescope. Unfortunately, there is no clear improvement by the multi time-step reconstruction because of the error of the wind estimation. With the oﬀ-line analysis with theon-sky measurements of the wave-front sensors, we improve the wind estimation
method and re-evaluate the multi time-step reconstruction. We conclude that themulti time-step reconstruction can reduce the WFE compared with the normal tomographic reconstruction if the decay ratio of the cross-correlation of WFS slopes
is more than 0.7 and there is no dome seeing.

iAcknowledgements
I would like to express my special appreciation and thanks to my supervisor
Masayuki Akiyama. He always encouraged me in the research and give me precious opportunities to research with overseas researchers. Your advice on bothresearch as well as on my career have been priceless. I would also like to thankShin Oya, who supports my research as well as my stay in Hawaii. I would also
like to thank Olivier Lardiére, Carlos Correia, Colin Bradley, and other RAVEN
members. This thesis is based on the work with them in the RAVEN project, andthey support my research. I would also like to thank staﬀ members of Subaru
telescope for their support.
I would like to thank the three referees of my thesis, Hideki Takami, Takashi
Murayama, Takashi Ichikawa. I received many helpful comments from them. I
would like to thank the students in the astronomical institute of Tohoku University.
A special thanks to my family. At the end I would like express appreciation to
my wife Aya who always support me.

ii

Contents
Abstracti

Acknowledgements

ii

Contents

iii

List of Figuresv

List of Tablesxii

1 Introduction to Wide-Field AO Systems and Tomographic Reconstruction
1
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Atmospheric Turbulence . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Principle of SCAO system . . . . . . . . . . . . . . . . . . . . . . . 8
1.4 Wide-Field Adaptive Optics . . . . . . . . . . . . . . . . . . . . . . 14
1.5 Tomographic Wave-Front Reconstruction . . . . . . . . . . . . . . . 21
1.6 Object of This Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2 Methodology
27
2.1 Tomographic Wave-Front Reconstruction . . . . . . . . . . . . . . . 27
2.2 Pre-Estimation Method for Required Parameters . . . . . . . . . . . 36
3 Analytical Evaluation
42
3.1 Analytical Tomographic Error . . . . . . . . . . . . . . . . . . . . . 42
3.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.3 Eﬀect of Lack of Information . . . . . . . . . . . . . . . . . . . . . . 45
3.4 Analytical Performance of Multi Time-Step Reconstruction . . . . . 47
3.5 Dependency on Time Diﬀerence ∆t . . . . . . . . . . . . . . . . . . 48
4 Numerical Simulation
51
4.1 Simulation Setting . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

iii

5 Laboratory Test with RAVEN
69
5.1 RAVEN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2 SLODAR and Wind Estimation . . . . . . . . . . . . . . . . . . . . 71
5.3 Multi Time-Step Reconstruction . . . . . . . . . . . . . . . . . . . . 75
5.4 Comparison with a numerical simulation . . . . . . . . . . . . . . . 76
6 On-Sky Evaluation with RAVEN
82
6.1 On-Sky Performance . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.2 Oﬀ-Line Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
7 Conclusion

106

A Mathematics
108
A.1 Minimum Variance Reconstructor and Analytical Tomographic Error108
A.2 The Scaling of Discrete Laplacian Matrix . . . . . . . . . . . . . . . 109
Reference

117

iv

List of Figures
1.1

1.2

1.3
1.4
1.5
1.6

The power spectrum (the left panel) and the structure function(the right panel) for phase fluctuations with the Kolmogorov model(the red solid lines) and the von Karman model as L0 =100 m (theblue dashed lines) and L0 =30 m (the green dotted lines). Fried
parameter is set to r0 =0.2 m. . . . . . . . . . . . . . . . . . . . . . 7
The seeing size for the Kolmogorov model (the red solid lines) from
Eq. (1.15) and the von Karman model as L0 =100 m (the blue dashedlines) and L0 =30 m (the green dotted lines) computed by Eq. (1.16).
Fried parameter is set to r0 =0.2 m at λ=0.5 µm. . . . . . . . . . . 8
The schematic view of a SCAO system. . . . . . . . . . . . . . . . . 9
The schematic images for the cone eﬀect and the spot elongation
due to the LGS. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
The schematic views of WFAO systems. . . . . . . . . . . . . . . . 15
Observed number density and the number density achieved by themulti-object spectroscopy with the future MOAO system at z ∼ 5(the left panel) and z ∼ 6 (the right panel). The gray curve indicates the observed cumulative luminosity function of high redshift
LBGs as a function of an apparent magnitude at J-band. The
observed luminosity functions are modeled with a Schechter function [7], and the parameters of the Schechter function are MU∗ V =−21.17, ϕ∗ = 0.74 × 10−3 , and α = −1.76 for z ∼ 5, and MU∗ V =−20.94, ϕ∗ = 0.50 × 10−3 , and α = −1.87 for z ∼ 6. The redshift
range for each redshift is assumed to 4.5 < z < 5.5 for z ∼ 5 LBGs
and 5.5 < z < 6.5 for z ∼ 6. The vertical lines indicate the limit
magnitude achieved by the future instrument for the point source(solid) and extended source (dotted) as an integration time of 10hours (red), 20 hours(blue), and 30 hours (green). The horizontal
lines represent the cumulative number density required to detect 10(solid) and 5 (dotted) LBGs over the FoR of 5 (magenta) and 10(black) arcminute diameter. . . . . . . . . . . . . . . . . . . . . . . 20

v

1.7

Geometric relation between positions of GSs and regions in theatmosphere measured by WFSs. Left: The solid orange circles show
the footprint of GS optical paths at the current time-step. There areareas covered by no or only one GS, referred to as uncovered areaor unoverlapped area. These areas cause significant tomographic
errors. Right: The blue arrows indicate wind direction at each
altitude. The GS footprints measured at a previous time-step move
due to the wind with time, which are dashed orange circles in the
figure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

2.1

Image of a MOAO model with multiple turbulence screens, LGSs,
and SH-WFSs. β and α indicate directions of science targets and
LGSs. The index i, j, k is used for the turbulence screens, LGSs
and SH-WFSs, and science targets and DMs, respectively. The goalof the tomographic reconstruction is to estimate DM command uβ
from the measured slope sα . . . . . . . . . . . . . . . . . . . . . . .
Fried geometry to define slopes for SH-WFS. . . . . . . . . . . . . .
Principle of SLODAR. The lights from two star with an angularseparation θ cross at altitude h with a spatial displacement of hθand propagate into [i, j] subaperture of WFS1 and [i + δi, æ; δj]subaperture of WFS2, respectively. The cross-correlation between
these subapertures have a correlation power proportional to theturbulence power at altitude h. . . . . . . . . . . . . . . . . . . . .
Examples of theoretical slope cross-correlation. The outer scale is
30 m. The only areas surrounding by the red lines are used for thefitting. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2
2.3

2.4

3.1

3.2

28
30

36

38

Footprints of the NGS optical paths (red) and the evaluated central
direction (blue) at altitude of 10 km for rast =20 arcseconds (left)
and rast =200 arcseconds (right). . . . . . . . . . . . . . . . . . . . . 44
The NGS asterism (red) and the evaluated directions. The blue plotis the central evaluated direction. The green plots are the outside
evaluated directions. . . . . . . . . . . . . . . . . . . . . . . . . . . 44

vi

3.3

3.4

4.1

4.2

4.3

The analytical geomatric error σgeo (top), the noise propagation
σnoise (middle), and the total tomgoraphic error σtomo (bottom) as
a function of the asterism radius rast . The all values are computedwith removing tip-tilt modes from the residual wave-fronts. The
red lines show the result with the single time-step tomographic reconstruction and the blue lines are the results of the multi time-stepreconstruction. The solid lines indicate the value evaluated in the
direction of the center of the FoR. The dashed lines are the results
averaged over outer directions. . . . . . . . . . . . . . . . . . . . . . 46
The analytical geomatric error σgeo (top), the noise propagation
σnoise (middle), and the total tomgoraphic error σtomo (bottom) as
a function of the asterism radius rast . The all values are computedwith removing tip-tilt modes from the residual wave-fronts. The
color and line type are same as Eq. (3.3) . . . . . . . . . . . . . . . 50
LGS asterisms used in the numerical simulation. A narrow asterism
is indicated with the red filled squares. The blue open squares showa wide asterism. The black solid lines show the 2.5 arcminutes and5 arcminutes radius from the center of a FoR. The dashed lines are
directions where the performance is evaluated. . . . . . . . . . . . . 52
Vertical profile of normalized intensity of Sodium layer used in thesimulation. This is approximated based on the lidar on-sky measurements [19]. The altitude and of sodium layer are 88 km and 10
km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Top and middle : Simulated SR maps within the 10 arcminutes
FoR computed by using the single time-step reconstructor (top) or
the multi time-step reconstructor with ∆t=50 ms (middle). The
right panels are the results with the narrow LGS asterism, and theleft panels are computed with the wide LGS asterism. The multitime-step reconstructor improves SR values over the FoR for bothof the asterisms. Bottom : SR profiles as a function of an angularseparation from the center direction of the FoR. The profiles areaveraged over directions shown as the dashed lines in Fig. 4.1. The
red filled squares show the averaged SR profile with the single timestep reconstructor, and the blue open squares is the result by themulti time-step reconstructor with ∆t=50 ms. The dashed blacklines show the radii of LGS positions. . . . . . . . . . . . . . . . . . 55

vii

4.4

4.5

4.6

4.7

5.1

The footprints of the LGS optial paths on the highest turbulence
screen at altitude of 16 km The left panal shows the footprints with
the narrow asterism and the right is for the wide asterism. The
dashed black circle indicates the radius of 150 arcseconds and 300
arcseconds from the center of the FoR. . . . . . . . . . . . . . . . .
Dependence of the multi time-step reconstructor on the ∆t is presented in Fig. 4.5. The vertical axis in Fig. 4.5 is the SR improvement ratio, kSR , which represents the ratio of the SR value achievedby the multi time-step reconstructor to the SR value from the single
time-step reconstructor and is averaged over the angular separationfrom the center of the FoV. The result with the narrow asterism isindicated as the red filled squares, and the blue open squares shows
the result with the wide asterism. The vertical dashed lines indicate the time duration that the decay ratio fdecay of the temporalcorrelation of SH-WFS measurements becomes 90 %, 80 %, and 70%, that is, t90 , t80 , and t70 , wind assuming the wind speed is 33 m
s−1 , which is maximum value in the wind model. . . . . . . . . . . .
SR improvement ratio, kSR , achieved by multi time-step tomographic reconstruction with diﬀerent wind speed errors e∥ (the filled
symbols) or wind direction errors e⊥ (the open symbols). The results with the narrow asterism are represented as the red squares
and the results with the wide asterism are the blue circles. The SR
ratio less than 1 means that the performance of the multi time-stepreconstructor is poorer than the performance achieved by the single
reconstructor due to the wind estimation error. . . . . . . . . . . .
Ensqured energy profiles with the narrow (left) and wide (right)
asterism. As the target intensity profile, the point source and extended source are considered. The extended source is assumed to
be an exponential profile with an eﬀective radius of 0.1 arcsecond.
The EE values are measured from the simulated PSF images (red)with a 50 mas box for the point source, and the exponential profile
convolved with the simulated PSF (blue) with with a 300 mas box
for the extended source. The solid lines indicate the results withthe single time-step reconstructor, and the dotted lines are the EE
profiles of the multi time-step reconstruction with ∆t=50 ms. . . .

56

58

59

60

Block diagram of the RAVEN system. Figure adopted from Lardière
et al. (2014) [32]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
viii

5.2
5.3
5.4
5.5

5.6

6.1

6.2

Expected CN2 fraction generated by the RAVEN CU. . . . . . . . .
Asterism of the GSs and the science targets for the laboratory testwith the RAVEN. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Profile of atmospheric turbulence of the RAVEN CU estimated by
the SLODAR. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Temporal correlation maps of slope at each altitude estimated bytomographic reconstruction. The values at the top of each image
indicate the time delay used to drive the temporal correlation. The
x and y-axes are δi and δj with a range of −9 ≤ δi, δj ≤ 9. One
pixel on the temporal-correlation corresponds to one subaperturesize, which is 0.8 m in the RAVEN . . . . . . . . . . . . . . . . . .
PSF images of each science channel taken in the lab-test with thesingle and the multi time-step reconstructors. The color scale islinear and aligned for each channel. The wavelength is 1650 nm(H-band). The size of each image is 0.5 arcsecond × 0.5 arcsecond,
and the pixel scale is 17.5 mas. . . . . . . . . . . . . . . . . . . . .

71
72
73

74

78

Images of the engineering fields. The red lines indicate the radiusof 1.75 arcminute and 1 arcminute from the center of the field.
The green points represent the NGSs position. The brightnessesof the NGSs are summarized in Table 6.1. The information of theengineering field is summarized in Table 6.1. These images comesfrom the Digital Sky Survey. The position of NGSs are plottedbased on the USNO-A2.0 catalog. . . . . . . . . . . . . . . . . . . . 83
The example of the SLODAR result dominated with a grund layer.
The upper left panel shows the NGS asterism and the dashed circle
indicates the angular separation of 60 arcseconds from the centerof the FoR. The upper right panel shows the CN2 profile estimatedby the correlation SLODAR and the bottom panel shows the corresponding correlation maps. The red arrows indicate the baselineof the guide stars, and the signals of the correlation are expected
to apper along these arrows. In this case, most of the turbulence
is concentrated in the ground layer and there are weak turbulencesat around 12 km. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

ix

6.3

6.4

6.5

6.6
6.7

6.8
6.9

Another example of the SLODAR result with a complex profile.
The turbulences occur over wide range of altitudes and the turbulences at altitude higher than 10 km have significant contribution
to the total turbulence power. On the correlation maps, there arestrong correlations along the base lines of the NGS pairs . . . . . .
Histogram of ground layer fraction measured by the covariance
SLODAR during on-sky observations from June 23th to July 2ndin 2015 (HST). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Top: Temporal-correlation maps with diﬀerent time delays, δt, ateach altitudes for the cases of Fig. 6.2. The x and y-axes are δi andδj with a range of −9 ≤ δi, δj ≤ 9. One pixel on the temporalcorrelation corresponds to one subaperture size, which is 0.8 m in
the RAVEN. The temporal-correlation maps are normalized by thepeak intensity of the map of δt=0 at each altitude. The color scale ofthe maps is set to [0 0.7]. Bottom: Wind speed (red) and direction
profile (blue) as a function of altitude from the telescope. The windspeed follows the left y-axis, and the wind direction follows the righty-axis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The same figure as Fig. 6.5 with the CN2 profile shown in Fig. 6.3. . .
The (a) FWHM, (b) SR, and (c) EE in 140 mas box measured fromthe PSF images obtained in the on-sky engineering observations.
These metrics are plotted as a function of a r0 . The red squresindicate the results of the single time-step reconstructor and theblue circles are the result of the multi time-step reconstructor. The
filled and open symbols indeciate the result of the science channel
1 and 2, respectively. In the enginnering observation, we set ∆t
to a duretion corresponding to 20 frames. In the most case, theframerate of the control for MOAO system is 150 Hz and ∆t is 133
ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The on-sky PSF images corrected by the single and multi time-step
recostructor. The size of PSF images is 1.5×1.5 arcsecond. . . . . .
The profiles of CN2 values (left panel), wind speeds (right panel),and wind directions (right panel) measured on sky. The field is
Fd38 and the asterism is same as Fig. 6.2. In the CN2 profile, morethan 80 % of the turbulence are concentrated in the ground layer
and there are weak turbulences at around 10 km. On the left panel,the detected wind speeds are almost zero at all altitudes. . . . . . .

x

86

88

89
90

91
93

93

6.10 The improved version of Fig. 6.5 . . . . . . . . . . . . . . . . . . . . 96
6.11 The improved version of Fig. 6.6 . . . . . . . . . . . . . . . . . . . . 97
6.12 Example of the phase distortion pattern averaged over each one
second at diﬀerent altitudes. . . . . . . . . . . . . . . . . . . . . . . 98
6.13 r0 profiles in the zenith direction on June 23–24th, 2015. The blue
plots indicate the r0 values measure by the correlation SLODAR
of RAVEN. The red plots are the r0 values measured by the autocorrelation fitting of RAVEN. The gray lines are r0 values measured
by the CFHT DIMM (http://mkwc.ifa.hawaii.edu/). . . . . . . . . 98
6.14 CN2 (h) profile in the zenith direction on June 23–24th, 2015. The
y-axis is altitudes from the telescope. . . . . . . . . . . . . . . . . . 99
6.15 Wind speed profile on June 23–24th, 2015. The y-axis is altitudes
from the telescope. . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
6.16 Wind direction profile on June 23–24th, 2015. The y-axis is altitudes from the telescope. . . . . . . . . . . . . . . . . . . . . . . . . 100
6.17 Asterism (top left), CN2 profile (top right), the wind profile (middle), and the tip-tilt removed residual WFE (bottom). The residual
WFEs are computed from the on-sky measurement of the SH-WFSs
with the single and multi time-step reconstruction. The multi timestep reconstructor is computed with the revised wind speeds anddirections and diﬀerent time delay. . . . . . . . . . . . . . . . . . . 103
6.18 The same figure as the Fig. 6.17 with a diﬀerent case. . . . . . . . . 104
6.19 Decay ratio for the ground layer in the first case shown in Fig.
6.17. The red and blue lines indicate the computed decay ratio
and a linear line for the moving component determined by the leastsquares method, respectively. . . . . . . . . . . . . . . . . . . . . . 105
6.20 Decay ratio for the ground layer in the first case shown in Fig. 6.18.
The red, blue and green lines indicate the computed decay ratio,a linear line for the moving turbulence, and a linear line for thedome seeing, respectively. The linear lines are computed by theleast squares method with the time delay range of [0 0.6] second forthe moving turbulence and [0.8 1.5] second for the dome seeing. . . 105
A.1 The. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

xiList of Tables
1.1

Parameters used in the estimation of the number of detected highredshift LBGs in 1.4.4. The dark count and readout noise is referenced from the TMT IRIS simulation [51]. . . . . . . . . . . . . . . 19

3.1

Parameters used for the analytical evaluation. . . . . . . . . . . . . 43

4.1
4.2
4.3

Parameters used in the numerical simulation. . . . . . . . . . . . . .
Atmospheric turbulence profile used in the numerical simulation. . .
Detected number of high redshift LBGs at z ∼ 5 and 6. These
are computed with the observed luminosity function from [7], thesimulated EE profile shown in Fig. 4.7, and the parameters used
in 1.4.4. The values in and out of the parenthesis indicates thedetected number with the narrow and wide asterism, respectively. .
The three rows from the top show the number of WFS measurements, phase points of the atmospheric turbulence layers, and thephase points on the aperture plane. The assumed parameters arelisted in Table 4.1. Other rows indicate the number of total elements, non-zero elements, and the number ratio of the non-zeroelements to the total elements of the matrices of Γα , Pα , and L. . .
Computation flow of the CGM. . . . . . . . . . . . . . . . . . . . .
Computation flow in the loop of the CGM. . . . . . . . . . . . . . .

4.4

4.5
4.6

xii

54
54

61

62
64
64

4.7

4.8
5.1
5.2
5.3
5.4

5.5

6.1

Computational time and the number of iteration for the convergenceof the CGM. The unit for the computational time is µs. While thevalue out of the parenthesis at the processes of 10 and 11 are assumed to compute the wavefront for 20 science directions by one
computer, the values in the parenthesis at these processes are assumed to compute the wavefront for 20 science directions by 20
computers, i.e, it is the computeion time for one science target direction. The tital computetinal times follow the same manner as
the values at the processes of 10 and 11. The number of iterationis computed with and without the diagonal scaling preconditioner,
where the values in the parenthesis are with the diagonal scaling
preconditioner. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Computational time of the processes in the loop of the CGM . . . . 68
Parameters of RAVEN. . . . . . . . . . . . . . . . . . . . . . . . . .
Atmospheric turbulence profile and wind profile generated by the
RAVEN CU. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Atmospheric turbulence profile and wind profile estimated by the
SLODAR and the wind estimation method. . . . . . . . . . . . . .
SR, kSR , and EE value of each science channel measured in thelab-test (the top table) compared to those predicted from the numerical simulation (the bottom table). The EE is defined with a
140 max box. There are two GS configuration used in the test, one
is only three NGSs and the other is three NGSs+LGS. The values in parenthesis in the lower part is the value accounting for the
implementation errors. . . . . . . . . . . . . . . . . . . . . . . . . .
SR, kSR , and EE value of each channel predicted from the numericalsimulation with taking account of the additional tip-tilt errors. The
additional tip-tilt errors are determined by the measurement on thelaboratory experiments. The values in a parenthesis in Table 5.5 isthe SR including the eﬀect of the implementation errors. . . . . . .

70
70
72

77

81

RA, Dec, and NGS magnitudes of the engineering field shown in
Fig. 6.1. The NGSs magnitude are referenced from the USNO-A2.0
catalog. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

xiii

Chapter 1
Introduction to Wide-Field AO
Systems and Tomographic
Reconstruction
1.1

BackgroundIn ground-based observations for astronomy, we suﬀer from the turbulence in the
Earth’s atmosphere, especially in the optical and near infrared (NIR) wavelengths.
Due to the atmospheric turbulence, the phase of the light propagating from space
to a telescope through the atmosphere is disturbed and the image produced by thetelescope is blurred. As a result, we lose the sensitivity and the spatial resolution
for the astronomical observation. Adaptive optics (AO) is a system which corrects the phase distortion due to the atmospheric turbulence. Now, most of large
ground-based telescopes, which have a primary mirror of 8–10 meters diameter,
have a classical AO system, called single-conjugate adaptive optics (SCAO). TheseSCAO systems successfully achieve the high spatial resolution close to the diﬀraction limit of these telescopes and allow us to observe fainter objects and smaller
structures in galaxies. AO systems are considered as an essential technique for
future extreme large telescopes (ELTs), which have a primary mirror of diameter
larger than 30 meters, to maximize their spatial resolution and sensitivity.
Classical SCAO systems have a critical limitation. It is a limitation that thecorrection of SCAO is eﬀective only within a limited angle from a guide star (GS).
Wide-field AO (WFAO) concepts, involving tomographic wave-front reconstruction to estimate the atmospheric volume, come out from requirements to apply
AO correction into a wide field and to increase the eﬃciency of AO assisted observations. The developments for WFAO system have already started and actually,some WFAO systems are working on the current large telescopes as a facility
1Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

2

instrument or demonstrator. However, there are many new techniques and diﬃculties for realizing WFAO systems on the future ELTs. One of the new techniques
not used in SCAO systems is tomographic wave-front reconstruction. In this thesis, we develop a tomographic wave-front reconstruction algorithm to expand thefield in which the WFAO correction is eﬀective, called field of regard (FoR), andto maximize the eﬃciency of WFAO systems.
We will start this chapter with explaining a basic theory of the atmosphericturbulence. Since the atmospheric turbulence is a non-stationary random phenomenon, we need to understand and characterize the atmospheric turbulence
statistically. This statistical theory of the atmospheric turbulence is an importantbasis to design and evaluate AO systems. We will then present the principle of
SCAO systems and their limitations. Then, we will outline three WFAO concepts
and the challenging points to realize these WFAO systems especially for future
ELTs. We will provide more introduction about the tomographic wave-front reconstruction. At the end of this chapter, we will summarize the object of thisthesis.

1.2

Atmospheric TurbulenceThis section is based on Hardy (1998) [24].

1.2.1

Kolmogorov Theory

Kolmogorov[30] investigated the statistical structure of turbulence, which is thevelocity fluctuations in a fluid medium. In his theory, when energy is injected to
the fluid medium, it forms a large length-scale disturbance. The length-scaleat which the energy injection occurs is referred to as outer scale, L0 . Then,the injected energy is broken down into smaller length-scale, forming an energy
cascade, until a small scale where the energy is dissipated by viscosity, called innerscale ℓ0 .
In Earth’s atmosphere, the main energy input is solar heating during the daytime, which causes turbulences on the length-scale of L0 by the local convection.
In the night time, the air masses of diﬀerent altitudes and temperatures are mixed
by wind shearing and the turbulences break down into smaller scale until the innerscale ℓ0 . This process causes temperature fluctuations in the atmosphere, which affect the density of the air and therefore result in refractive index fluctuations. The
phases of light rays from space are disturbed by the refractive index fluctuations

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

3in the atmosphere. In the remainder of this thesis, we use a term of atmosphericturbulence for the turbulence of the refractive index in the atmosphere.

1.2.2

Structure Function

In order to describe a non-stationary random variation of the atmospheric variables, such as velocity, temperature, and refractive index, structure function[30]Df for the non-stationary random variable f at two diﬀerent points is defined as
Df (ρ) = ⟨[v(r) − v(r + ρ)]2 ⟩,(1.1)

where r = (x, y, z) is a vector of a three-dimensional coordinate, ρ is a spatial
separation between two points, and ⟨⟩ represents an ensemble average. The structure function corresponds to the intensity of the fluctuations in f (r) over scales
equal to or less than ρ. The turbulences of the atmospheric variables are a nonstationary random variables, that is, their mean value varies over time. However,in fact, this change occurs rather smoothly and therefore the structure function
Df can be considered a stationary random function, if ρ is not too large. The
turbulences of the atmospheric variables are also considered to be isotropic andhomogeneous in small scales. Therefore, the structure functions of these variables
only depend on a separation distance ρ = |ρ|. The structure function Df can be
represented by the correlation function Bf as
Df (ρ) = 2[Bf (0) − Bf (ρ)],(1.2)

where Bf (ρ) = ⟨f (r)f (r + ρ)⟩ and Bf (0) is a mean-square value of the variable f .
Kolmogorov found that a two-thirds power law for the velocity structure function Dv (ρ) ∝ ρ2/3 . Obukhov [41] showed that Kolmogorov’s law can be applied
to the temperature structure function, i.e. DT (ρ) = CT2 ρ2/3 , where CT is the temperature structure constant. Similarly, the refractive index structure function is
given by
DN (ρ) = CN2 ρ2/3 ,(1.3)
where CN2 is the structure constant of the refractive index with a unit of m−2/3 .
This constant CN2 represents the strength of the atmospheric turbulence and varies
with altitude and time. The vertical profile of this constant, CN2 (h), is commonlyused to characterize the condition of the atmospheric turbulence. The CN2 profile
is measured at many places. From the profiling of the CN2 with balloons [37], it isChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

4

found that the CN2 profile has peaks at several altitudes and therefore the atmospheric turbulence can be considered as the multiple discrete layers. In good sites
for astronomical observations, it is known that more than 60% of the turbulence
is concentrated in the ground layer [16], which is the turbulence lower than 0.5km, and additional turbulence layers are often within 4–10 km.
Tatarski [47] showed that the three-dimensional power spectrum of refractive
index can be expressed asΦN (κ) = 0.033CN2 κ−11/3(1.4)

where κ is a spatial wave number κ = 2π/ρ. This expression is only valid within
the range of L0 > ρ > ℓ0 , which is called inertial range.

1.2.3

The Near-Field Approximation

In the remainder of this thesis, we use the near-field approximation for the lightpropagation through the atmosphere turbulence, in which the diﬀraction is neglected and light rays propagate straight through the atmospheric turbulence.
Under this approximation, the phase distortion of the wave-front propagating toan aperture can be considered as the sum of phase distortions along the optical
path of the light ray in the atmosphere. This approximation is valid for weakturbulence found at good sites for astronomical observations.
With a characteristic size of the turbulence ℓ and a wavelength λ, the firstorder diﬀracted angle is given as θ = λ/ℓ. Diﬀracted rays interfere when the travel
distance is comparable to L ≈ ℓ/θ. From the turbulence profile at the top of theMaunakea measured by Thirty Meter Telescope project [16], the characteristic sizeof the turbulence ℓ is 0.25 m at altitude of 16 km in the worst case, which assuming
75 %ile model and low observation elevation of 30 degrees. This corresponds toa travel distance L = 130 km, which is much larger than the maximum of theatmospheric turbulence altitude. Therefore, the interference eﬀect is negligible atgood sites such as the top of Maunakea and the eﬀect due to the atmosphericturbulence can be described based on the near-field approximation.

1.2.4

Eﬀect of Turbulence

Let us consider that there is nlayer turbulence layers. In the near field approximation, the shift of the phase due to the refractive index fluctuations of i-thChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

5

turbulence layer at altitude hi with a thickness of δhi is∫ hi +δhi
ϕi (x) = k(1.5)

dz n(x, z),hiwhere ϕi (x) is a phase aberration at a two-dimensional coordinate x on i-th turbulence layer, k is a wave number k = 2π/λ, and n(x, z) is a refractive index
minus 1.
The total aberration due to nlayer layers is given simply by∑nlayerϕ(x) =(1.6)

ϕi (x).

i=1

The phase structure function for the Kolmogorov turbulence is represented bysubstituting Eq. (1.5) into Eq. (1.1) as f = ϕ and using Eq. (1.4) as(Dϕ (ξ) = 6.88

ξ
r0

)5/3
(1.7)

,

where ξ is a separation distance on a two-dimensional plane and r0 is Fried parameter [17] in unit of meter as
[

]−5/3

∫
2

r0 = 0.423k sec ζ

dhCN2 (h)

(1.8)

,

where ζ is a zenith angle. Fried parameter r0 is also referred to as a turbulencecoherence length and is the diameter over which a mean-square value of the phase
distortion corresponds to 1 rad2 at a wavelength of λ. Fried parameter r0 isconvenient measure of the strength of turbulence. The typical size of r0 is 0.15m at the visible length of 0.5 µm in good sites at night. Fried parameter varies
with wavelength as λ1.2 . This means that shorter wavelength is more aﬀectedby the turbulence compared with longer wavelength at the same CN2 . Similarly,
a coherence angle, which represents an angular separation within which a meansquare value of the phase diﬀerence between two wave-fronts propagating diﬀerent
path in the atmosphere corresponds to 1 rad2 at a wavelength of λ, is defined as
[

]

∫
2

8/3

θ0 = 2.914k (sec ζ)

dhCN2 (h)h5/3

.

(1.9)

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

6

This coherence angle is referred to isoplanatic angle. The isoplanatic angle isrepresented by r0 as
r0
θ0 = 0.314(cos ζ)(1.10)h̄
where h̄ is a weighted mean turbulence height,∫
dhCN2 (h)h5/3 .

h=(1.11)The power spectrum for the phase fluctuation is given by−5/3 −11/3Φϕ (κ) = 0.023r0

1.2.5

κ

.

(1.12)

Limitation of the Kolmogorov Model

As stated in 1.2.2, the Kolmogorov turbulence model is valid only within theinertial range, that is, for L0 > ξ > ℓ0 . In the astronomical imaging, the innerscale has only negligible eﬀect because the turbulence power is weak at the innerscale, but the outer scale has a major impact when the size of a telescope aperture
D is comparable to the outer scale. The turbulence with the spatial scale larger
than the aperture diameter of the telescope contributes the tip-tilt modes of thephase distortion over the aperture. This tip-tilt modes causes image motion on
the focal plane. In the case of long exposure, The size of point spread function(PSF), which is equal to the image of a point-like source, depends on this motion.
If the atmospheric turbulence follows the von Karman model and the outer scaleis finite, the PSF size aﬀected by the atmospheric turbulence become smaller than
the size expected by the Kolmogorov model.
The modified model including the outer scale was proposed by von Karman.
Using the von Karman model, the power spectrum Eq. (1.12) and the structure
function Eq. (1.7) for phase fluctuations are given as [9][−5/3Φϕ (κ) = 0.023r0(
Dϕ (ξ) = 0.172

r0
L0

κ2 +

(

2π
L0

)2 ]−11/6,()−5/3 [)]( )5/62πξ
2π 5/6ξK5/6
1−
Γ(5/6) L0
L0(1.13)(1.14)

where Γ is the gamma function and K denotes the modified Bessel function of thesecond kind. Fig. 1.1 shows the power spectrum and the structure functions of thephase fluctuations with the Kolmogorov model and the von Karman model. In

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

7

Spatial scale, ξ [m]
108

100

10

1

0.1

106

Kolmogorov
106

von Karman, L0 =30m

104
102
100
10−2

102
101
100

10−6

10−1

0.1

1

10

von Karman, L0 =30 m

103

10−4

10−8
0.01

von Karman, L0 =100 m

104
Dϕ (ξ) [rad2 ]Φϕ (κ) [m2 ]Kolmogorov
105

von Karman, L0 =100m

100

10−2
0.01

0.1

1

10

100

Spatial scale, ξ [m]Spatial wave number, κ [1/m]Figure. 1.1: The power spectrum (the left panel) and the structure function (the rightpanel) for phase fluctuations with the Kolmogorov model (the red solid lines)
and the von Karman model as L0 =100 m (the blue dashed lines) and L0 =30m (the green dotted lines). Fried parameter is set to r0 =0.2 m.

both figures, the von Karman model shows a flatter spectra in spatial scales larger
than the size of the outer scale.
The full width at half maximum (FWHM) is a diameter at which the intensityof PSF is a half of its peak value. With assuming the Kolmogorov model, FWHM
of long-exposure PSF, θKolmogorov , at a wavelength of λ is computed numerically bya fourier transform of the optical transfer function exp[−0.5Dϕ (λf )] [48], where fis the spatial frequency,
λ(1.15)
θKolmogorov = 0.98 .
r0Tokobinin [48] simulated the eﬀect of the outer scale on the seeing size numerically.
The FWHM of the PSF for the von Karman turbulence is approximated as√
θvonKarman ≈1 − 2.183

(

r0
L0

)0.356
θKolmogorov .

(1.16)

This approximation is obtained by fitting to the numerical results and shows good
agreement with numerical results for L0 /r0 > 22. Fig. 1.2 shows the seeing size
diﬀerence between the Kolmogorov and the von Karman model. The seeing sizeof the von Karman model becomes much smaller because there is no turbulence
which spatial scale is larger than L0 in the von Karman model. This feature ismeasured at 6.5 m Magellan telescopes [13].
A dome seeing, which occurs in a telescope dome, do not follow the Kolmogorov
turbulence. The dome seeing is diﬃcult to be modeled because it depends on thedome and telescope structure, temperature, and observation site. This dome seeingChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

8

0.55
Kolmogorov
von Karman, L0 =100m

FWHM of seeing [arcseconds]

0.5

von Karman, L0 =30m
0.45
0.4
0.35
0.3
0.25
0.2
0.15
0.5

1

1.5

2

2.5

wavelength λ [µm]Figure. 1.2: The seeing size for the Kolmogorov model (the red solid lines) from Eq.
(1.15) and the von Karman model as L0 =100 m (the blue dashed lines) andL0 =30 m (the green dotted lines) computed by Eq. (1.16). Fried parameteris set to r0 =0.2 m at λ=0.5 µm.

turbulence increases the turbulence power at ground and sometimes has a large
eﬀect on astronomical observations.

1.2.6

Taylor Hypothesis

The time evolution of the atmospheric turbulence can be expressed by two eﬀects :the boiling due to the energy cascade and the shifting due to wind. The first eﬀect
can be considered to be neglected within a short time scale and this assumption iswell known as Taylor’s hypothesis of frozen flow. In this condition, the atmosphericturbulence can be modeled simply by multiple turbulence layers moving due towinds with time. Furthermore, the frozen flow assumption allow us to predict thetime evolution of the atmospheric turbulences, i.e., the spatial displacement of aturbulence layer with wind v during time duration ∆t is v∆t. Many studies verify
the frozen flow assumption on-sky and suggest that the time scale in which thefrozen flow assumption is valid is range from a few ms to 100 ms, which depends
on conditions [21, 42, 46].

1.3

Principle of SCAO system

1.3.1

Single-Conjugate Adaptive Optics

Fig. 1.3 shows a schematic view of a SCAO system. In principle, an AO system
consists of a guide star (GS), wave-front sensor (WFS), controller, and deformableChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

Science
TargetAtmospheric Turbulence

Natural
Guide Star

9

Telescope
Collimator
Deformable
MirrorBeam
Splitter
Controller

Science
Instrument
Wave-front
SensorFigure. 1.3: The schematic view of a SCAO system.

mirror (DM). The process of AO correction is basically performed by the following
steps. First, the light from a GS is observed by a WFS to measure the wave-front
distortion in the direction of the GS. Second, a shape of a deformable mirror iscomputed from the measurement in order to compensate the wave-front distortion
when the DM reflects the light. Third, DM shape is controlled and the correction isapplied to the wave-front of the light from a science target. In most of the current
AO system, the wave-front distortion is measured in the optical wavelength andthe correction is done in NIR wavelength. This is because the AO correction works
better in longer wavelength than shorter wavelength, which is more sensitive tothe residual of the correction of the atmospheric turbulence. The optical lightfrom the GS is split by a beam splitter and fed to WFS.
Guide star
GS is required as wave-front reference. Although the best case is that sciencetarget itself is a bright point source and can be used as a GS, the sciencetarget is usually a faint and extended object, such as an extra-galaxy, andtherefore we have to find a bright star around the science target, which canbe used as a GS, it is called natural guide star (NGS). In this case, theChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

10

AO correction is performed based on the wave-front distortion in the NGS
direction. If the science target is far from the NGS, the light of the sciencetarget propagates through diﬀerent path from the path of the NGS at highaltitude in the atmosphere, as shown in Fig. 1.3. This diﬀerence makes theperformance of the SCAO correction worse. Therefore, the performance of
the SCAO correction decreases with the angular separation from NGS. The
permissible separation from the NGS is characterized by the isoplanatic angle θ0 , which is roughly 10 arcseconds in H -band even at good observingsites. The limitation is known as angular anisoplanatism. Due to the angular anisoplanatism, the shape of corrected PSF varies increasingly with theangular separation from the NGS. In addition, suﬃciently bright stars are
not available in all parts of the sky, and therefore, the sky-coverage of NGS
based AO system is very limited. Related to the sky-coverage, however, we
can improve by using artificial guide star, that is, laser guide star (LGS).
Wave-front sensor
The most used WFS is Shack-Hartmann WFS (SH-WFS). The SH-WFS
divides the image of the telescope aperture into the sub-apertures by amicro-lens-array and measures a wavefront gradients in each subaperture.
The wave-front propagating into SH-WFS can be reconstructed from themeasurement of the local gradients. This process is called wave-front reconstruction.
Deformable mirror
The correction is done by a DM. In general, tip-tile modes, which are thestrongest mode, are corrected by another tip-tilt mirror to relax the requirements for the stroke of the DM. The size of the element of a DM determines
the smallest spatial size which we can correct by an AO system. The turbulence with the spatial scale below the element size of the DM can not becorrected by the system, known as fitting error.
Control
For the real-time compensation, the AO loop needs to be controlled fasterthan the coherence time of the atmospheric turbulence [24] given byτ0 = 0.134/fG ,(1.17)

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

11

where fG is the Greenwood frequency given by
[
2

∫ ∞

fG = 0.102k sec(ζ)]3/5
CN2 (h)v(h)5/3 dh.

(1.18)

0where v(h) is a wind speed at altitude h. The typical value of the coherencetime τ0 is a few ms even at good observation sites. It means that the required
frequency for the AO control is order of hundreds Hz. Therefore, the required
computation time for the wave-front reconstruction is around 1 ms.
In order to optimize and stabilize the performance of AO system, we have toconsider the temporal control of AO system. The SCAO system is performedin closed loop, in which the WFS is put behind of a DM and measures
the residual wave-front distortion corrected by the DM, that is, there is
feedback from WFS into DM. The simplest controller is an integrator, butmore optimized controls are proposed. This temporal control is not in thescope of this thesis, but it is very important for an AO system.

1.3.2

Laser Guide Star AO System

LGS can improve the sky-coverage of an AO system. It is an artificial star createdby a laser beam. There are two kinds of the LGS system : Rayleigh LGS andsodium LGS. The first one uses Rayleigh scatter in the atmosphere and is createdat altitude of 10–20 km. The other one is created in the sodium layer at altitude
of around 90 km by exciting the state of sodium with a laser at wavelength of 589
nm. In this thesis, we focus on the sodium LGS, which is mainly used in the AO
systems of the current large telescopes.
Although using LGS allows us to create GS anywhere we want, there are somediﬃculties and limitation of LGS.
Insensitivity for the overall tip-tilt modes
Even if we use a LGS, another NGS is required to measure overall tip-tilt
modes, which is called tip-tilt NGS (TT-NGS). Since the light of LGS propagates through the atmosphere twice, going up from a laser launch telescopeand going down from the LGS, the eﬀect of the overall tip-tilt due to the atmospheric turbulence is cancelled out and therefore we can not measure theoverall tip-tilt modes by LGS. The TT-NGS can be fainter than the NGS.
In the case of NGS, the limit-magnitude is determined by the size of thesubapertures of SH-WFS. On the other hand, since the overall tip-tilt canChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

12

Science
TargetLGS

LGS

SodiumLayer

Laser Launch
Telescope

Aperture
Plane
SH-WFS
CameraTelescope

Cone Effect

Spot Elongation

Figure. 1.4: The schematic images for the cone eﬀect and the spot elongation due to the
LGS.

be measured from the image motion in only one aperture, the magnitude
limit of TT-NGS is determined by the size of the telescope aperture. In
addition, since overall tip-tilt modes have relatively large isoplanatic angle
compared with high order modes, the angular anisoplanatism for the overall
tip-tilt modes, known as tilt anisoplanatism, has weaker dependence on theseparation than angular anisoplanatism. Therefore, the permitted separation of TT-NGS from the LGS is larger than the one for NGS. Thus, using
the LGS can improve the sky-coverage compared with NGS AO system eventhough we need another TT-NGS.
Cone eﬀect
The LGS is created at finite altitude and its light propagates in the atmosphere with a conical optical path as shown in the left panel of Fig. 1.4. As
a result, the volume of the atmospheric turbulence measured by the LGS is
diﬀerent especially at high altitudes from the turbulence volume in which
the light of the science target propagates. This result in the deteriorationof the SCAO performance. This is known as cone eﬀect or focal anisoplanatism. The eﬀect increases with the size of the telescope aperture and will
be severe in the future ELTs.
Spot elongation of sodium laser on SH-WFS
From on-sky lidar measurements of sodium layer [5], it is observed that thesodium layer has a vertical thickness of roughly 10 km. Due to this thickness,

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

13

the angular size of the LGS spot is elongated in radial direction on the SHWFS as shown in the right panel of Fig. 1.4. The degree of the elongation,
θel , increases with the horizontal distance on the aperture plane between thesubaperture and the laser launch telescope, rsa , and is approximated asθel ≈rsa σNa
h2Na(1.19)

where hNa is the altitude of the sodium layer and σNa is the thickness of thesodium layer. If it is assumed that hNa =90 km, σNa =10 km, zenith direction,and the center launch system, in which the laser is launched from the back of
the secondary mirror, the elongation at the outermost subaperture is roughly
1 arcseconds for a 8 m telescope and 4 arcseconds for a 30 m telescope. Thiselongation causes spot position estimation errors for outer subapertures [19].
In reality, the vertical profile of the sodium layer has complex structure andvaries with time, and the spot elongation pattern also varies.

1.3.3

Performance Metrics

There are diﬀerent metrics to evaluate the performance of AO correction. One is awave-front error (WFE) which is computed from WFS telemetry data. Others are
FWHM, Strehl ratio (SR), and ensquared energy measured from observed PSF.
Wave-front error
WFE is defined as the root mean square (rms) value of the residual wavefront distortion. If the residual phase distortion of the wave-front is denotedby φ, the WFE σ is expressed as
√
σ=

1
S

∫
(φ(x) − φ̄)2 dx(1.20)

where S is the telescope pupil area and φ̄ is an average of the residual
distortion over the telescope pupil.
Full width at half maximum
FWHM is a diameter at which the intensity of PSF is a half of its peakvalue. The FWHM aﬀected by the atmospheric turbulence is represented
by Eq. (1.15) and Eq. (1.16). It is noted that the FWHM is often measured
by fitting two-dimensional gaussian or moﬀat function to the observed PSF
for convenience, but, to be exact, the gaussian and moﬀat function are notChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

14

good approximation of the PSF corrected by an AO system because the AO
corrected PSF has two components, which are a corrected bright core anduncorrected diﬀuse halo.
If there is no turbulence, the PSF is determined by the diﬀraction due to optics, which is called diﬀraction-limited, and have diﬀraction pattern, known
as airy disk for circular aperture. The radius of the innermost dark ring of
the diﬀraction-limited PSF is given byθdif = 1.22

λD

(1.21)

where D is aperture diameter. The FWHM of the diﬀraction-limited PSF isgiven by 1.03λ/D. The FWHM of the uncorrelated PSF is represented by
Eq. (1.15) or Eq. (1.16) depending the turbulence model.
Strehl ratio
Strehl ratio (SR) is the ratio of peak intensity of observed PSF to the peakintensity of the diﬀraction-limited PSF, has a value from 0 to 1. The goal of
AO system is maximizing the SR by the correction. The SR can be estimatedfrom the WFE by using Mareshal approximation [24],SR ≈ e−σ .
2(1.22)

Although this expression is very useful, it is accurate to a couple of percent
for rms errors of 1/10λ.
Ensquared energy (EE)The ensquared energy (EE) is a fraction of intensity within a defined square,
which is used to evaluate AO performance for spectroscopic observation by
setting a square size to the width of a slit.

1.4

Wide-Field Adaptive Optics

While SCAO systems with both NGS and LGS are well-established and provides
near diﬀraction-limited spatial resolution in NIR wavelength on the current large
telescopes, there are scientific demands to apply AO correction to wide FoR to
increase an observation eﬃciency or/and PSF uniformity over the wide FoR. As
mentioned in the previous section, the size of FoR of a SCAO system is limitedby the angular anisoplanatism. In order to overcome the angular anisoplanatism,

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

FoR

FoR

15

Each FoVLGSLGS

Conjugated
plane

WF
S

High-altitude
DM

WFS

WFS

WF S

Low-altitudeDM

optical conjugation

optical conjugation

Conjugated
plane

S

WF

WF

S

DM

Science
Instrument

Science
Instrument(a) MCAO(b) MOAO
FoR

LGS

S

WF

WF

S

WFS

optical conjugation

Ground-Layer
DM

Science
Instrument(c) GLAO

Figure. 1.5: The schematic views of WFAO systems.

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

16

three WFAO concepts with using multiple LGSs are proposed: multi-conjugate
adaptive optics (MCAO), multi-object adaptive optics (MOAO), and ground-layer
adaptive optics (GLAO). The schematic views of these systems are shown in Fig.
1.5. Each WFAO system is design for diﬀerent purposes. In this section, we
introduce principle, advantage, disadvantage, and current/future systems of each
WFAO concepts briefly.

1.4.1

Multi-Conjugate Adaptive Optics

MCAO system was proposed by Beckers [6]. MCAO system requires the tomographic wave-front reconstruction to reconstruct a three-dimensional structure of
the phase distortion due to the atmospheric turbulence with using multiple LGSs.
Then, MCAO system compensates for the phase distortion three-dimensionally bymultiple DMs based on the reconstructed phase distortion. The number of DMs
is usually 2–3, which are put in series and conjugated at diﬀerent heights in theatmosphere (the blue planes in the Fig. 1.5(a)). MCAO system provides a uniform
correction over a wide field of view (FoV) of 1–2 arcminutes. This is a main advantage of MCAO system and good for the precise photometry and astrometry. In
addition, using multiple LGSs can reduce the cone eﬀect of LGSs, which however
is not limited to a MCAO system. This is important for future ELTs, where thecone eﬀect will be severe. As shown in Fig. 1.5(a), the measurement of wave-fronts
is done in closed-loop in MCAO system, in which WFSs are put after DMs andmeasure the residual wave-front distortions.
MCAO system suﬀers from generalized fitting error due to the discrete numberof DMs [43]. Consider the case of three turbulence layers and two DMs conjugatedto diﬀerent altitude from the turbulence layers, as Fig. 1.5(a). Even if the tomography perfectly estimates the turbulence layers, it is impossible to provide theperfect correction for all direction by the two DMs. This is known as the generalized fitting error. Increasing the number of DMs is a solution of this problem,
but it requires additional optical surfaces, which decreases throughput, and complexity of the control. Another diﬃculty associated with a MCAO system is therequirement of large DMs to compensate the phase distortions at high altitudes.
These diﬃculties limit the size of FoV size of a MCAO system to the diameter of
1–2 arcminutes, which is smaller than other WFAO concepts.
At this moment, one MCAO system works as a facility instrument for GeminiSouth 8 meters telescope, named the Gemini Multi-conjugated adaptive optics
System (GeMS) [40, 44]. The GeMS computes the tomography by using 5 Na

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

17

LGSs and 1–3 TT-NGSs, which measure unseen modes of the LGS, and corrects
the phase distortion over the wide FoV more than 1 arcminutes by three DMs
conjugated to 0, 4.5, 9 km. Another MCAO system is developed for the Large
Binocular Telescope (LBT) [26]. MCAO systems are also designed as first light
instruments for future ELTs, for example, NFIRAOS [27] for Thirty Meter Telescope (TMT) and MAORY [14] for European ELT (E-ELT). The NFIRAOS on
TMT will do the tomographic wave-front reconstruction using six LGSs and provide the best correction for 30 arcseconds diameter FoV and the mild correction
for 2 arcminutes FoV by three DMs.

1.4.2

Multi-Object Adaptive Optics

MOAO system was proposed originally by Hammer et al. [23] and also requires
the tomographic wave-front reconstruction. MOAO system provides the optimized
correction for multiple small patches in the large FoR of 5–10 arcminutes diameter. As shown in Fig. 1.5(b), there are multiple science pick-oﬀ arms on sciencetargets. After each pick-oﬀ arm provides the light from the science targets foreach science channel, a DM contained in each science channel makes the optimalcorrection in its science direction. This parallel approach can increase an eﬃciencyof the observation. In addition, since each DM in a MOAO system corrects thephase distortion in one direction, there is no generalized fitting error in MOAO
systems, and the size of each DM can be comparable to the size of DM in SCAO
systems. Therefore, MOAO systems can realize larger FoR, which is considered
5–10 arcminutes diameter, compared with MCAO systems. The multiplicity of
MOAO systems is suited for multi-object observations of galaxies in high redshift
universe because these high redshift galaxies have a small apparent angular size
and a low number density.
However, MOAO is still considered as very challenging system. The mostchallenging point of a MOAO system is open-loop control. MOAO system corrects
the phase distortion in diﬀerent directions from WFSs and therefore, we need theopen-loop control. In the open-loop control, since there is no feedback fed to the
WFSs from the DMs, it is required excellent calibration and stability of all DMs
and WFSs compared with closed-loop AO systems. Moreover, open-loop WFSs
observe the raw phase distortions not corrected by DMs. It means that WFSs
with the large dynamical range are required. We have to construct the opticaland mechanical design of future MOAO systems with considering how we calibrate
and control the system. The complexity of MOAO system is the other challenging

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

18

point. Because the FoV of each science path is small and the optical design ofeach science path can be simple, there can be science channels as many as 10 infuture MOAO systems [3, 45]. Therefore, and the number of moving elements,
such as DMs and pick-oﬀ arms, is much larger than other WFAO concepts.
In order to identify technical readiness and diﬃculties associated with an
MOAO system, there are several projects to develop an MOAO demonstrator
and to test it on-sky. Here, we introduce two projects : CANARY [18] for William
Herschel Telescope (WHT) of 4.2 meters diameter and RAVEN [32] for Subaru
telescope of 8 meters diameter. The CANARY is the first on-sky demonstratorof an MOAO system, which works with three NGSs, four Rayleigh LGSs, andone on-axis science channel in the phase B, and is developed as a pathfinder for
EAGLE, which is MOAO system on E-ELT [45]. RAVEN is the first MOAO technical and science demonstrator on 8-m class telescope. RAVEN is developed andtested for the Subaru telescope on Maunakea in Hawaii. RAVEN provides MOAO
correction for two science targets simultaneously by using three NGSs and one Na
LGS installed in the Subaru telescope. We will explain the system of the RAVEN
later.

1.4.3

Ground Layer Adaptive Optics

GLAO system can be considered as a variation of the MCAO system with only
one DM conjugated to the ground, as shown in Fig. 1.5(c). As mentioned above,
more than half of the atmospheric turbulence are caused in the ground layer.
GLAO provides the correction only for the ground layer and improves the seeing
rather than achieves near diﬀraction-limited resolution. The ground layer is acommon component to wide range of directions. Therefore, GLAO system provides
moderate correction over the wide FoV, which is possibly more than 10 arcminutes
diameter. This is useful to survey observations over large field. Mostly, an adaptive
secondary mirror is considered as a DM for a GLAO system, then we can perform
GLAO correction without any additional optical component. In a GLAO system,a tomographic reconstruction is not necessary since the ground layer turbulence
can be computed by averaging measurements from multiple WFSs.
A GLAO system is working in LBT [22], named ARGOS. This is a GLAO
system with three Rayleigh LGSs and an adaptive secondary mirror. ARGOS
provides 0.3 to 0.4 arcsecond images in J, H and K with GLAO correction from
0.7 to 1 arcsecond uncorrected seeing. Furthermore, other large telescopes plan a
GLAO system with adaptive secondary mirror.

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

19

Value
Telescope aperture

30 m

Total system throughput

30 %

Dark

0.006 e− /pixel/s

Readout

5 e− /pixel/readIntegration Time10, 20, and 30 hours

FoR

5, and 10 arcminutes diameter

AO performance

50 % EE

Aperture for EE

50 mas box for point source
300 mas box for extended source

Signal to noise ratio (SNR)
for detection

5

Spectral Resolution

3000

Observed wavelength

J-band (1250 nm)

Tab. 1.1: Parameters used in the estimation of the number of detected high redshift
LBGs in 1.4.4. The dark count and readout noise is referenced from the TMTIRIS simulation [51].

1.4.4Importance of Wide Field

The key parameters of WFAO systems are the performance of AO correction andthe FoR size. We show the importance of these parameters with respect to multiobject spectroscopy of lyman-break galaxies (LBGs) in high redshift universe,z > 5.
From deep imaging surveys by the Hubble Space Telescope (HST), many candidates of galaxies in high redshift universe z > 5 are detected by the lyman-breakmethod [7]. It is known that these high redshift LBGs are star-forming galaxies. In order to understand the star formation in the early universe, multi-object
spectroscopy of rest-frame UV continuum, emission lines, and absorption lines forlarge number of high redshift LBGs is required. Such a observation is consideredas the main science of a MOAO instrument in the future ELTs [1]. As the low
number density of high redshift LBGs, the large FoR is required to maximize themultiplicity of the MOAO for high redshift LBGs.
Here, we estimate the number of high redshift LBGs expected to be detected bya multi-object spectroscopy with a MOAO system on the 30 m aperture telescopeand compare the achieved number density with the observed luminosity function.
The parameters used in this estimation are summarized in Table 1.1. The total

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction
M1600 [mag]M1600 [mag]

10−3

−23.5

Number Density [Mpc−3 ]Number Density [Mpc−3 ]−23 −22.5 −22 −21.5 −21 −20.5 −20 −19.5

FoR=5′ φ10−420

FoR=10′ φ
10−5

10−6

−23 −22.5 −22 −21.5 −21 −20.5

−20

10−3

FoR=5′ φ10−4

FoR=10′ φ
10−5

10−6
25

25.5

26

26.5

27
27.5
mJ [mag]

28

28.5

29

25

25.5

26

26.5

27
27.5
mJ [mag]

28

28.5

29

Figure. 1.6: Observed number density and the number density achieved by the multiobject spectroscopy with the future MOAO system at z ∼ 5 (the left panel)and z ∼ 6 (the right panel). The gray curve indicates the observed cumulative luminosity function of high redshift LBGs as a function of an apparent magnitude at J-band. The observed luminosity functions are modeledwith a Schechter function [7], and the parameters of the Schechter function
are MU∗ V = −21.17, ϕ∗ = 0.74 × 10−3 , and α = −1.76 for z ∼ 5, and
MU∗ V = −20.94, ϕ∗ = 0.50 × 10−3 , and α = −1.87 for z ∼ 6. The redshift
range for each redshift is assumed to 4.5 < z < 5.5 for z ∼ 5 LBGs and5.5 < z < 6.5 for z ∼ 6. The vertical lines indicate the limit magnitudeachieved by the future instrument for the point source (solid) and extended
source (dotted) as an integration time of 10 hours (red), 20 hours(blue),and 30 hours (green). The horizontal lines represent the cumulative number
density required to detect 10 (solid) and 5 (dotted) LBGs over the FoR of
5 (magenta) and 10 (black) arcminute diameter.

system throughput, including the telescope, the AO system, and the science instrument, is assumed to be 30 %. We consider a point source and an extended
source as a intensity profile of LBGs. For the point source, the performance ofthe MOAO system is assumed to be 50 % EE in a 50 mas box over the FoR. For
the extended source, the MOAO performance is 50 % EE in a 300 mas box tomaximize the signal to noise ratio (SNR). The background and transmittance ofthe sky are referenced from the measurement by GEMINI observatory1 . The darkcount and readout noise of the detector are referenced from [51]. The observation is performed at J-band, which is close to the redshifted wavelength of therest-frame UV continuum. The criteria for the detection is set to SNR=5.
Fig. 1.6 shows a relation between the observed cumulative luminosity function
and the number of LBGs detected by the multi-object spectroscopy with the future
MOAO system at z ∼ 5 (the left panel) and z ∼ 6 (the right panel). The gray curve
indicates the observed cumulative luminosity function of high redshift LBGs. The
1

http://www.gemini.edu/sciops/telescopes-and-sites/observing-condition-constraints/irbackground-spectraChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

21

vertical lines indicate the limit magnitude achieved by the future instrument forthe point source (solid) and extended source (dotted) as an integration time of 10hours (red), 20 hours(blue), and 30 hours (green). The horizontal lines represent
the cumulative number density required to detect 10 (solid) and 5 (dotted) LBGs
over the FoR of 5 (magenta) and 10 (black) arcminute diameter. For example, atz ∼ 5, an intersection point of the red vertical solid line and the magenta horizontal
solid line is above the observed cumulative luminosity function of LBG, and itmeans that the future instrument can not detect the rest-frame UV continuum of
10 LBGs at z ∼ 5 within the FoR of 5 arcminutes diameter with a SNR of 5 andan integration time of 10 hours, if the LBGs are point sources. On the contrast,
since an intersection point of the red vertical solid line and the black horizontal
solid line is below the observed cumulative luminosity function, the rest-frame UV
continuum of 10 LBGs at z ∼ 5 can be detected by the future instrument within
the FoR of 10 arcminutes diameter with a SNR of 5 and an integration time of 10hours. Similarly, if the LBG is extended object, while 5 LBGs can not be detectedwith a 5 arcminutes FoR and the integration time of 30 hours, it will be possible
with a 10 arcminutes FoR. This trend is also seen at z ∼ 6. These estimates showthat in order to maximize the multiplicity for the high redshift LBGs, expanding
the FoR without reducing the AO performance is important.

1.5

Tomographic Wave-Front Reconstruction

In this section, we provide more introduction about the tomographic wave-front
reconstruction, which is main interest of this thesis. As mentioned in subsection
1.4.4, the key parameters of WFAO systems are the performance of AO correction
and the FoR size. The AO performance and the size of FoR depend strongly on
the accuracy of the tomographic wave-front reconstruction.
The tomographic wave-front reconstruction has already been demonstrated onsky by some instruments and actually it works on sky[32, 40]. However, the estimation error due to the tomographic wave-front reconstruction causes a significant
WFE, called tomographic error, and the tomographic error is a one of major error term in WFAO systems. Therefore, improving the tomographic wave-front
reconstruction directly results in the enhancement of the WFAO performance.
In this section, we will explain the principle of the tomographic wave-front
reconstruction briefly. Then, we will explain the problem associated with the tomographic wave-front reconstruction, which limits the the performance and theChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

22

FoR size of WFAO system. Finally, we will explain the computational complexity, which is the other challenging point especially for future ELT-scale WFAO
systems.

1.5.1

Tomography Step and Fitting Step

The goal of the tomographic wave-front reconstruction is reconstructing a threedimensional structure of the phase distortion due to the atmospheric turbulence
and computing the shapes (or commands) of DMs from the reconstructed phase
distortion. The first part of the tomographic reconstruction is often referred to asa tomography step and the second part is a fitting step.
The measurements from multiples WFSs contain the information of the phase
distortion in corresponding GS directions. In order to overcome the angular anisoplanatism, we have to consider the optical path diﬀerence at high altitude depending on the direction. This is achieved by converting the measurements from WFSs,
which probes the phase distortion in each direction, into the phase distortion ateach altitudes using the tomographic technique. In the case of SH-WFSs, thisprocess is equal to solving an inverse problem for a linear equation between themeasurements from SH-WFSs and the phase distortion due to the atmosphericturbulence as
s = Hϕ + η(1.23)
where s and η are a vector of the measurements and measurement noises of multiple SH-WFSs, ϕ is a vector of phase values on a discrete grid of points on theatmospheric turbulence layers, and H is an influence matrix between measurement s and the phase distortion ϕ. The influence matrix H is created basedon the GSs direction, the altitudes and strengths of the atmospheric turbulencelayers, and SH-WFS parameters. The altitudes and strengths of the atmosphericturbulence layers vary with time. Therefore, it is required to measure these information during observation and update the influence matrix H to optimize theperformance of the tomographic wave-front reconstruction.
The fitting step depends on with the WFAO concepts. In MCAO system,
the correction is performed by the series of DMs, which conjugated at diﬀerent
altitudes. The commands of DMs are determined to provide the uniform correction
over the FoV. On the other hand, in MOAO system, the phase at each DM-plane
is optimized for each direction and the fitting step is performed for each DMs inparallel.

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

23

In a simple case, these steps are performed with one matrix-vector multiplication as u(t) = Es(t), where u is a vector of the commands of DMs and E isa matrix performing the tomography step and the fitting step. The matrix E iscalled reconstructor or reconstruction matrix. With following the update of theinfluence matrix H, the reconstruction matrix E needs to be updated.

1.5.2

Lack of Information in Tomographic Reconstruction

The performance of the tomographic wave-front reconstruction depends strongly
on the number of GSs. However, the number of GS is limited in astronomical
AO systems. From the limited number of GSs, we can not get information of thephase distortion in a volume of the atmosphere corresponding to the FoR fromthe WFS measurements. This lack of information limits the performance of thetomographic wave-front reconstruction and also makes it diﬃcult to expand theFoR size of WFAO systems.
Fig. 1.7 shows the schematic view of this problem. There are three turbulencelayers and three GSs. The footprints of the optical path of GSs on the turbulencelayers are indicated with orange circles. In the figure, the some parts of turbulencelayers are not covered by the GS optical paths. Namely, we can not get informationof the phase distortion of this uncovered area from WFSs. The uncovered area
causes a significant tomographic error mostly on the outer area of the FoR. Thisis one of the lack of information in the tomography.
The other lack of information is caused by the area referred to as unoverlapped
area. In Fig. 1.7, most of area on the atmospheric turbulence at high altitude iscovered by only one GS optical path. The WFSs probes the integrated phase
distortion in its GS direction. In other words, the phase distortions from multiple
turbulence layers are degenerated in the WFS measurement. In order to solve thisdegeneracy and reconstruct the phase distortion on each atmospheric turbulence
screens by the tomographic wave-front reconstruction, the information from multiple directions is necessary, that is, all areas should be covered by two or more GS
optical paths. This degeneracy due to the unoverlapped area results in a significant
tomographic error even in the direction within the GS asterism.
In particular, the tomographic error due to these areas becomes more severe
when one considers to expand a FoR of WFAO systems with the wide configurationof GSs. In addition, using LGS increases the geometric error compared with using
NGS due to its conical optical path. This means that the tomographic error willChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

24

Field of Regard	
Guide Star	

Atmospheric	
Turbulence	Meta-pupil	

Uncovered area	
Unoverlapped area	
GS footprint at a 	previous time-step	

Telescope	

GS footprint at a 	current time-step	

Telescope	

Figure. 1.7: Geometric relation between positions of GSs and regions in the atmospheremeasured by WFSs. Left: The solid orange circles show the footprint of
GS optical paths at the current time-step. There are areas covered by noor only one GS, referred to as uncovered area or unoverlapped area. These
areas cause significant tomographic errors. Right: The blue arrows indicate
wind direction at each altitude. The GS footprints measured at a previous
time-step move due to the wind with time, which are dashed orange circles
in the figure.

be more severe in WFAO systems for the future ELT with multiple LGSs because
the cone eﬀect gets larger with increasing the size of aperture.
Some approaches are proposed to overcome this issue[2, 11]. These methods
are based on the frozen flow hypothesis. Under the frozen flow assumption, thetime evolution of the atmospheric turbulence can be computed as the shift due
to the wind. In this case, as shown in the right panel of Fig. 1.7, the areas in theatmospheric turbulence measured at previous time-steps also moves with time.
This means that the estimation and the measurement at previous time-steps canbe considered as the information of displaced areas at the current time-step. Ammons et al. [2] propose the multi-layer shift-and-average technique. They averaged
the reconstructed phase distortions over multiple time-steps for each turbulencelayer with considering the shift due to the wind. Their method can reduce thetomographic error due to the uncovered and unoverlapped areas successfully in the
numerical simulation. Moreover, predictive linear quadratic gaussian (LQG) control probably can reduce the tomographic error due to the lack of information[11].
The predictive LQG is proposed mainly for reducing the time lag error. In the
LQG control, the estimations at multiple time-steps are integrated. Therefore,
the LQG may works not only for the time lag error, but also for the tomographic
error due to the lack of information.

Chapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

25

However, It is noted that in the methods above, an estimate at each time-stepis already aﬀected by uncovered and unoverlapped areas. In order to solve the lackof information directly, we should use information from multiple time-steps in the
tomographic reconstruction, not after the tomographic reconstruction such as themulti-layer shift-and-average technique and the LQG.
Estimating wind speeds and directions at diﬀerent altitudes is essential for
our tomography method. Several methods to estimate wind speeds and directions
from measurements of Shack Hartmann WFS (SH-WFS) are proposed and tested
by using on-sky measurements[21, 42]. These methods successfully detect the windat multiple altitudes and start to be applied to the real-time operation.

1.5.3

Computation Complexity

The computation complexity of the tomographic reconstruction is the other problem especially for future ELT-scale WFAO system. For the real-time compensationof the phase distortion, we have to compute the tomographic reconstruction within
1 ms as mentioned in 1.3.1. Furthermore, the tomographic reconstructor should
be updated every from a few ten seconds to a few minutes to follow the change ofatmospheric condition.
For example, consider MOAO systems on a current large telescope with 8 m
diameter and a future ELT with 30 m diameter. The number of GSs and science channels are assumed to be 5 and 10, respectively. The size of subaperturesof SH-WFSs and elements of DMs are 0.5 m, that is, there are 16times16 subapertures in a SH-WFS and 16×16 elements in a DM for the MOAO system on
the current MOAO system and there are 60×60 subapertures and elements forthe future MOAO system. The total number of measurements from 5 SH-WFSs
is 2560 and 36000, respectively, where it is noted that each subaperture measures local gradients along x and y direction. The total number of estimates of
the tomographic reconstruction, i.e. the number of elements of 10 DMs, is 2560
and 36000, respectively. The size of the tomographic reconstructor for the future
ELT-scale MOAO system is roughly 200 time larger than the size for the MOAO
system for the current large telescopes. The complexity for computing the tomographic reconstructor with simply inversing the influence matrix, H, scales as
O(n3 ), where n is a dimension of the reconstructor matrix. Also, the complexity ofthe matrix-vector multiplication for the tomographic reconstruction scales O(n2 ).
Thus, the complexity of the tomographic reconstruction increase largely with thesize of the aperture. In the ELT-scale WFAO systems, we can not achieve theChapter 1. Introduction to Wide-Field AO Systems and Tomographic
Reconstruction

26

required control speeds with the standard computation method. Therefore, more
eﬃcient methods are required.
There are a lot of eﬀort to reduce the computation time of the ELT-scale tomographic reconstruction. In the ELT-scale tomographic reconstruction, we can notcompute the tomographic reconstructor directly because it is too large. Therefore,
the iterative methods, such as a conjugate gradient method, are proposed, whichsolves the inverse problem for the tomographic reconstruction iteratively at every
time-step. Ellerbroek [15] proposed the sparse matrix technique, which uses thesparseness of the matrix and can reduce the computation time of the matrix-vector
multiplication. The most of matrices in the tomographic reconstruction are known
as a sparse matrix, and therefore, the sparse technique allows us to reduce the computation time for the iterative method. Also, the conjugate gradient method with
an eﬃcient preconditioning technique [49] and a warm start technique [36] is proposed to reduce the number of iterations. Recently, the distributed Kalman filter
is developed as a fast non-iterative spatially invariant controller [38]. In addition,
the general purpose computing on GPU (GPGPU) can be powerful computing
engines of those algorithms.

1.6

Object of This Thesis

In this work, we will develop a new algorithm for the tomographic wave-front
reconstruction, which reduces the tomographic error due to the lack of information and to extend the size of FoR of WFAO systems without reducing the AO
performance, and the method to estimate wind speeds and directions at multiple altitudes, which is essential for the new tomography algorithm. The basic
idea of the new tomographic reconstruction algorithm is to solve the tomographic
wave-front reconstruction by using measurement from multi time-steps based on
the frozen flow assumption. The main goal is to demonstrate and evaluate theperformance of the new tomographic reconstruction by analytical method, numerical simulation, laboratory test, and on-sky observation. Also, we will discuss the
GPGPU capability for the future-ELT scale tomographic reconstruction.

Chapter 2
Methodology
2.1

Tomographic Wave-Front Reconstruction

Fig. 2.1 shows a schematic of geometry for a MOAO system with multiple LGSs
and SH-WFSs. The goal of tomographic reconstruction is reconstructing DM
commands of each science direction from the measurements of multiple SH-WFSs.
We consider that the atmospheric turbulence consists of Nlayer thin layers locatedat diﬀerent altitudes. The tomography is computed from Nlgs LGSs for Nobj
science objects. The directions of each LGS and each science target are indexed by
α = (αx , αy ) and β = (βx , βy ), respectively. The command, or voltage, provided
to the elements of DM for kth science object is indicated by a column vector uβk .
In the near-field approximation, the phase distortion of the aperture-planewave-front φθ (x, t) in the science target direction βk at time t is represented as
the sum of the phase distortion caused by the turbulence layers,∑

Nlayerφβk (x, t) =ϕi (x + hi βk , t),

(2.1)

i=1

where x is two-dimensional spatial coordinate vector x = (x, y), hi is an altitude of
ith turbulence layer, and ϕ(x, t) is a phase distortion at a point x on ith turbulencelayer at time t. In the case of LGS in direction αj at altitude hNa , the scaling of
the optical footprint at high altitudes due to the cone eﬀect should be consideredas)
Nlayer∑ ( hNa − hi
ϕi x
+ hi αj , t .
(2.2)
φαj (x, t) =h
Na
i=1
In the later discussion, wave-fronts and turbulence layers are discretized spatially and denoted by column vectors, this is, φ(t) for wave-fronts and ϕ(t) for27

Chapter 2. Methodology

28

FoR
TargetLGS

Meta-pupil

Wind

Turbulent
Layer

Wind

Turbulent
Layer
Telescope
Aperture

LGS
Footprint
TargetFootprint

Wave-front
at pupil plane
SH-WFS
MeasuredSlope

SH-WFS
DM

Figure. 2.1: Image of a MOAO model with multiple turbulence screens, LGSs, and SHWFSs. β and α indicate directions of science targets and LGSs. The index
i, j, k is used for the turbulence screens, LGSs and SH-WFSs, and sciencetargets and DMs, respectively. The goal of the tomographic reconstruction
is to estimate DM command uβ from the measured slope sα .

turbulence layers. We define a ray-tracing matrix Pθi , which extracts phase distortion within a optical path footprint in the direction θ from i-th atmosphericturbulence by using a bilinear interpolation. Using the ray-tracing matrix, Eq.
(2.1) and Eq. (2.2) are rewritten with considering all phase points in the apertureas∑

Nlayerφβk (t) =Pβik ϕi (t) = Pβk ϕ(t),(2.3)

Pαi j ϕi (t) = Pαj ϕ(t),

(2.4)

i=1
Nlayerφαj (t) =∑i=1

where Pαi j is a ray-tracing matrix including the scaling eﬀect due to the coneeﬀect, ϕ(t) = [ϕT1 (t) · · · ϕTNlayer (t)]T , and Pβk and Pαk are a concatenation of all
submatrices Pβik and Pβik , respectively.

Chapter 2. Methodology

29

We define the relation between the DM commands reconstructed by the tomography and the wave-front made by the DM for kth science direction asφ̂βk (t) = Nβk ûβk (t)(2.5)

where the hat symbol indicates estimated quantities and Nβk is a interaction
matrix for kth DM. The interaction matrix is determined experimentally.
The phase variance of the aperture-plane residual wave-front for science direction βk is defined as(2.6)
ϵβk = ⟨||φβk − φ̂βk ||2 ⟩.
In MOAO system, a tomographic reconstructor is determined to minimize theresidual phase variance ϵβk [10],Eβk = argmin ϵβk ,

(2.7)

Eβk

where Eβk is the tomographic reconstructor for the direction βk and ⟨⟩ indicatesensemble average over time. This reconstructor is known as minimum variancereconstructor [15], which provides the solution to minimize the variance of theresidual phase at on aperture-plane. The DM commands uβk is reconstructedfrom the measurements sα by the reconstructor Eβk as uβk = Eβk sα .

2.1.1

Measurement Model

The slope sαj (t) measured by the j-th SH-WFS aiming the direction αj at time
t is defined assαj (t) = Γαj φαk (t) + ηj (t) = Γαj Pαj ϕ(t) + ηj (t),(2.8)

where Γαj is a discrete phase-to-slope operator which converts phases into slopesand ηj (t) is a column vector of the noise in measurements of j-th SH-WFS.
We use Fried geometry to define Γαj . In the Fried geometry, slopes for x andy direction in a subaperture are defined by surrounding 4 phase points (Fig. 2.2)
as


φ
i,j
( )
(
)

x

si,j
−1 1 −1 1  φi+1,j 
1

(2.9)
=

2dsa −1 −1 1 1 syi,jφ i,j+1 
φi+1,j+1

Chapter 2. Methodology

30

Figure. 2.2: Fried geometry to define slopes for SH-WFS.
x(y)where φi,j and si,j are a phase value and x(y)-slope at a discrete point indexed iand j. Eq. (2.8) is concatenation of Eq. (2.9) for all subapertures.
Concatenating Eq. (2.8) of Nlgs SH-WFSs, we can obtain a equation connecting
the phase distortion on Nlayer atmospheric layers and the slope provided by Nlgs
SH-WFSs as follows.

sα1 (t)..
 = Γα Pα ϕ(t) + η(t)sα (t) = (2.10).

sαNlgs (t)where Γα is a block diagonal matrix, Γα = diag[Γα1 , · · · , ΓαNlgs ], Pα is a raytracing matrix concatenating Pαj for 1 ≤ j ≤ Nlgs , and η(t) is a column vector
including measurement noises from all SH-WFSs.

2.1.2

Minimum Variance Reconstructor

The minimum variance reconstructor can be obtained by minimising Eq. (2.6). As
shown in Appendix A.1 , the final formula is given as−1 −1 T T −1
Eβk =(NβTk Nβk )−1 NβTk Pβk (PαT ΓTα Σ−1ηη Γα Pα − Σϕϕ ) Pα Γα Σηη=Fβk Eϕ .

(2.11)(2.12)

where Σηη is a noise covariance matrix, Σϕϕ is a covariance matrix of the phase
distortion of the atmospheric turbulence, Eϕ is a estimation matrix, and Fβk is afitting matrix. If it is assumed that the measurement noise from jth WFS is zero2
and noises of all WFS subapertures aremean gaussian noise with a variance of ση,jindependent of each other, the noise covariance matrix becomes a block diagonalChapter 2. Methodology

3122
matrix as Σηη = diag[ση,1I, · · · , ση,NI]. Furthermore, under the assumption
lgs
that all WFSs have same noise with a variance ση2 , the noise covariance matrix
becomes Σηη = ση2 I and the estimation matrix Eϕ can be rewritten as−1 T TEϕ = (PαT ΓTα Γα Pα + ση2 Σ−1ϕϕ ) Pα Γα .

(2.13)

Then, the DM commands for science direction βk (t) is reconstructed from measurements as ûβk (t) = Eβk sα (t) = Fβk Eϕ sα (t).
The term of ση2 Σ−1
ϕϕ works as a regularization matrix, which smooths the tomographic errors of high spatial frequencies due to the measurement noise. The
eﬀect of the regularization matrix is scaled by the balance between the measurement noise ση2 and the turbulence power CN2 . The phase covariance Σϕϕ can becomputed based on the turbulence powers, altitudes, and the Kolmogorov or the
von Karman power spectra, and is a block diagonal matrix including the phasecovariance matrix at each altitude as Σϕϕ = diag[Σϕϕ,1 , · · · , Σϕϕ,Nlayer ].
This reconstructor is constructed based on the model of the turbulence structure and SH-WFSs, which is known as model-based reconstructor. The turbulence
altitude, CN2 values, measurement noise are required to create the model-based
reconstructor. We will explain the method to estimate the CN2 profile in . In addition, the DM commands is estimated from the measured slopes at one time-stepby this reconstructor. Therefore, Eβk is referred to as single time-step tomographic
reconstruction in the remainder of this thesis.

2.1.3

Discrete Laplacian Approximation

As mentioned, the oﬀ-line and on-line computation complexity of the tomographic
reconstructor is one of the problem for the future ELT-scale WFAO systems. The
sparse matrix technique can reduce the complexity. However, the phase covariancematrix Σϕϕ is a dense matrix and can not be incorporated into the sparse matrix
computation.
The sparse approximation is proposed by Ellerbroek [15], which assumes theturbulence follows κ11/3 ≈ κ4 . Under the approximation, the inverse of the phasecovariance matrix for ith turbulence layer can be rewritten by using Laplacian orcurvature operator as
2 TT(2.14)
Σ−1
ϕϕ,i ≈ Li Li = ci L̄i L̄iwhere L̄i is a discrete normalized Laplacian matrix, ci is a scaling constant for L̄i ,and Li = ci L̄i . The derivation of the scaling constant is given in Appendix A.2

Chapter 2. Methodology

32

based on [20]. The discrete Laplacian for all turbulence layers is a block diagonal
matrix similar to Σϕϕ .
While this approximated regularization matrix is fully sparse and can be incorporated into eﬃcient sparse-matrix techniques, it is known that the approximatedmatrix does not regularize tip-tilt modes and other modes that are curvaturefree[35]. The unregularized modes, mainly tip-tilt modes, may aﬀect the tomographic error especially in the case with LGSs, where tip-tilt modes can not bemeasured by LGS.

2.1.4

Multi Time-Step Tomographic Wave-Front Reconstruction

In this section, we will introduce our idea to reduce the tomographic error due
to the uncovered and unoverlapped areas by increasing information of the atmospheric turbulence using the measurements at both of the previous and thecurrent time-steps simultaneously. As mentioned in Chapter 1, the frozen flow
assumption allows us to consider the evolution of the atmospheric turbulence as
the movements of the turbulence layers due to the winds. Under this assumption,
it can be considered that the areas measured by WFSs at a previous time-step
shift due to the wind with time, which is indicated in Fig. 1.7. Thanks to the spatial displacement due to the wind, the areas corresponding to the measurements
at the previous time-step cover the uncovered and the unoverlapped areas at thecurrent time-step. An important point is using multi time-step measurements inthe tomography.
If an atmospheric turbulence layer moves with a constant speed v = (vx , vy ),the movement of the phase distortion pattern with a position x at ∆t previous
time-step is written by
ϕ(x, t − ∆t) = ϕ(x + v∆t, t).

(2.15)

Using Eq. (2.15), we can rewrite Eq. (2.2) and Eq. (2.4) into equations between theaperture-plane wave-front at the previous time-step t−∆t and the phase distortiondue to the turbulence layers at the current time-step t as()
hNa − hiϕi xφαj (x, t − ∆t) =+ hi αj + vi ∆t, t ,
hNa
i=1
∑

Nlayerφαj (t − ∆t) = Pα∆tϕ(t),j(2.16)(2.17)

Chapter 2. Methodology

33

where Pα∆t
is a ray-tracing matrix considering the movement of the phase distorj
tion pattern within the optical path footprint of LGS in direction αj during theduration time ∆t. Similarly, using Eq. (2.17) allows us to create a measurement
model connecting measurements at the previous time-step (t − ∆t) with the phase
distortion due to the atmospheric turbulence at current time-step t,sα (t) = Γα Pα∆t ϕ(t) + η(t)(2.18)

for j = 1, . . . , Nlgs .
where, Pα∆t is a concatenation of all submatrices Pα∆tj
In order to use the measurements from both of the current and previous timesteps simultaneously to increase information and to improve the performance ofthe tomography, the final measurement model is a concatenation of Eq. (2.10)
and Eq. (2.18). The reconstructor can be given as well as the single time-stepreconstructor.
−1
T
T
−1 T
T
Γ̃Tα(∆t) Σ̃−1Ẽβk (∆t) =Fβk (P̃α(∆t)ηη(∆t) Γ̃α(∆t) P̃α(∆t) + L L) P̃α(∆t) Γ̃α(∆t) Σ̃ηη(∆t)(2.19)=Fβk Ẽϕ(∆t)(2.20)

where
[

]
[
]
[
]
Pα
Γα 0
Σηη 0
P̃α(∆t) =, Γ̃α(∆t) =, Σ̃η(∆t) =
.
Pα∆t
0 Γα
0 Σηη

(2.21)

The wave-front in the science direction is derived from the measurements atthe current and the previous time-step as φ̂βk (t) = Ẽβk (∆t) s̃α(∆t) (t), in whichs̃α(∆t) (t) = [sα (t)T sα (t − ∆t)T ]T . We refer to this reconstructor as multi timestep tomographic reconstruction. Although the multi time-step reconstructorẼβk (∆t) uses measurements only from two time-steps, it is easily to expand thereconstructor to use multi time-step measurements more than two time-step, forexample, Ẽβk (∆t1 ,∆t2 ,...,∆tn ) . The wind speed and direction at each altitude shouldbe measured before the reconstruction. We will discuss the method to estimate
the wind speed and direction at each altitude in .

2.1.5

Time Diﬀerence ∆t

The time diﬀerence between the current and the previous time steps, ∆t, is an
important parameter for the multi time-step reconstructor. Although we assume

Chapter 2. Methodology

34

the frozen flow assumption, in reality, the time scale in which the frozen flow
assumption is valid is limited. In order to improve the performance of the tomographic reconstruction by the multi time-step reconstructor, ∆t should be chosen
to be shorter than the time scale within which the frozen flow assumption holds.
There are several studies to investigate the time scale of the frozen flow. Schöck
et al. [46] investigated the temporal decay of the intensity of the auto-correlationof the WFS slope measured on sky, and found that a value t50 , which is the times
for which the cumulative auto-correlation decreases to 50 % of its initial value,
is of order of 50–100 ms. If we use this time scale, t50 , as an index for ∆t, ∆tshould be equal to or shorter than t50 . However, it is noted that it is not clear
how the decay ratio of the auto-correlation and t50 relate to the performance ofthe multi time-step reconstructor. Also, the result from Schöck et al. consider
the total auto-correlation decay and they did not investigate the temporal decayof each turbulence screen. In order to understand the eﬀect of the frozen flow
assumption, more on-sky experiments are required. We will discuss the influenceof the frozen flow assumption on the multi time-step reconstructor in Chapter 4,again.

2.1.6

Iterative Computation for ELT-scale Tomography

In the case of WFAO systems on the current large telescopes, since we can computethe reconstructor matrix directly with inversing matrices of Eq. (2.11), the DMcommands are estimated by multiply the measured slope by the reconstructor.
However, on future ELT-scale WFAO system, computing Eq. (2.11) directly is
hard because the matrices are too large. Therefore, iterative method is used tocompute the tomographic reconstruction. In this thesis, we use the conjugate
gradient method (CGM) for the tomographic reconstruction on the ELT-scale
WFAO system.
The CGM solves an equation of Ax = b for x iteratively, where x is a column vector of unknown estimate, A is a coeﬃcient matrix, and b is a known
measurements. The coeﬃcient matrix A should be a positive definite symmetricmatrix. In the case of the tomographic wave-front reconstruction, we need the
CGM for both of the estimation step and the fitting step. In the estimation step,an equation to be solved by the CGM is given as−1−1(PαT ΓTα CηηΓα Pα + LT L)ϕ̂(t) = PαT ΓTα Cηηsα (t)(2.22)

where ϕ̂(t) is a column vector of the phase distortion of the atmospheric turbulent

Chapter 2. Methodology

35

layers and this is a x in this case. We use the sparse approximation for the phasecovariance matrix. Then, in the fitting step, we solve the following equation by
the CGM for the DM commands ûβk .
(NβTk Nβk )ûβk (t) = NβTk Pβk ϕ̂(t)(2.23)

The convergence of the CGM is decided by the relative error ||b − Axk ||/||b||,where xk is the estimates at k-th iteration.
Each iteration of the CGM consists of matrix-vector multiplications. Thanks
to the sparse approximation for the phase covariance matrix, all matrices in the
tomographic reconstruction are sparse. Therefore, we use sparse matrix computation for all matrix-vector multiplications in the CGM to reduce the computation
time of each iteration.
The other thing to consider for reducing computation time is reducing thenumber of iteration of the CGM. The number of the iteration of the CGM depends on the condition number of a coeﬃcient matrix A. The condition number,
κ(A), is used to measure how sensitive a equation is to error of measurements,b, and is defined as κ(A) = σmax (A)/σmin (A), where σmax (A) and σmin (A) are a
maximum and a minimum singular value of a matrix A, respectively. If condition
number is large, the estimate is aﬀected strongly by the measurement error and thenumber of the iteration required for the convergence increases. In general, before
computing the CGM, we multiply the equation by a matrix M , which reduces thecondition number of the matrix, so-called pre-condition matrix, as M Ax = M b or
M AM T (M T )−1 x = M b in order to reduce the number of iteration. It is notedthat M A or M AM T should be a symmetric matrix to use the CGM. Several
pre-conditioners are proposed for the CGM of the tomographic reconstruction. In
this thesis, we use a diagonal-scaling matrix as the pre-conditioner matrix, whichtransforms the diagonal elements of a matrix to 1. This is known as the simplest
pre-condition matrix.
Since the atmospheric turbulence may change smoothly with time, the estimate
at the one time-step before is a good initial condition for the CGM at current
time-step. This technique is called warm start, and can reduce the number of
iteration largely. We also use this technique out computation for future ELT-scale
tomographic reconstruction.

Chapter 2. Methodology

36
Star 2	

Star 1	

Atmosphericturbulent layer	Spatial cross-correlation map	

Telescope
aperture plane	

SH-WFS 1

SH-WFS 2	

Figure. 2.3: Principle of SLODAR. The lights from two star with an angular separation θcross at altitude h with a spatial displacement of hθ and propagate into [i, j]
subaperture of WFS1 and [i + δi, æ; δj] subaperture of WFS2, respectively.
The cross-correlation between these subapertures have a correlation power
proportional to the turbulence power at altitude h.

2.2

Pre-Estimation Method for Required Parameters

For implement the tomographic reconstruction, we need knowledge of the measurement noise, turbulence height and strength, this is, CN2 profile. In addition,wind speed and direction at each altitude are also required to perform the multitime-step reconstruction. Since these quantities vary with time, we should estimated these parameters and update the reconstructor based on the new estimatesin each minute. In this section, we describe the method to estimate the required
parameter from the measurement of multiple SH-WFSs.

2.2.1

SLODAR for Estimating the CN2 Profiles

Slope detection and ranging (SLODAR) is a well-studied method to estimates
the CN2 (h) profile with the measurement from two SH-WFSs. The SLODAR was
proposed originally by [50], and is updated by some studies [8, 12]. The basic
idea of SLODAR is shown in Fig. 2.3. Consider that there are two stars withan angular separation of θ and there is a turbulent layer at an altitude h. The
lights from the two stars cross at altitude h, and reach a telescope aperture witha spatial displacement of hθ. Therefore, if we observe these stars with diﬀerent
SH-WFSs and computes a spatial cross-correlation of the measurements from two
SH-WFSs, there is a peak in the spatial cross-correlation map, corresponding toChapter 2. Methodology

37

the spatial oﬀset of hθ. The intensity of the peak is proportional to the strengthof the turbulent layers, that is, CN2 value. In the case of multiple turbulent layers,
there are several peaks corresponding to each turbulent layers and measuring thelocation and intensity of each peak provides us the vertical profile of CN2 (h). In thisstudy, we use a method proposed by [8]. They determined the CN2 (h) by fitting thetheoretical spatial cross-correlation maps to the observed cross-correlation map.
We denote a measured slope in x direction of a subaperture indexed [i, j] at
x[1]
x[2]
time t for the first star as si,j (t) and the x slope for the second star as si+δi,j+δj (t),where δi and δj are an oﬀset of subapertures. In general, the CN2 estimation by
SLODAR is computed with using 1–3 minutes time series of the slope. We subtract
tip-tilt modes over the aperture from the measured slope to avoid the influence
from the vibration or tracking error. In the slope space, this is equal to subtract
the spatial average over subapertures from the measurement.
s′x[1]i,j (t) = si,j (t) −
x[1]

∑

1

nvalid valid i,jx[1](2.24)

si,j (t)where s′ is the tilt subtracted slope and nvalid is a number of the illuminatedsubaperture. We also subtract a temporal average in the slope measurement fromthe measured slope before computing the spatial cross-correlation. The average
subtracted slope are computed as′′x[1]i,js′x[1](t) = si,j (t) −

⟨

⟩

′x[1]
si,j (t)(2.25)

where s′′ is a slope after the tilt and average pattern subtraction and ⟨⟩ is atemporal average over 1–3 time series.
The spatial cross-correlation with subaperture oﬀset of [δi, δj] is defined as
⟨∑
obs
Cxx
(δi, δj) =⟩′′x[1]′′x[2]s(t)s(t)
i+δi,j+δjvalid i,j i,jN (δi, δj),(2.26)

where ⟨⟩ indicates averaging over time and N (δi, δj) is the number of the pair withthe separation [δi,δj]. The cross-correlation of y–y, x–y, and y–x can be defined
obs
obs
obs
similar to Eq. (2.26). However, since the cross-correlations, and Cyx
, Cxy
by Cyyobs
obs, have a very week correlation, we use only
and Cyx
between x and y slopes, Cxy
obs
obs
obs
obs
, can
and Cyyfor the SLODAR. The cross-correlation map, Cxx
and Cyy
Cxx
obs
obs
be obtained to compute Cxx and Cyy for all possible subaperture oﬀset (δi, δj).
The theoretical cross-correlation maps are computed for each altitude bin based

Chapter 2. Methodology

38

・・・	

・・・	

Figure. 2.4: Examples of theoretical slope cross-correlation. The outer scale is 30 m.
The only areas surrounding by the red lines are used for the fitting.

on the turbulence model, as shown in Fig. 2.4. The optimal altitude resolution dh
is given as dh = dsa /θ, namely the spatial oﬀset is an integral proporional to thesubaperture size dsa . In the case of WFAO with multiple GSs, we can compute the
SLODAR for all GS pairs. Since each GS pair has a diﬀerent angular separation θ,
the optimal resolution in the height direction is also diﬀerent. To avoid averaging
profiles with diﬀerent resolution, we set to a same altitude resolution for all GS
pairs. The maximum altitude measured by SLODAR also depends on the angularseparation of two stars, which is given by hmax = (nsa − 1)dsa /θ, where nsa is anumber of subaperture across a telescope pupil. Thus, smaller angular separation
gives us the information at higher altitudes with lower altitude resolution. It isnoted that the estimates of the CN2 value at higher altitude becomes noisier. The
high altitude turbulence makes the correlation peak at large spatial oﬀset [δi, δj].
The number of subaperture pair with the large oﬀset is limited compared with
the number of the subaperture pair with small oﬀset. In addition, if we haveGSs more than two, the high altitude is covered by GS pairs with large angularseparation. Therefore, the estimated CN2 value becomes statistically noisy. In thelater discussion, we remove the correlation between largest spatial oﬀset from thefitting to avoid the influence from the noisy estimates.
Here, we consider the case of npair pair of GSs. The theoretical spatial crosstheo[k]theo[k]correlation maps of altitude hi for kth GS pair are denoted as Cxx
and Cyy,
where the theoretical correlation maps are normalized so that CN2 (hi ) = 1. The
obs[k]
observed spatial cross-correlation maps for kth pair are represented as Cxx
andobs[k]
Cyy . If the turbulence at multiple altitudes are independent from each other,
the observed cross-correlation maps can be represented by a linear combination

Chapter 2. Methodology

39
theo[k]of the theoretical spatial correlation maps for diﬀerent altitudes as Cxx≈∑
obs[k]. The coeﬃcient ai for the normalized theoretical correlation map isi ai Cxx
2
equal to CN value at altitude hi . Considering all GS pairs, we can estimate CN2
profile by minimizing the following equation.
∑

npair

min
ai

k

])
[ theo[k]([ obs[k] ]
n
∑
C
(h
)
C
xx
xx
i
ai
Wk
theo[k]obs[k] −
(hi )
Cyy
Cyy
i

(2.27)

where W k is a masking matrix for the spatial cross-correlation for kth GS pair to
extracts the area along the direction of the GS baseline on the cross-correlation
map, which is the red area in Fig. 2.4. Since the peaks due to the turbulencelayers lie along the direction of the GS baseline, we use only the areas for the
SLODAR. The minimization of Eq. (2.27) is done with using non-negative leastsquare method [34] because CN2 value should be more than 0 at all altitudes. Thismethod provides better solution than the normal least square method without any
constraint.

2.2.2

Estimation of SH-WFS measurement noise

The measurement noise can be estimated by a method similar to SLODAR [12].
On , we use the auto-correlation of the slope measured by SH-WFSs, which isdefined as⟨∑⟩′′x[1]′′x[1]s(t)s(t)
i+δi,j+δji,j i,jAobs.
(2.28)xx (δi, δj) =N (δi, δj)
The auto-correlation has a peak at the center, [δi, δj] = [0, 0], as shown in theleft panels of Fig. 2.4. The central peak is enhanced by the measurement noise
since the noise has a correlation power with no spatial oﬀset for the auto correlation. Therefore, fitting the theoretical auto-correlations, Atheo
and Atheoxx
yy , toobs
obs
the observed auto-correlations, Axx and Ayy , with removing the central point,[δi, δj] = [0, 0], from the fitting and measuring the diﬀerence at the center provide
us an estimate of the noise averaged over all subapertures. This fitting also gives
us estimate of total turbulence of CN2 value including high altitude turbulence
which can not be measured by the SLODAR. Comparing the total CN2 value measured by the auto-correlation fitting with the value from the SLODAR provides
the residual turbulence power at high altitude not covered by SLODAR.
We note that this noise measurement gives us the averaged noise and can not
be applied for the LGS system since the noise of each subaperture is not same dueto the spot elongation and the outer subaperture has systematically a larger noise.

Chapter 2. Methodology

40

In order to optimise the tomographic reconstruction for the WFAO system withmultiple LGSs, the noise diﬀerence due to the spot elongation should be taken
into account, especially for future ELTs, in which the spot elongation will be more
severe.

2.2.3

Wind Speeds and Directions Estimation

In this section, we present the method to estimate wind speeds and directions
at multiple altitudes by using the result of tomographic reconstruction. If there
are Nl atmospheric layers, the tomographic reconstruction provides the phase
distortion on Nl atmospheric turbulence layers, ϕ̂(t) = [ϕ̂T1 (t) · · · ϕ̂TNl (t)]T . The
reconstructed phase distortion pattern on the atmospheric layers at each altitude
moves with time due to the wind. Tracking this movement allow us to estimate
wind speeds and directions at multiple altitudes separately. In order to detect thismovement due to wind, we use temporal correlation.
After reconstructing the phase distortion at each altitude by the tomographic
reconstruction, we extract aperture-size wavefront at each altitude in the directionof the center of a GS asterism from the reconstructed phase distortion because thesize of the turbulence layers are too large to compute temporal correlation. In
addition, it is better to use only areas reconstructed accurately for estimating
wind speeds and directions. Although the GS directions has smaller integrated
WFE or better SR than other direction because the tomographic reconstruction
is computed based on the measurements from the GS direction, it is possible thatthe most areas at high altitude in GS directions are covered by only own opticalpath, as shown in the left panel of Fig. 1.7, and aﬀected by the degeneracy due tothe unoverlapped area. This is because the reconstructed phase distortion at each
altitude in the GS direction can be cancelled each other to reduce total WFE inthis direction. Therefore, the reconstructed phase distortion in the GS directions
may not reflect the real phase distortion pattern, and using GS direction for thewind estimation is not optimal. As shown in Fig. 1.7, the direction of the centerof the GS asterism are covered by multiple GS footprints, and the reconstructed
phase distortion at each altitude is more accurate compared to the other directions.
We define a matrix Pci as the cropping matrix, which extracts the wave-front in
the direction of the center of the GS asterism from the estimated phase distortionϕ̂i of the ith turbulent layer.
ϕ̂ic (t) = Pci ϕ̂i (t).

(2.29)

Chapter 2. Methodology

41

where ϕ̂ic (t) is the wave-front at time t extracted from ith turbulence layer. Before
computing the temporal correlation, the extracted wave-front is converted into theslope by a discrete phase-to-slope operator Γc .
ŝic (t) = Γc ϕ̂ic (t),(2.30)

where ŝic (t) is a reconstructed slope vector corresponding to a phase distortion ofi-th atmospheric turbulence layer in the direction of the center of the GS asterism
at time t. A peak of a slope correlation is sharper than the peak of a wavefront
correlation, and can be detected easier than the peak of the wavefront correlation.
This is because the slope measured by SH-WFS corresponds to a diﬀerential of thephase distortion and the diﬀerential play a role of a high-pass filter. Therefore,
the small-scale and weak phase distortion is more weighted than large-scale strong
distortion by converting to slope, and the movement of the distortion pattern ismade clearer on the temporal correlation.
Temporal correlation of the reconstructed slope at i-th turbulence layer can becomputed by the same way as SLODAR.
⟨∑
obs[i]
Txx(δp, δq, δt) =

⟩

′′x,i
′′x,i
p,q ŝp,q (t)ŝp+δp,q+δq (t + δt)

N (δp, δq).

(2.31)

where δt is the time delay for the temporal correlation. The temporal correlation
obs[i]
map Txx can be obtained by computing Eq. (2.31)for all possible pairs of δp andδq. The correlation for y-direction slope is represented as similar to Eq. (2.31). In
order to increase signal to noise, the correlations of x and y-direction slopes are
obs[i]
obs[i]
averaged, T obs[i] = (Txx + Tyy ).
The movement of the correlation peak on the correlation map T obs[i] with δtcorresponds to the wind speed and direction of ith turbulence layers. We use thecentroid algorithm to detect the peak. When a subaperture size of SH-WFS is dsa
and the peak position is (δp, δq), the wind speed is computed by vx = dsa δu/δt
and vy = dsa δv/δt. By computing vx and vy and averaging it over diﬀerent δt, thewind speed and direction at each altitude can be estimated.

Chapter 3
Analytical Evaluation
3.1

Analytical Tomographic Error

The tomographic error in direction βk is defined as σβk = (ϵβk tomo /np )1/2 , where
ϵβk is the residual phase variance at an aperture plane defined in Eq. (2.6) and
np is the number of phase point at aperture plane. Here, we assume that the allsubaperture has a same measurement noise variance ση2 . In this case, the covariancematrix of the measurement noise Σηη become ση2 I, where I is an identity matrix.
We also assume that the DM interaction matrix Nβk for βk direction is an identity
matrix for simplicity, which means the DM command is equal to the phase causedby the DM. Using this assumption and Eq. (A.2), the tomographic error can be
represented asσβ2 k ,tomo =[]T
Tr Rpiston (Pβk − Eβk Γα Pα )Σϕϕ (Pβk − Eβk Γα Pα )T Rpiston
np
[
]
T
T
Tr Rpiston Eβk Eβk Rpiston
2
+ ση
np

=σβ2 k ,geo + σβ2 k ,noise .

(3.1)

where Rpiston is a matrix removing a piston mode from the residual wave-front. If
you would like to compute the tomographic error without piston and tilt modes,
it is achieved by the replacing Rpiston to Rtip-tilt Rpiston , where Rtip-tilt is a matrix
removing tip-tilt. This tomographic error can be computed by the matrix calculation, so-called analytical tomographic error. The first term σβ2 k ,geo represents the
tomographic error determined by the geometry such as GS configuration, turbulence structure, and reconstructor, and we refer to σβ2 k ,geo as geometric error. The
influence from the uncovered and the unoverlapped areas are included in the firstterm. The second term σβ2 k ,noise is a noise propagation of the measurement noiseση . The purpose of the multi time-step reconstructor is to reduce the geometric42

Chapter 3. Analytical Evaluation

43

Parameters

ValuesDiamter of an telescope aperture

30 m

Fried parameter r00.156 mouter scale L0

30 m

Altitudes of the turbulence layers

0 km, 5 km, and 10 km

Fraction of CN2

50, 25, and 25 %Wind speeds

5 m s−1 , 10 m s−1 , and 20 m s−1

Wind directions

90◦ , 0◦ , 0◦

Number of NGSs

3

NGS asterism radius rast

20–200 arcseconds

WFS

Shack-Hartmann WFS
with 30×30 subapertures

Diameter of subaperture

1m

Measurement noise

0.01 arcsecond

Tab. 3.1: Parameters used for the analytical evaluation.

error.

3.2

Model

The model for the analytical evaluation is as follows. We consider a simple MOAO
model on a 30 m circular aperture telescope with three NGSs. The NGSs are put
on a ring of a radius rast , which is a range of 20–200 arcseconds, indicated as the redplots in Fig. 3.2. We evaluate the tomographic error in the direction of the centerof FoR and in the outer directions, shown as the blue and green plots in Fig. 3.2.
Fig. 3.1 shows a footprints of the NGS optical paths and the evaluated directions at
10 km as rast = 20 arcseconds (the left panel) and rast = 200 arcseconds (the rightpanel). Fried parameter r0 =0.156 m and the outer scale is 30 m, which are typicalvalue at Maunakea. The atmospheric turbulence consists of three turbulence layers
at 0 km, 5 km and 10 km. The CN2 fractions are [50 %, 25 %, 25 %], the windspeeds are [5 m s−1 , 10 m s−1 , 20 m s−1 ], and the wind directions are [90◦ , 0◦ , 0◦ ]at 0, 5, and 10 km, respectively. The size of subaperture of SH-WFSs is 1 m on
the aperture-plane and all subapertures have a same measurement noise with avariance of ση2 . The spatial sampling of discrete grids on the turbulence layers andthe aperture-plane is 1 m. The parameters used for this analysis are summarized
Table 3.1. Over this chapter, the temporal lag error and the fitting error are notChapter 3. Analytical Evaluation

Footprints of NGS optial paths
Footprint of the central direction
Footprints of the outer directions

30
20

20

10

10

0

0

−10

−10

−20

−20

−30

−30
−30 −20 −10

0
10
X [m]

20

Footprints of NGS optial paths
Footprint of the central direction
Footprints of the outer directions

30

Y [m]

Y [m]

44

30

−30 −20 −10

0
10
X [m]

20

30

Y [arcseconds]

Figure. 3.1: Footprints of the NGS optical paths (red) and the evaluated central direction (blue) at altitude of 10 km for rast =20 arcseconds (left) and rast =200
arcseconds (right).

200

NGSs	

100

Central evaluateddirection	

0Outer evaluateddirections	

−100
−200
−200

−100
0
100
X [arcseconds]

200

Figure. 3.2: The NGS asterism (red) and the evaluated directions. The blue plot is
the central evaluated direction. The green plots are the outside evaluateddirections.

Chapter 3. Analytical Evaluation

45

included in this analysis and we assume that we have the perfect knowledge of the
turbulence altitudes, powers, and wind speeds and directions.

3.3

Eﬀect of Lack of Information

First we demonstrate the influence from uncovered and unoverlapped area. Fig.
3.3 show the geometric error σgeo , the noise propagation σnoise , and the total tomographic error σtomo as a function of the asterism radius rast . The red solid lines
in Fig. 3.3 indicate the result with the errors of the single time-step reconstructor
evaluated in the central direction and the red dotted lines are the errors of the single time-step reconstructor averaged the values evaluated in the outer directions.
Those errors are computed with removing tip-tilt errors.
When rast is small, the geometrical error of the central direction (solid red line)
is small because the central direction is covered densely by the multiple footprintsof the NGS optical path in the atmosphere as the left panel of Fig. 3.1. In the outerdirections, the geometric error (dashed red line) is large due to the uncovered area
at small rast . Then, as the asterism radius increases, the geometric error in thecentral direction also increases rapidly. At rast = 200 arcseconds, the geometricerror in the central direction is 538 nm, which is ten times larger than the value
at rast = 20 arcseconds. This is due to the unoverlapped area. At rast = 200
arcseconds, all area of the footprint in the central direction is covered by the NGS
optical paths, as shown in the right panel of Fig. 3.1. However, some parts offootprints in the central direction are covered by only one footprint of the NGS.
The unoverlapped area increases with the asterism radius and results in the large
geometric error. Although the geometric error in the outer directions decreases
because of decreasing the uncovered area, the geometric error is still large, which
is 408 nm, since most of the footprints in the outer directions is unoverlapped
area. The noise propagation, which are indicated by the red lines in Fig. 3.3(b),is much smaller than the geometric error for both directions. Therefore, the total
tomographic error of the single time-step reconstructor in Fig. 3.3(c) is dominatedby the geometric error.
This result suggests the diﬃculty of expanding the FoR size of the tomographic
WFAO systems. In order to expand the FoR size, we have to measure and reconstruct a large volume of the atmospheric turbulence corresponding to the large
FoR by expanding the NGS asterism. However, it involves the large tomographic
error as we show. This eﬀect become more severe with using LGSs because ofChapter 3. Analytical Evaluation

46

700
Single,Center

Single,Outer

Multi,Center

Multi,Outer

Tip-tilt removed σgeo [nm]

600
500
400
300
200
100
0
20

40

60

80

100

120

140

160

180

200

rast [arcseconds](a) Tip-tilt removed geometric error σgeo

Tip-tilt removed σnoise [nm]350Single,Center

Single,Outer

Multi,Center

Multi,Outer

300
250
200
150
100
50
0
20

40

60

80

100

120

140

160

180

200

rast [arcseconds](b) Tip-tilt removed noise propagation σnoise

Tip-tilt removed σtotal [nm]700Single,Center

Single,Outer

Multi,Center

Multi,Outer

600
500
400
300
200
100
0
20

40

60

80

100

120

140

160

180

200

rast [arcseconds](c) Tip-tilt removed tomographic error σtomo

Figure. 3.3: The analytical geomatric error σgeo (top), the noise propagation σnoise (middle), and the total tomgoraphic error σtomo (bottom) as a function of theasterism radius rast . The all values are computed with removing tip-tilt
modes from the residual wave-fronts. The red lines show the result withthe single time-step tomographic reconstruction and the blue lines are theresults of the multi time-step reconstruction. The solid lines indicate thevalue evaluated in the direction of the center of the FoR. The dashed lines
are the results averaged over outer directions.

Chapter 3. Analytical Evaluation47

its smaller footprints at high altitudes due to the cone eﬀect. This is a criticalproblem at future ELT-scale WFAO system.

3.4

Analytical Performance of Multi Time-Step Reconstruction

The blue lines in Fig. 3.3 show the results of the multi time-step reconstructor
with using the measurement from two time-steps. We set the time diﬀerence ∆t
for the previous time-step to 50 ms. This value is determined by considering
temporal de-correlation of the atmospheric turbulence. In [46], they investigated
the temporal decay of the auto-correlation intensity of the slope measured by a
SH-WFS and defined the time t50 , which is the times for which the intensity of theauto-correlation decreases to 50 % of its initial value. In Fig. 3.3(a), the diﬀerenceof the geometric error in the central direction between the single and the multitime-step reconstruction is small as the asterism radius is small. This meansthat the single time-step reconstructor is not aﬀected by the unoverlapped area inthe central direction with the small asterism radius. In the outer directions, thegeometric error of the multi time-step is smaller than the error of the single timestep reconstruction for the outer directions, but is much larger than the geometricerror of the multi time-step reconstruction for the central direction. Here, we use
the measurement at 50 ms previous time-step. The spatial displacement due tothe wind between the current and the previous time-steps is 0.25 m, 0.5 m, and
1.0 m at 0 km, 5 km, and 10 km, respectively. Since these oﬀset is not enoughto fill the uncovered area in the outer directions, the geometric error in the outerdirections are not improved so much by the multi time-step reconstruction.
As the asterism radius large, the diﬀerence of σgeo can be seen clearly. The
multi time-step reconstructor can reduce dramatically the geometric error in boththe central and outer directions. The geometric error in the outer directions arecomparable to the value in the central direction. This result demonstrates thatusing the measurements from multi time-steps is eﬀective to solve the degeneracy
due to the unoverlapped area. The σgeo of the multi time-step reconstructor is
158 nm at rast =200 arcseconds and is much smaller than the value of the single
reconstructor.
The noise propagations of the multi time-step reconstruction is a bit larger than
the noise propagation of the single time-step reconstructor in both the central andouter directions. At large asterism radius, the error due to the noise propagation

Chapter 3. Analytical Evaluation48

is roughly comparable to or larger than the geometric error. As result, the total
tomographic error is aﬀected by the both eﬀects, as shown in Fig. 3.3(c). However, the tomographic error is much smaller than the value of the single time-stepreconstructor. As the asterism radius is 200 arcseconds, the multi time-step reconstructor achieves the tomographic error σgeo ∼220 nm. This is corresponding to
the tomographic error in the central direction of the single time-step reconstructor
at the asterism radius of rast =80 arcseconds. Although this result is based on the
optimal condition, where the turbulence altitudes, powers, and wind speeds anddirections are given, the multi time-step reconstructor has a possibility to reduce
the geometric error and expand the FoR size dramatically.

3.5

Dependency on Time Diﬀerence ∆t

Fig. 3.4 shows the geometric error σgeo , the noise σnoise , and the total tomographic
error σtomo as a function of the time diﬀerence ∆t for the multi time-step reconstructor. Those errors are also computed with removing tip-tilt errors and we use
two time-step measurements for the multi time-step reconstruction. The asterism
radius is assumed 200 arcseconds and therefore the main problem is the degeneracy due to the unpverlapped area. The used colors are same as Fig. 3.3. The
horizontal red lines in Fig. 3.4 are results from the single time-step reconstructor,
which doesn’t depend on the time diﬀerence ∆t.
In Fig. 3.4(a), an interesting point is that the multi time-step reconstructor
can reduce the geometric error compared to the single time-step reconstructor
even with very small ∆t such as 10 ms. The spatial displacement due to the windduring 10 ms are 0.05 m, 0.1 m, and 0.2 m at 0 km, 5 km, and 10 km, respectively.
These displacement is much smaller than the size of subaperture of 1 m. While thegeometric error of the multi time-step error decreasing with ∆t, it looks saturated
as ∆t =60–80 ms, where the spatial displacements are 0.4 m, 0.8 m, 1.6 m ateach altitude, respectively. This results suggests that the spatial displacement
equal to the subaperture size may be enough to solve the degeneracy due to theunoverlapped area.
With respect to the noise propagation, the multi time-step reconstructor isaﬀected more by the measurement noise compared with the single reconstructor
as shown in Fig. 3.4(b). The WFE due to the noise propagation has a similar trend
to the Fig. 3.4(a) and is comparable to the WFE due to the geometric error.
As a result, the total tomographic error of the multi time-step reconstruction isChapter 3. Analytical Evaluation49

aﬀected by both σgeo and σnoise . However, it is much smaller than the tomographic
error of the single time-step reconstructor. This diﬀerence of σtomo becomes larger
with ∆t.

Chapter 3. Analytical Evaluation

50

700

Tip-tilt removed σgeo [nm]600Single, Center

Single, Outer

Multi, Center

Multi, Outer

500
400
300
200
100
0
20

40

60

80

100

120

140

160

180

∆t [ms](a) Tip-tilt removed geometric error σgeo

Tip-tilt removed σnoise [nm]

350
300

Single, Center

Single, Outer

Multi, Center

Multi, Outer

250
200
150
100
50
0
20

40

60

80

100

120

140

160

180

∆t [ms](b) Tip-tilt removed noise propagation σnoise

Tip-tilt removed σtomo [nm]

700
600

Single, Center

Single, Outer

Multi, Center

Multi, Outer

500
400
300
200
100
0
20

40

60

80

100

120

140

160

180

∆t [ms](c) Tip-tilt removed tomographic error σtomo

Figure. 3.4: The analytical geomatric error σgeo (top), the noise propagation σnoise (middle), and the total tomgoraphic error σtomo (bottom) as a function of theasterism radius rast . The all values are computed with removing tip-tilt
modes from the residual wave-fronts. The color and line type are same as
Eq. (3.3)

Chapter 4
Numerical Simulation
4.1

Simulation Setting

In this chapter, we show the results based on the numerical simulation with assuming ELT-scale MOAO system and investigate the expected performance of themulti time-step reconstruction on these future ELT. We assume a circular aperture
with 30 m diameter. The MOAO system has a FoR of 10 arcminutes diameterand does the tomographic reconstruction with eight sodium LGSs at altitude of
88 km. Two LGS asterisms are used in this simulation, which are shown in Fig.
4.1. One is a narrow asterism indicated as open squares in Fig. 4.1, and the otheris a wide asterism indicated as filled squares. As mentioned in Chapter 1, tip-tilt
modes in wave-front can not be sensed by the LGS and therefore, TT-NGSs are
required in a FoR to measure the tip-tilt modes. In this simulation, however, we
assume that we can measure the low-order modes of the phase distortion fromthe LGSs for simplicity. The cone eﬀect and the spot elongation on SH-WFS areconsidered in the simulation. The vertical sodium-layer profile is approximated toreproduce the lidar measurements [19] as

Na(h) =)({
(h−hNA )2(for |h − hLGS | < σNa
exp − 2σ2Na(for |h − hLGS | > σNa

0

(4.1)

where hNA = 88 km is an altitude of a sodium layer from the telescope, and σNa =5
km is a half width of the sodium-layer. The vertical profile of the sodium layer
is presented in Fig. 4.2. We use SH-WFSs with 60×60 subapertures and DMs
for science targets with 60×60 elements. The spatial sampling of atmosphericturbulence is 0.05 m. The WFE is also evaluated with this spatial sampling tosimulate the fitting error. We assume that the MOAO system is controlled with800 Hz and there is two frames delay between the measurement by WFSs andcorrection by DMs, that is, 2.5 ms delay on a frame rate of 800 Hz.
51

Chapter 4. Numerical Simulation

52

300

Y [arcsecond]

200
100
0
−100
−200
−300
−300 −200 −100 0
100
X [arcsecond]

200

300

Figure. 4.1: LGS asterisms used in the numerical simulation. A narrow asterism is indicated with the red filled squares. The blue open squares show a wide
asterism. The black solid lines show the 2.5 arcminutes and 5 arcminutes
radius from the center of a FoR. The dashed lines are directions where theperformance is evaluated.

1
0.9
Nomalized intensity

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
80

82

84

86

88

90

92

94

96

Altitude [km]Figure. 4.2: Vertical profile of normalized intensity of Sodium layer used in the simulation. This is approximated based on the lidar on-sky measurements [19].
The altitude and of sodium layer are 88 km and 10 km.

We assume the top of Maunakea as the observation site in the simulation anduse a seven-layer model used in [4] for the CN2 profile. This model is created basedon image-quality measurements from the Subaru Observatory [39] and diﬀerential
image motion monitor (DIMM) and multiple aperture scintillation sensor (MASS)
measurements by Thirty Meter Telescope Project at Maunakea [16]. The modeltakes into account for additional dome seeing in the Subaru telescope, thus theground layer has a strong turbulence power, which is 60 % of the total turbulence
power. Fried parameter r0 is 50 %ile value of the seeing measurements, r0 = 0.156

Chapter 4. Numerical Simulation

53

m, and an outer scale is assumed to be 30 m. The CN2 values at each altitude
are summarized in Table 4.2. We assume a gaussian model for the wind speeds oneach layer based on [24]. The assumed wind speeds and directions are also listed
in Table 4.2. The performance is evaluated at wavelength of 1650 nm (H-band).
For the computation of the tomographic reconstruction, we use the iterative
CGM with the simple diagonal scaling preconditioner. Also, we use Tesla K40
GPU for the parallel computation of the tomographic reconstruction. The codes
for GPGPU are made with CUDA 6.0, which is C++ like language for GPGPU.
However, all results shown in this chapter are computed with assuming that thecomputation time is negligible. With respect to the computational eﬃciency withGPGPU, we will discuss it at end of this chapter.

4.2

Simulation Results

4.2.1

Performance of Multi Time-step Wave-Front Reconstruction

First, we simulated an optimal performance of the single and the multi timestep reconstructors for both of the narrow and the wide LGS asterisms with anassumption that the turbulence profile, wind speeds and directions are knowna priori. The results are illustrated in Fig. 4.3. SR values are measured fromsimulated point spread function (PSF) images by comparing the peak intensitiesof the PSF images to the peak intensity of a diﬀraction-limited PSF image createdby the simulation. All PSF images from the simulation are normalized by the total
intensity with in a 1 arcsecond box. The top and the bottom panels in Fig. 4.3
show the simulated maps of the Strehl ratio (SR) across the FoR of 10 arcminutes
diameter with the single (top) and the multi time-step (bottom) reconstructors.
We use ∆t of 50 ms for the multi time-step reconstructor as well as Chapter 3. The
bottom panels in Fig. 4.3 show SR profiles as a function of an angular distance fromthe center of the FoR. The profiles are computed by averaging over six directions
shown as dashed lines in Fig. 4.1.
With the narrow asterism, the single time-step reconstructor, which is indicated as red filled squares in the left panels of Fig. 4.3, achieves the averaged SR
values of ∼0.5 within the inner LGS radius at 75 arcseconds from the center. The
averaged SR value decreases slowly with an angular separation from 75 arcseconds
to 150 arcseconds because of increasing the unoverlapped area as shown in the left

Chapter 4. Numerical Simulation

54

Parameters

ValuesDiameter of aperture

30 m

FoR

10 arcminutes

Zenith angle

0◦

Number of LGSs

8

Height of LGSs

90 km

Number of turbulence layers

7

r0 at 500 nm

0.156 m

L0

30 m

WFS

SH-WFS

Number of WFS subapertures

60×60

Number of DM elements

60×60

Frame rate

800 HzTime lag2 frames

count per a subaperture

700 electrons

readout noise

3 electrons

evaluate wavelength

1650 nm (H-band)

Tab. 4.1: Parameters used in the numerical simulation.

Altitude∫

CN2 dhFraction wind speed

wind direction[km][10−14 m1/3 ]of CN2[m s−1 ][degree]

16

2.734

0.0826

7.0

0

8

2.264

0.0684

33.0

45

4

2.879

0.0869

19.7

90

2

1.233

0.0372

11.6

135

1

1.074

0.0325

9.0

180

0.5

3.190

0.0963

8.0

225

0

19.7370.59607.0

270

Total r0 [m]0.156

Tab. 4.2: Atmospheric turbulence profile used in the numerical simulation.

Chapter 4. Numerical Simulation

0.
15

0.2
5

0.25

0.4

Strehl Ratio

0.
05

0.2
0.15

0.0
5

0.1

−200

0.25

0.0
2

0.2

0.2
5

0.2

0.5

0.3
−100
25
0.

−200

0.35

5
0.0

2
0.0

0.1
5
0.0

0.3
5

5
0.3

0.4

0

5
0.1

0.5
0. 0
2 .3

0.3

100

0.2
5

0.4

5
0.0

0.3

0.5

−100

0.6
0.1
5

0.5
0.4

0.7

0.0
5

0.2
5

200
0.25

0.1

0.5

0

300

0.6

0.4

100

0.1
5
0.0

Y [arcsecond]

200

0.7

0.02

Y [arcsecond]

2
0.0
0.05
0.1
0.2
0.3

Strehl Ratio

300

55

0.1

0.15

0
100

200

0.0
4

100

200

0.6

0.4

0.48
38
0.

0.2

−100

0.1

−200

0

−300
−300 −200 −100

0.3

0.18

0.2
8

0.5

0.38

38
0.
0.4
8

0.2
8

0

0.7

0.
18

0.2

0.48
0.3
8

300

0.28

0.38

0.18

0

Strehl Ratio

0.
08

Y [arcsecond]100

X [arcsecond]

0.1
0

100

200

300

200

250

X [arcsecond]

0.7

0.7

0.6

0.6

0.5

0.5
Strehl Ratio

Strehl Ratio

300

0.2
8

8
0.0

0

Strehl Ratio

0.32
0.2
2

0.3

2
0.1

−300
−300 −200 −100

200

0.38

0.38

0.52

0.
04

200

0.5
0.4

8
0.1

0.18

0.3
2 0.42

2
0.1

22
0.

−200

0.6

2
0.6

0
−100

0.12
0.22
2
0.2 .32
0
2
0.
4
52
0. 2
0.5

300

0.7

0.0
4

42
0.

100

2
0.1

Y [arcsecond]

200

100

08
0.

4
0.0

0

X [arcsecond]

X [arcsecond]
300

0

−300
−300 −200 −100

300

0.
08

0

0.
28

−300
−300 −200 −100

0.4
0.3
0.2

0.4
0.3
0.2

Single0.1

Single

0.1

Multi

Multi

0

0
0

50

100

150

200

Angular sparation [arcsecond]

250

300

0

50

100

150

300

Angular sparation [arcsecond]

Figure. 4.3: Top and middle : Simulated SR maps within the 10 arcminutes FoR computed by using the single time-step reconstructor (top) or the multi timestep reconstructor with ∆t=50 ms (middle). The right panels are the resultswith the narrow LGS asterism, and the left panels are computed with thewide LGS asterism. The multi time-step reconstructor improves SR values
over the FoR for both of the asterisms. Bottom : SR profiles as a function ofan angular separation from the center direction of the FoR. The profiles areaveraged over directions shown as the dashed lines in Fig. 4.1. The red filledsquares show the averaged SR profile with the single time-step reconstructor,
and the blue open squares is the result by the multi time-step reconstructor
with ∆t=50 ms. The dashed black lines show the radii of LGS positions.

56

40

40

30

30

20

20

10

10
Y [m]

Y [m]

Chapter 4. Numerical Simulation

0

0

−10

−10

−20

−20

−30

−30

−40

−40
−40 −30 −20 −10 0 10
X [m]

20 30

40

−40 −30 −20 −10 0 10
X [m]

20 30

40

Figure. 4.4: The footprints of the LGS optial paths on the highest turbulence screen
at altitude of 16 km The left panal shows the footprints with the narrowasterism and the right is for the wide asterism. The dashed black circle
indicates the radius of 150 arcseconds and 300 arcseconds from the centerof the FoR.

panel of Eq. (4.4). In the outer region than 150 arcseconds radius, the averaged
SR values decrease steeply because there is no LGS and these areas are aﬀected byuncovered area. The multi time-step reconstructor with ∆t =50 ms, represented as
the blue open squares, can increase the SR value over the FoR. The improvementof the averaged SR value is approximately 0.1 at 150 arcseconds, even though
there is no uncovered area at r ≤ 150 arcseconds. This suggests that the single
time-step reconstructor is aﬀected by the degeneracy due to the unoverlapped areaeven with the narrow asterism.
In the case of the wide asterism, illustrated in the right panels of Fig. 4.3, theSR value at the outer region than 200 arcseconds radius is better than the valueof the narrow asterism because the wider area is covered by the LGSs as shown in
the right panel of Eq. (4.4). In the contrast, the larger separation between LGSs
causes a larger unoverlapped area at high altitudes. Hence, the SR value at inside
area becomes smaller than the narrow asterism, and there are valley areas between
LGSs on the SR map (the left top panel of Fig. 4.3). From the bottom right panelof Fig. 4.3, the averaged SR value is 0.24 with the single time-step reconstructor
at maximum. This is aﬀected strongly by unoverlapped area. Using the multitime-step reconstructor with ∆t = 50 ms can double the averaged SR over the
FoR for the wide asterism. At 300 arcseconds from the center, the SR value fromthe multi time-step reconstructor with the wide asterism is 3 times larger than
the result of the narrow asterism.

Chapter 4. Numerical Simulation

4.2.2

57

Dependency on Time Diﬀerence ∆t and Frozen Flow
Assumption

The performance of the multi time-step reconstructor depends on the spatial displacement between the areas measured by the current and the previous time-steps
at each altitude, which we refer as d(h) in the paper. This displacement d(h)corresponds to the movement of the atmospheric turbulence layer at an altitude
h during ∆t, d(h) = v(h)∆t, if the wind speed is v(h) at an altitude h. The
dependence of the multi time-step reconstructor on the ∆t is presented in Fig.
4.5. The vertical axis in Fig. 4.5 is the SR improvement ratio, kSR , which represents the ratio of the SR value achieved by the multi time-step reconstructor
to the SR value from the single time-step reconstructor and is averaged over the
angular separation from the center of the FoV. The multi time-step reconstructor
achieves larger kSR for the wide asterism than the narrow one. This suggests thewide asterism is aﬀected more by the unoverlapped area than the narrow asterism.
However, the trend of the dependence of the kSR on ∆t looks similar for the both
asterisms. This suggests that the dependence on ∆t doesn’t depend the asterismof GSs. The multi time-step reconstructor improves the SR even with small ∆t.
There is 23 % increase in the SR value as ∆t=20 ms for the narrow asterism.
The improvement ratio kSR is maximized at ∆t ∼ 100 ms, where kSR ∼1.8 and∼2.1 for the narrow and the wide asterism, respectively. The displacements d with
∆t=100 ms are 3.3 m and 0.7 m for the fastest and slowest layer in the model.
While, as ∆t is larger than 100 ms, kSR is almost no or very week dependency on∆t. This suggests that only small displacement between the areas measured at
the current and the previous time-step is enough to solve the degeneracy due tothe unoverlapped area.
As mentioned before, the time scale of the frozen flow assumption is valid
only within a short time. After the time scale of the frozen flow assumption,
the phase distortion pattern of the atmospheric turbulence varies. If the timescale of the frozen flow is smaller than ∆t for the multi time-step reconstructor,
the measurement at previous time step can not be used for the information of thecurrent time step because the phase distortion pattern is no longer kept. Guesalaga
et al. [21] also investigated the time evolution of the atmospheric turbulence byusing spatio-temporal cross-correlations of the measurements from multiple SHWFSs installed in the GeMS, which is a MCAO system on the Gemini South.
They found that the decay ratio of the correlation intensity for an individual
layer, fdecay , decreases with the time delay for the correlation, ∆t, and depends

Chapter 4. Numerical Simulation

58

2.6

t70

2.2

t80

t90

2.4

kSR

2
1.8
1.6
1.4
1.2

narrow

wide

1
0

50

100

150

200

250

300

∆t [ms]Figure. 4.5: Dependence of the multi time-step reconstructor on the ∆t is presented inFig. 4.5. The vertical axis in Fig. 4.5 is the SR improvement ratio, kSR ,which represents the ratio of the SR value achieved by the multi time-stepreconstructor to the SR value from the single time-step reconstructor and isaveraged over the angular separation from the center of the FoV. The resultwith the narrow asterism is indicated as the red filled squares, and the blue
open squares shows the result with the wide asterism. The vertical dashedlines indicate the time duration that the decay ratio fdecay of the temporalcorrelation of SH-WFS measurements becomes 90 %, 80 %, and 70 %, thatis, t90 , t80 , and t70 , wind assuming the wind speed is 33 m s−1 , which is
maximum value in the wind model.

on the distance travelled by the layer, dtravel . In their paper, the decay ratio areapproximated as a function of ∆t and v.
fdecay = (−0.157v − 0.365)∆t + 1(4.2)

By using this equation, we can estimate ∆t to keep the decay ratio fdecay to 90 %,80 %, and 70 %, that is, t90 , t80 , and t70 , respectively. The dashed vertical lines inFig. 4.5 show t90 , t80 , and t70 with assuming the fastest wind speed in the model,
v = 33 m s−1 . The time t70 is 54 ms for the fastest layer and 200 ms for the slowest
layer with wind speeds of 7 m s−1 . In order to keep the decay ratio of more then
70 % for all altitudes, ∆t should be equal to or shorter than 54 ms. It is noted
again that the eﬀect of the decay ratio on the performance of the multi time-stepreconstructor is still not quantitative. On-sky experiments and validations arenecessary to evaluate the influence of the frozen flow assumption.

Chapter 4. Numerical Simulation

59

2

kSR

1.5

1

0.5
ek , narrowe⊥ , narrow

0
−10

−5

ek , widee⊥ , wide
0

5

10

Wind Noise ek or e⊥ [m s−1 ]Figure. 4.6: SR improvement ratio, kSR , achieved by multi time-step tomographic reconstruction with diﬀerent wind speed errors e∥ (the filled symbols) or winddirection errors e⊥ (the open symbols). The results with the narrow asterism are represented as the red squares and the results with the wide asterism
are the blue circles. The SR ratio less than 1 means that the performance ofthe multi time-step reconstructor is poorer than the performance achievedby the single reconstructor due to the wind estimation error.

4.2.3

Eﬀect of Uncertainty of Wind Speeds and Direction

We discuss an eﬀect from uncertainties of wind speed and direction. The errorsof wind speed and direction result in an error of the spatial displacement between
the areas measured by the current and the previous time-step. Here, we represent
the wind error as e = (e∥ , e⊥ ). The first one, e∥ , is an error parallel to thewind direction and referred to as wind speed error. The second, e⊥ , is an error
perpendicular to the wind direction and referred to as wind direction error. Fig.
4.6 shows the improvement ratio kSR achieved by the multi time-step reconstructor
comparing to the single time-step reconstructor with diﬀerent errors of wind speed
and direction, with both of the narrow (red symbols) and wide (blue symbols)
asterisms. We assume the same wind error for all atmospheric layers and ∆t=50
ms for the multi time-step reconstructor. The curve of kSR with the wind direction
error e⊥ is symmetric with respect to the point e⊥ = 0. The multi time-stepreconstructor has the advantage comparing to the singe time-step reconstructor,
kSR > 0, as −7 m s−1 ≤ e⊥ ≤ 7 m s−1 for the condition we use for both asterisms.
This corresponds that the error of the spatial displacement is less than 0.35 m
and that the error of angle of wind direction is less than 12 degrees for the fastest
layers with 33 m s−1 wind and 45 degrees for the slowest layer with 7 m s−1 wind.

Chapter 4. Numerical Simulation
5 arcminute FoR

0.8

Point, Single
Point, Multi
Extended, Single
Extended, Multi

0.7

Ensuqured Energy

0.6

10 arcminute FoR

0.8

Point, Single
Point, Multi
Extended, Single
Extended, Multi

0.7

Ensuqured Energy

60

0.5
0.4
0.3
0.2
0.1

0.6
0.5
0.4
0.3
0.2
0.1

0

0
0

50

100

150

200

250

Angular Separation [arcsecond]

300

0

50

100

150

200

250

300

Angular Separation [arcsecond]

Figure. 4.7: Ensqured energy profiles with the narrow (left) and wide (right) asterism.
As the target intensity profile, the point source and extended source areconsidered. The extended source is assumed to be an exponential profile
with an eﬀective radius of 0.1 arcsecond. The EE values are measured fromthe simulated PSF images (red) with a 50 mas box for the point source,
and the exponential profile convolved with the simulated PSF (blue) with
with a 300 mas box for the extended source. The solid lines indicate theresults with the single time-step reconstructor, and the dotted lines are theEE profiles of the multi time-step reconstruction with ∆t=50 ms.

For the wind speed error e∥ , the multi time-step reconstructor aﬀected more by
negative e∥ , which is an error in opposite direction to the wind, than positivee∥ . This is because e∥ is negative, the covered areas are shorted than if e∥ ispositive. The allowable range for e∥ is −5.5 m s−1 ≤ e∥ ≤ 8 m s−1 for the wide
asterism and −5 m s−1 ≤ e⊥ ≤ 7 m s−1 for the narrow asterism. The dependenceof the multi time-step reconstructor on the wind errors depends on the turbulence
power, wind speed, direction and ∆t for the multi time-step reconstructor. If the
turbulence powers are stronger, the influence from the uncertainty of wind speed
and direction is larger. Using larger ∆t also makes the tomographic error due tothe wind uncertainty larger because the spatial displacement error gets larger.

4.2.4

Detected Number of High Redshift LBGs

We evaluate the gain in multiplicity by the multi time-step reconstruction withrespect to the detected number of high redshift LBGs. The parameters used in
the evaluation are same as Table 1.1, except for the AO performance. As the
AO performance for the point source, we use the EE values measured from PSF
images simulated in the numerical simulation. The extended source is assumed
to be an exponential intensity profile with an eﬀective radius of 0.1 arcseconds.
The EE value of the extended source is estimated from the the exponential profile
convolved with the simulated PSF image.

Chapter 4. Numerical Simulation

61

5 arcminutes 10 arcminutes

5 arcminutes 10 arcminutes

single3.83 (0.83)

5.99 (2.22)

single1.26 (1.12)

4.29 (4.31)

multi

6.40 (2.96)

11.60 (8.06)

multi

1.35 (1.26)

4.63 (4.78)(a) z ∼ 5, Point source, 10 hours(b) z ∼ 5, Extended source, 20 hours

5 arcminutes 10 arcminutes

5 arcminutes 10 arcminutes

single1.71 (0.32)

2.59 (0.84)

single0.13 (0.11)

0.41 (0.41)

multi

3.00 (1.30)

5.22 (3.43)

multi

0.14 (0.13)

0.46 (0.48)(c) z ∼ 6, Point source, 30 hours(d) z ∼ 6, Extended source, 30 hours

Tab. 4.3: Detected number of high redshift LBGs at z ∼ 5 and 6. These are computedwith the observed luminosity function from [7], the simulated EE profile shownin Fig. 4.7, and the parameters used in 1.4.4. The values in and out of theparenthesis indicates the detected number with the narrow and wide asterism,
respectively.

The simulated EE values at J-band as a function of angular separation fromthe center are plotted in Fig. 4.7. The time diﬀerence for the multi time-stepreconstruction is ∆t = 50 ms. In the figure, the EE values for the point source,which is defined with a 50 mas box and represented by the red plots in Fig.
4.7, have a similar trend to the SR values at H-band shown in Fig. 4.3, but theabsolute values are diﬀerent due to the diﬀerence in the simulated wavelength.
In the case of the extended object, the EE profiles, defined with a 300 mas, box
are flatter than the EE profiles for the point source, and the diﬀerence between
the narrow and wide asterism is very small. This is because the eﬀective radiusof 0.1 arcsecond is larger than the FWHM of simulated PSF, and as a result,the EE values are dominated by the original intensity profile rather than the
convolved PSF. Therefore, the diﬀerence between the single and multi time-stepreconstruction is also small for the extended object.
Using the EE values shown in Fig. 4.7 and the observed luminosity function
from [7], we compute the detected number of LBGs at z ∼ 5, 6, which are summarized in Table 4.3. The values in and out of the parenthesis are the detected
number with the narrow and wide asterism, respectively. If the high redshift LBGs
are point sources, the detected number of LBGs gets 3 times larger by the multitime-step reconstruction with the FoR of 10 arcminute diameter compared withthe single time-step reconstruction with the 5 arcminutes FoR for both redshift
ranges. In the case of the extended object, as expected from the EE profiles shownin Fig. 4.7, there is almost no diﬀerence in the detected number of LBGs between

Chapter 4. Numerical Simulation

62

5 arcminute10 arcminuteSingleMulti

SingleMultis45,24890,49645,24890,496

φ29,76859,536

29,76859,536

ϕ40,616

41,77760,672

62,0591.35×109

5.39×109

1.35×109

5.39×109

180,992

361,984

180,992

361,984

0.013 %

0.007 %

0.013 %

0.007 %

9

9

2.49×109

1.81×10

3.70×109

658,800

1,318,960

658,628

1,318,788

0.055 %

0.053 %

0.036 %

0.036 %

99

1.75×109

3.68×10

3.85×109

200,999

206,774

300,911

307,8160.012 %

0.012 %

0.008 %

0.008 %

Γα1.21×10Pα1.65×10
L

Tab. 4.4: The three rows from the top show the number of WFS measurements, phase
points of the atmospheric turbulence layers, and the phase points on the aperture plane. The assumed parameters are listed in Table 4.1. Other rows indicate the number of total elements, non-zero elements, and the number ratio of
the non-zero elements to the total elements of the matrices of Γα , Pα , and L.

the single and multi time-step reconstruction because of the extent of object itself. The size of the FoR is more important than the EE values for the extend
case. Furthermore, the multi-object spectroscopy of the extend LBGs at z ∼ 6 isdiﬃcult even if we use the MOAO system and the future ELT.

4.2.5

Computation Eﬃciency with GPGPU

Until the previous section, we assume that the computational complexity of thetomographic wave-front reconstruction is negligible. In reality, however, it is thecritical issue in the ELT-scale tomographic reconstruction. At end of this chapter,
we present our result regarding the computation eﬃciency with GPGPU for the
ELT-scale tomographic reconstruction.
As mentioned in Section 2, the tomographic wave-front reconstruction is solved
iteratively by the CGM with the diagonal-scaling matrix. We use the GPU computation for this calculation with a graphic board of Tesla K40. The code is writtenwith C++ for CPU computation and CUDA with a version of 6.0 for the GPU
computation. The CUDA is the C-like programming code for GPU developed

Chapter 4. Numerical Simulation

63

by the NVIDIA. We use GPGPU packages of CUSPARSE for the sparse matrix
computations and CUBLAS for the vector computations.
The basic parameters and the turbulence profile used for this analysis are thesame as that shown in Table 4.1 and Table 4.2. We used the narrow asterism witha FoR of 5 arcminute diameter and the wide asterism with a FoR of 10 arcminutediameter, which are shown in Fig. 4.1. The numbers of all elements of the vectors
and the matrices are summarized in Table 4.4. For the matrices, the number andthe sparse rate, which is the number rate of non-zero elements, are also listed.
The all matrices have the sparse rate less than 1 %. We assume that there are 20
science targets.
The process of the CGM is summarized in Table 4.5 and Table 4.6. First, themeasured slopes are transferred to the memory of the GPU from the memory of
the CPU because the GPU needs to compute in own memory space. The datatransfer between CPU–GPU relatively takes time compared with the data transfers
between CPU–CPU or GPU–GPU.
The matrix A is a coeﬃcient matrix as
A = M (PαT ΓTα Γα Pα + ση2 LT L)M T ,(4.3)

where a matrix M is a diagonal scaling matrix. In order to keep the A to bea symmetric matrix, A is multiplied by M from the both the left and rightsides. In sparse matrix-vector multiplication, the computation time depends on
the number of the non-zero elements in a matrix. If Eq. (4.3) is pre-computed
before the CGM, the number of the non-zero elements gets larger than the sum of
the non-zero elements of the matrices in Eq. (4.3). As a result, using pre-computed
A for a matrix-vector multiplication increases the computation time comparedwith multiplying a vector by the matrices in Eq. (4.3) in each a matrix-vector
multiplication. Therefore, we do not pre-compute A. Since all required matrices
are transferred to the GPU oﬀ-line and kept in the GPU during the control, thetransfer time for the matrices are not included in this analysis.
The convergence is checked at the process 7 and 17 in Table 4.5 and Table 4.6.
The criteria for the convergence is determined by the relative error with respect
to the norm of the right-hand-side vector in Eq. (2.22) and the coeﬃcient a (i.e, if
ci < ac0 , we consider that the CGM convergences). We set a to 10−3 empirically.
If a is smaller values than 10−3 , there is almost no improvement in the residual
WFE and the required number of iteration increases. If a is larger than 10−3 , thenumber of iteration decreases but the residual WFE increases.

Chapter 4. Numerical Simulation

64

1.

Transfer sα (t)To GPU from CPU2.

b = PαT GTα sα (t)the right-hand-side vector

3.

ccriteria = a|b|compute criteria for the convergencea is a relative convergence criteria

4.

r0 = b − Ax0

initial residual5.

p0 = r0

initial search direction6.

c0 = |r0 |7.

Check convergence8.

f or i = 1 . . . Nmax

Loop
12–19 in Table 4.5

loop
end
9.

ϕ̂(t) = M T xest

xest is an estimates computed in the loop of step 8.

10.

φˆβ (t) = Pβ ϕ̂(t)Compute the wavefront in the science direction β

11. Transfer φˆβ (t)To CPU from GPU

Tab. 4.5: Computation flow of the CGM.

Process in the loop12.

qi−1 = Api−1

13.

γi−1 = pi−1 · pi−114.

xi = xi−1 + γpi−1

update the estimate

15.

ri = r0 + γqi−1

update residual

16.

ci = |ri |17. Check convergence18.

δ = ci /ci−119.

pi = ri + δpi−1

update search direction

Tab. 4.6: Computation flow in the loop of the CGM.

Chapter 4. Numerical Simulation

65

The computational time for each process of the CGM is summarized in Table
4.7 and Table 4.8. We estimate the computational time with the single and multitime-step tomographic reconstruction for the narrow and wide asterisms. The ∆t
for the multi time-step reconstruction is 50 ms, and the one previous measurement
is used for the multi time-step reconstruction.
Although the area of the FoR for the wide asterism is 4 times larger thanthe FoR area for the narrow asterism, the diﬀerence in the computational time
between the narrow and wide asterisms is only 60 µs in and out of the loop. Also,
the multi time-step reconstruction uses the measurements at 2 frames, but thecomputational complexity of the multi time-step reconstruction is smaller thanthe double of the computational time of the single time-step reconstructor. These
results show that the GPU computation is eﬃcient for the large problem.
The computationally heaviest part is the matrix-vector computation in theprocess 2, 4, and 12. The process 12 in the loop takes 459–699 µs, depending on
the case, and accounts for 80–85 % of the total computation time in one loop.
The vector-vector operations in process 3, 6, 13–16, and 19 take a few tens µs,
and these are much faster than the matrix-vector computations. If we assume
that the wavefront in all 20 science directions are computed by one computer,the computational time at the process 10, in which the wavefronts in the science
directions are computed from the reconstructed atmospheric turbulence layers,
takes 168 µs. If the wavefront in each science direction is computed parallel by adiﬀerent computer, it takes only 30 µs.
The data transfer time from CPU to GPU in the process 1 is the same forthe all cases. This is because the data transferred in each computation is only
the measurement at the current time-step, if the previous measurements are savedin the GPU memory. This data transfer takes 165 µs and is smaller than thematrix-vector computation, but accounts for around 10 % of the total time out of
the loop. The final data transfer from GPU to CPU in the process 11 takes 211 µs
for the system with one computer and 28 µs for the case of multiple computers.
The total computation times out of the loop are 1409–1760 µs for the caseof one computer, and the total times in the loop are 579–831 µs, depending on
asterism and reconstructor. The total computation time including the iteration is
96–162 ms, even if the diagonal scaling is used, and these are much larger than
the required computational time of 1 ms. It is expected that using the warm
start, in which the estimate at one previous time-step are used as the initial of
the CGM at the current time-step, can reduce the number of iteration largely.
However, if the number of iteration is one, the computational time of 1 ms canChapter 4. Numerical Simulation66

not be achieved with the current computational time. In order to achieve therequired computational time, accelerating the matrix-vector multiplication anddata transfer between CPU and GPU are necessary, that is, the GPU itself shouldbe improved.
In this section, we assume that the required computational time is 1 ms to minimize the temporal-lag error. However, the required computational time should
be determined with considering not only the temporal-lag error, but also a guidestar brightness and a WFS measurement noise. If the guide stars are faint and theeﬀect of the measurement noise is larger than the temporal-lag error, the performance of the MOAO is improved by controlling the MOAO system more slowly.A
simulation comparing the eﬀect of the temporal lag-error and the measurement
noise is required to determined the computational time to be achieved.

Chapter 4. Numerical Simulation

67

Narrow
process

single

Wide
multi

1.

single

multi

165

2.

249

343

251

358

3.

57

67

61

57

4.

484

676

508

724

5.

9

9

8

9

6.

36

37

39

39

7.

Check convergence

8.

loop

9.

30

30

30

30

10.

168(37)167(36)168(37)167(35)11.

211 (28)

Total time in one loop

579

773

620

831

Total time out of the loop1409(1095)

1705(1391)

1441(1121)

1760(1445)

257(165)224(165)302(203)227(194)

number of
required iteration

Tab. 4.7: Computational time and the number of iteration for the convergence of theCGM. The unit for the computational time is µs. While the value out of theparenthesis at the processes of 10 and 11 are assumed to compute the wavefront
for 20 science directions by one computer, the values in the parenthesis at theseprocesses are assumed to compute the wavefront for 20 science directions by 20
computers, i.e, it is the computeion time for one science target direction. The
tital computetinal times follow the same manner as the values at the processesof 10 and 11. The number of iteration is computed with and without the
diagonal scaling preconditioner, where the values in the parenthesis are withthe diagonal scaling preconditioner.

Chapter 4. Numerical Simulation

68

Narrow

Wideprocess

singlemulti single

multi

12.

459

652

486

699

13.

13

39

39

40

29

29

32

34

34

34

35

35

14.
15.
16.
17.

Check convergence

18.

0

0

0

0

19.

19

19

20

22

Tab. 4.8: Computational time of the processes in the loop of the CGM

Chapter 5Laboratory Test with RAVEN
5.1

RAVEN

RAVEN is a MOAO demonstrator installed and tested on the Subaru telescopeat Maunakea [32], which is the first MOAO demonstrator at 8 m-class telescope.
The RAVEN is designed not only to demonstrate the technical readiness associatedwith MOAO systems, but also to verify scientific capability of MOAO systems.
The RAVEN can apply MOAO correction simultaneously for two science targets
within a field of 2 arcminute diameter. The tomographic reconstruction is computed using the measurements from d by 3 NGSs within a FoR of 3.5 arcminutediameter and 1 LGS attached to the Subaru telescope [25] located at the centerof the FoR.
There are 4 open-loop SH-WFSs (OL-WFS) for 3 NGSs and 1 LGS in the
RAVEN. The measurements from the OL-WFSs are used for the tomographic
wave-front reconstruction, SLODAR, the wind estimation, and the tracking of
NGSs. There are two science channels each containing a DM to provide AO
correction into its science target. Each science channel also has SH-WFS behind
the DM and it is called closed-loop WFS (CL-WFS). This CL-WFS is used forthe calibration of the system, measuring a raw or residual wave-front in science
direction, providing SCAO mode correction, and tracking the science targets. All
SH-WFSs consist of 10×10 subapertures and an Andor iXon 860 EMCCD. Each
subaperture has a diameter of 0.8 m on the aperture plane and a 4.8 arcseconds
FoV. The DMs are an ALPAO DM with 11×11 actuators. The light rays comingfrom NGSs and science targets are picked oﬀ by pick-oﬀ arms and fed to each
OL-WFS or science channel. The light fed to science path is corrected by a DM.
Then, the light divided into optical light and near infrared light by a beam splitter.
The optical light is fed into the CL-WFS and the infrared light is sent into the
Subaru Infrared Camera and Spectrograph (IRCS)[29]. The RAVEN also has a69

Chapter 5. Laboratory Test with RAVEN

70

Parameters

ValuesDiameter of aperture

8m

Number of NGSs

3

Number of LGS

1

FoR for NGS

3.5 arcminutes diameter

FoR for science target

2 arcminutes diameter

WFS

SH-WFS

Number of subapertures

10×10

Number of DM elements

11×11

Number of science targets

2

Maximum Frame rate

250 HzTab. 5.1: Parameters of RAVEN.

∫

wind speed [m s−1 ]CN2 dhFraction of[km][10−14 m1/3 ]

CN2

x

y

10.5

5.961

0.180

17.0

0.0

5.5

7.418

0.224

6.0

0.0

0

19.737

0.596

0.0

5.68

Total r0 [m]0.156

AltitudeTab. 5.2: Atmospheric turbulence profile and wind profile generated by the RAVEN CU.

Figure. 5.1: Block diagram of the RAVEN system. Figure adopted from Lardière et al.
(2014) [32].

Chapter 5. Laboratory Test with RAVEN

71

1

2 fraction
CN

0.8

0.6

0.4

0.2

0
0

2

4

6

8

10

12

Altitude [km]2 fraction generated by the RAVEN CU.
Figure. 5.2: Expected CN

near infrared camera in the optical bench for the laboratory experiment without
the IRCS. These parameters are listed in Table 5.1, and Fig. 5.1 show the block
diagram of the RAVEN system.
RAVEN has a calibration unit (CU) on the optical bench, which simulates multiple GSs, multiple science targets, atmospheric turbulence layers, and a telescope[33]. The CU is used to calibrate and test the AO system. The atmospheric turbulences generated by the CU consist of three layers. Two high layers at altitude of
10 km and 5 km are simulated by phase screens, which can be rotated to simulates
a movement of the turbulence layers due to the wind. The ground layer at 0 km
is generated by a calibration DM with 17×17 actuators. The designed parametersin the CU are summarized in Table 5.2 and Fig. 5.2. It should be noted that thephase distortion patterns generated by the CU are periodic, and the period is 40 s
for the 5 km layer and 19 s for the 10 km layer. Also, the atmospheric turbulence
generated by the CU is 100 % frozen flow layer.

5.2

SLODAR and Wind Estimation

First, we test the SLODAR and the wind estimation method with the CU. An
asterism is shown in Fig. 5.3. The NGSs are at approximately 60 arcsecond radius
and the LGS is located at the center of the FoR. The brightness of the NGS is set
to R∼8 mag, which is very bright and the WFS measurement error is minimal.
Fig. 5.4 and Table 5.3 show a turbulence profile estimated by the SLODAR.
We use the slope for one minute from three OL-WFSs for the SLODAR andthe wind estimation. The estimated profile has three major peaks at the almostChapter 5. Laboratory Test with RAVEN

72

80

90"

NGS3

60
40

Y [arcsec]

★

60"
30"

20

★

NGS1

0

Ch2

✩
★
★ LGS

Ch1

-20
-40

★

NGS2

-60
-80
-80

-60

-40

-20

0

20

40

60

80

X [arcsec]

Figure. 5.3: Asterism of the GSs and the science targets for the laboratory test with the
RAVEN.

∫

wind speed [m s−1 ]CN2 dhFraction of[km][10−14 m1/3 ]

CN2

x

y

10

4.329

0.169

15.4

-0.6

8

0.896

0.035

8.6

-0.2

6

8.222

0.321

6.3

0.1

0

12.140

0.474

0.1

5.0

Total r0 [m]0.182

AltitudeTab. 5.3: Atmospheric turbulence profile and wind profile estimated by the SLODAR
and the wind estimation method.

Chapter 5. Laboratory Test with RAVEN

73

0.8
0.7

r0 = 0.182m

Fraction of CN2

0.6
0.5
0.4
0.3
0.2
0.1
0
0

2

4

6

8

10

12

Altitude [km]Figure. 5.4: Profile of atmospheric turbulence of the RAVEN CU estimated by the SLODAR.

same altitudes as the designed profile in Fig. 5.2. The smallest turbulence at
8 km in the estimated profile may come from the turbulence layers at 5.5 km
or/and 10.5 km due to the mismatch of the altitude between the SLODAR andthe designed turbulence layers. The largest diﬀerence between the designed andmeasured profiles is the CN2 value at ground layer, and the estimated value of
12.140×10−14 m1/3 is 60 % of the designed value. As a result, the total power
measured by the SLODAR is smaller than the designed value, i.e., the measured
r0 of 0.182 m is larger than the designed r0 of 0.156 m. The Cn2 of the other highlayers are consistent between the measured and the designed values.
This larger r0 than the designed value is already reported in the previous
estimation in the lab on RAVEN, which is estimated by a diﬀerent SLODAR
method from the one presented in this thesis. The diﬀerence between the two
method is that, while the method presented in this thesis uses the correlation
map of the slope form multiple SH-WFSs, which is defined in Eq. (2.26), the othermethod uses the covariance matrix of the slopes. The covariance matrix is definedas⟨⟩∑′′x[1]′′x[2]
obs(i, j, i′ , j ′ ) =Covxx
si,j (t)si′ ,j ′ (t) ,(5.1)
valid i,j,i′ ,j ′where s′′ is a slope after the tilt and average pattern subtraction defined in Eq.
(2.25). Both methods fit the theoretical correlation or covariance matrix to themeasured one to estimate CN2 (h). We refer to the SLODAR presented in thisthesis as a correlaton SLODAR. The other SLODAR method is referred to as acovariance SLODAR.
The phase distortion at ground layer in the CU is generated by the calibration

Chapter 5. Laboratory Test with RAVEN

74

60ms

120ms

180ms

0ms

60ms

120ms

180ms

0ms

120ms

240ms

360ms

0ms

120ms

240ms

360ms

0km

6km

8km

10km

0ms

Figure. 5.5: Temporal correlation maps of slope at each altitude estimated by tomographic reconstruction. The values at the top of each image indicate the timedelay used to drive the temporal correlation. The x and y-axes are δi andδj with a range of −9 ≤ δi, δj ≤ 9. One pixel on the temporal-correlation
corresponds to one subaperture size, which is 0.8 m in the RAVEN

Chapter 5. Laboratory Test with RAVEN

75

DM. Since the actuator size of the calibration DM is roughly 0.47 m, the DM
can not reproduce the phase distortion with spatial scales smaller than 0.47 m.
We conclude that the actuator size causes the discrepancy and that the estimatedvalue reflects the real turbulence profile on the DM better than the designed value.
By using the estimated profile, we evaluate the wind speed and direction ateach layer. First, we compute the tomographic reconstruction to reconstruct thephase distortion pattern at each altitude. Then, the temporal correlation map of
the estimated slope are calculated with diﬀerent delay times, δt. The temporalcorrelation maps are illustrated in Fig. 5.5. The delay time for the lower layers at
0 km and 6 km are set to larger than the higher layers (Fig. 5.5). Generally, thewind speed at low layer is considered slower than the wind speed at high altitude.
In order to get clear movement of the temporal correlation peak, longer delay timeis preferred for the low altitude layers. For all layers detected by the SLODAR,
the peaks on the temporal correlation maps are detected. The estimated windspeeds and directions, which are summarized in Table 5.3, are consistent with thedesigned value. The error of the wind speed is less than 2 m s−1 for all layers, which
is an acceptable wind error computed by the numerical simulation in Section 3.
With this wind estimation error, the multi time-step reconstructor will work withthe future MOAO system (see Fig. 4.6).

5.3

Multi Time-Step Reconstruction

We present the result of the laboratory test of the multi time-step reconstructionin Fig. 5.6 and the top table of Table 5.4. Fig. 5.6 shows the PSF images of two
science channels taken by RAVEN infrared camera with the MOAO correction.
The pixel scale of the PSF images in Fig. 5.6 is 17.5 mas, and the size of the images
are 0.5 arcsecond × 0.5 arcsecond. The maximum wind speed simulated by the
CU , which is 17.0 m s−1 , is almost half of the maximum wind speed in the modelthat we used for the numerical simulation in Section 3, and the time durationthat the decay ratio fdecay of the temporal correlation of SH-WFS measurements
is more than 0.7 with the wind speed of 17.0 m s−1 is ∼100 ms. Therefore, we
use 100 ms for ∆t of the multi time-step reconstructor. The asterism used for
the MOAO experiment in the laboratory is the same as the asterism used for the
SLODAR, which is shown in Fig. 5.3.
The performance of the MOAO is evaluated by SR, EE, and improvement
ratios of the SR and the EE at wavelength of 1650 nm (H-band). First, weChapter 5. Laboratory Test with RAVEN

76

perform a background subtraction to the original image observed by the sciencecamera in RAVEN. Then, the PSF image of each science channel is cropped with 1
arcsecond × 1 arcsecond box from the original image and normalized by the total
intensity within the cropped box. In parallel, we create a diﬀraction-limited PSF
image by a simulation with a pixel scale of 1 mas and normalize the diﬀractionlimited PSF image within a 1 arcminute box. We estimate a SR by comparing
the peak intensity in the observed PSF image with that in the diﬀraction-limited
PSF image. Accurately, this SR is an EE within a 1 mas box , but the 1 mas box
is small enough compared with the FWHM of the PSFs and it can be considered
as a SR. The EE is computed within a 140 mas box, which is a size of the slit of
the slope. The measured values are summarized in the top table of Table 5.4
Visually, the PSFs become sharper and have higher peak intensity by the multitime-step reconstruction compared to the single time-step reconstruction for bothof the asterisms, as shown in Fig. 5.6. The EEs achieved by the single time-stepreconstructor are 0.268–0290 without the LGS and 0.349–0.365 with the LGS.
The improvement ratio of the EE achieved by the multi time-step reconstructionis 1.16–1.18 without LGS and 1.05 with LGS. The SR of the single time-stepreconstructor without the LGS is 0.072–0.087 and these values are roughly consistent with the value of 0.10 measured previously in the laboratory experiment on
RAVEN [32]. On the other hand, the SR of the single case with the LGS is 0.125–
0.150, which is much lower than the previous value of 0.23. The improvement ratiosof the SR achieved by the multi-time step reconstructor are 1.32–1.51 without the
LGS and 1.22–1.25 with the LGS. These improvements show the multi time-stepreconstruction works in the lab-test. The improvements in the SR and the EE are
better without the LGS than with the LGS. In this laboratory experiment, since
the part of the tomographic error due to the uncovered and unoverlapped area issolved with the LGS and the tomographic error would be dominated by the othercomponents (e.g. the open loop error, the temporal lag-error, the optical error,and so on) rather than the uncovered and unoverlapped area. Therefore, the improvement by the multi time-step reconstruction becomes smaller relatively thanthe case without the LGS.

5.4

Comparison with a numerical simulation

We compare the results of the lab-test on the multi time-step reconstructor toa numerical simulation. We assume the parameters of the RAVEN system andthe turbulence profile generated by the RAVEN CU, listed in Table 5.1 and TableChapter 5. Laboratory Test with RAVEN

Lab-test

EE (140mas)

reconstructor
single
multi

77

GSs
3NGSs

Improvement ratio

SR

Ch1

Ch2

Ch1

Ch2

0.268

0.290

0.072

0.087

0.313

0.343

0.095

0.131

1.16

1.18

1.32

1.51

single

3NGSs

0.349

0.365

0.125

0.150

multi

+LGS

0.367

0.385

0.153

0.187

Improvement ratio

1.05

1.05

1.22

1.25

Simulation
reconstructor
singlemulti

EE (140mas)
GSs

3NGSs

Improvement ratio

Ch1

Ch2

SR
Ch1

Ch2

0.334(0.294) 0.346(0.303)

0.251(0.211)

0.267(0.224)

0.406(0.348) 0.424(0.362)

0.361(0.302)

0.384(0.321)

1.21(1.18)

1.23(1.19)

1.44(1.43)

1.44(1.43)

single3NGSs

0.443(0.379)

0.430(0.367)

0.398(0.333)

0.386(0.324)

multi

+LGS

0.477(0.402)

0.477(0.402)

0.465(0.390)

0.465(0.389)

Improvement ratio

1.07(1.06)

1.11(1.10)

1.17(1.17)

1.20(1.20)

Tab. 5.4: SR, kSR , and EE value of each science channel measured in the lab-test (thetop table) compared to those predicted from the numerical simulation (thebottom table). The EE is defined with a 140 max box. There are two GS
configuration used in the test, one is only three NGSs and the other is three
NGSs+LGS. The values in parenthesis in the lower part is the value accounting
for the implementation errors.

Chapter 5. Laboratory Test with RAVEN
Single Time-Step (3NGSs)

Ch1 : EE=0.268
SR=0.072

Ch2 : EE=0.290
SR=0.087

78
Multi Time-Step (3NGSs)

Ch1 : EE=0.313
SR=0.095

Ch2 : EE=0.343
SR=0.131

Single Time-Step (3NGSs+LGS)Multi Time-Step (3NGSs+LGS)

Ch1 : EE=0.349
SR=0.125

Ch1 : EE=0.367
SR=0.153

Ch2 : EE=0.365
SR=0.150

Ch2 : EE=0.385
SR=0.187

Figure. 5.6: PSF images of each science channel taken in the lab-test with the single andthe multi time-step reconstructors. The color scale is linear and aligned foreach channel. The wavelength is 1650 nm (H-band). The size of each image
is 0.5 arcsecond × 0.5 arcsecond, and the pixel scale is 17.5 mas.

5.2. The coordinates of the GSs and the targets are also set to the same as thecoordinates used in the lab-test, illustrated in Fig. 5.3. The EE values and theimprovement ratio, kSR , predicted by the simulation are summarized in the lower
part of Table 5.4. There is a good agreement between the kSR measured in thelab-test and derived by the simulation, therefore the performance improvementof the multi time-step reconstructor is as good as predicted for both asterisms
with and without LGS. It is also consistent between the laboratory test and thesimulation that the improvement by the multi time-step reconstruction without
LGS is larger than the improvement without the LGS.
However, the absolute values of the EE and the SR measured in the laboratory experiment are much smaller than the values predicted by the simulation forboth of the single and the multi time-step reconstructors. This diﬀerence may
be due to the implementation errors not considered in the simulation, which arecalibration error, DM control error, and optical aberration. In the calibration of
RAVEN, there are two steps [28]. The first one is to find the rotation matrices,which compensate the misalignment between the OL-WFSs and the CL-WFSs.
In particular, in RAVEN, OL-WFSs are rotated respect to the CL-WFSs. The
second one is to determine the command matrices to control the science DMs.
The rotation and command matrices are determined experimentally using the CU
and the CL- and OL-WFSs. Therefore, these include the noise.
The calibration error results from the uncertainties of these matrices. The DM

Chapter 5. Laboratory Test with RAVEN

79

control error is regarding to the accuracy and long term stability of the DM, i.e.,
how the DM can reproduce the demanded wave-front. The optical aberration is astatic error caused by the misalignment in the system. The optical aberration in
the science path can not be seen from the OL-WFSs, and therefore, this aberration
can not be corrected by the tomography. This aberration is called non-common
path aberration (NCPA). We set to the reference shape of the DMs to remove thestatic optical aberration [31], but the residual aberration may aﬀect the MOAO
performance. Anderson et al. (2012) [4] simulates the performance of the RAVEN
and estimates the error budget of the AO correction of the RAVEN system. They
estimate a wavefront error (WFE) caused by the implementation errors of the
RAVEN to be 107 nm.
In order to evaluate the contribution from the implementation errors, we calculate the total WFE as a quadrature sum of WFE due to the implementation andtomographic errors. The WFE due to the tomographic error is estimated from theSR derived from the simulation by using an approximation, SR=exp [−σ 2 ], whereσ = 2πWFE/λ and λ=1650 nm. It is noted that this expression is valid for σ 2 < 1
rad2 (this is SR< 0.37), and our result can be biased when SR is lower than 0.37.
The SR, including the tomographic error and the implementation errors, can becomputed with this approximation from the total WFE. Since there is no approximation connecting to EE and WFE, we assume that the absolute loss of the EE
due to the implementation errors is same as the absolute loss of the SR. The performance predicted by the simulation, accounting for the implementation errors,
are listed in the bottom table of Table 5.4 as the values in a parenthesis. While thesimulated EE including the implementation errors is slightly higher than the EE
values measured in the laboratory experiment, the simulated SR including the implementation errors is still much higher than the measured SR. This suggests that
there are additional errors, which make a PSF peak blurred and aﬀects mostly SR
values. The improvement ratios does not vary with the implementation errors.
As the additional errors, we examine the eﬀect of the residual tip-tilt components. Since the tip-tilt errors change the position of the PSF, the fluctuation of
the tip-tilt errors aﬀects strongly the SR of the long-exposure PSF. On the other
hand, the EE is expected to be less sensitive to the tip-tile error than the SR because the EE is evaluated within a 140 mas box which is larger than the FWHM
of the corrected PSF. In the simulation, the standard deviations of the residual
tip-tilt errors are roughly 6.5 mas for both of the x and y-direction. The residual
tip-tilt error in the laboratory experiment can be estimated from the CL-WFS
measurements, and the standard deviations of the residual tip-tilt errors is 15 mas

Chapter 5. Laboratory Test with RAVEN

80

for x-direction and 20 mas for y-direction in the laboratory experiment. These
are much larger than that of the simulation and suggest that the results in thelaboratory experiment are aﬀected more by the tip-tilt errors than the result of thesimulation. The sparse approximation of the inverse of the phase covariance matrix, L, may contribute the uncorrected additional tip-tilt errors in the laboratory
test because the sparse approximation can not regularize the tip-tilt modes [35].
The standard deviations of the tip-tilt modes of the uncorrected wave-front are 55
mas in the simulation and 50 mas in the laboratory experiment. Therefore, theadditional vibration from the RAVEN does not aﬀect the result in the laboratory
experiment.
In order to evaluate the contribution from the tip-tile errors, we add the tip-tilt
errors to the corrected wave-front in the simulation so that the standard deviationof the tip-tilt errors in the simulation becomes equal to that of the laboratory
experiment. The re-simulated values with the additional tip-tilt errors are listed
in Table 5.5. The values in a parenthesis in Table 5.5 is the SR including the eﬀectof the implementation errors. As expected, while the SR is aﬀected strongly on the
additional tip-tile errors, there is almost no change in the EE due to the additional
tip-tile errors. The simulated SRs including the eﬀect of the uncompensated tiptilt and the implementation errors are much closer to the measured values than theoriginal simulated SRs in Table 5.4, but these are still higher than the measured
values.
In this section, we compare the result in the laboratory with the numericalsimulation and discuss the eﬀect of the implementation errors and residual tip-tile
errors to explain the lower SRs and EEs in the laboratory experiment than that of
the simulation. With considering these eﬀects, the diﬀerence between the result of
the laboratory experiment and the simulation gets smaller, but the SRs and EEsof the laboratory test are still lower than the values predicted by the simulation.
In order to understand what makes this diﬀerence, more experiments to examine
the implementation errors in detail are required.

Chapter 5. Laboratory Test with RAVEN

Simulation
reconstructor
singlemulti

81

EE (140mas)
GSs

3NGSs

Improvement ratio

SR

Ch1

Ch2

Ch1

Ch2

0.326

0.337

0.137(0.116)

0.149(0.126)

0.397

0.413

0.194(0.164)

0.204(0.173)

1.22

1.23

1.42

1.37

single

3NGSs

0.432

0.418

0.213(0.180)

0.211(0.179)

multi

+LGS

0.466

0.464

0.250(0.212)

0.247(0.209)

Improvement ratio

1.08

1.11

1.17

1.17

Tab. 5.5: SR, kSR , and EE value of each channel predicted from the numerical simulation with taking account of the additional tip-tilt errors. The additional
tip-tilt errors are determined by the measurement on the laboratory experiments. The values in a parenthesis in Table 5.5 is the SR including the eﬀectof the implementation errors.

Chapter 6
On-Sky Evaluation with RAVEN
6.1

On-Sky Performance

The on-sky engineering and science observations with RAVEN on the Subaru
telescope are successfully performed on May and August 2014 and June 2015. We
test the multi time-step reconstruction and the wind estimation method on sky
during the third observation on June 2015. In this chapter, we present the onsky performance of the SLODAR, the wind estimation, and the multi time-stepreconstruction based on the data from the third on-sky observation.

6.1.1

Engineering Field

The engineering fields of the third run are shown in Fig. 6.1. The green points
in the figure indicate the position of the NGSs and its magnitude at R-band arelisted in Table 6.1. All NGSs are brighter than R=12.6. The limit magnitude of
OL-WFSs in RAVEN is around 14 mag at R-band. Target objects are selected inradius of 1 arcminute from the center.

6.1.2

SLODAR

The performance of the SLODAR of RAVEN in the first and second on-sky observations is reported in [32]. The result of the SLODAR in [32] is computed by thecovariance SLODAR method. The covariance and correlation SLODAR are implemented in RAVEN and provide the almost same results, but the computation
time of the covariance SLODAR is a few seconds faster than the correlation SLODAR, where the covariance SLODAR takes around 5 seconds for its computationand the correlation SLODAR takes around 10 seconds. Therefore, the covariance
SLODAR is usually used during on-sky observations.

82

Chapter 6. On-Sky Evaluation with RAVEN

Fd38

83

Fd327
1	

2	

1	

2	

3	

3	

(a) Fd38(b) Fd327

Fd406

Fy1
1	

1	

2	
3	

2	

(c) Fd406

3	(d) Fy1

Figure. 6.1: Images of the engineering fields. The red lines indicate the radius of 1.75
arcminute and 1 arcminute from the center of the field. The green points
represent the NGSs position. The brightnesses of the NGSs are summarizedin Table 6.1. The information of the engineering field is summarized in Table
6.1. These images comes from the Digital Sky Survey. The position of NGSs
are plotted based on the USNO-A2.0 catalog.

Chapter 6. On-Sky Evaluation with RAVEN

84

Field IDRA

Dec

Field ID

RA

Dec

Fd38

15:15:30+5:15:57

Fd38

19:45:00+33:14:41

NGS

Rmag

NGS

Rmag

1

12.2

1

11.0

2

12.5

2

10.6

3

10.2

3

11.0

(a) Fd38(b) Fd327

Field IDRA

Dec

Field ID

RA

Dec

Fd38

20:11:28+35:51:58

Fd38

19:45:26+33:25:32

NGS

Rmag

NGS

Rmag

1

11.4

1

10.9

2

12.6

2

10.6

3

12.0

3

11.4

(c) Fd406(d) Fy1

Tab. 6.1: RA, Dec, and NGS magnitudes of the engineering field shown in Fig. 6.1. The
NGSs magnitude are referenced from the USNO-A2.0 catalog.

Fig. 6.2 shows an example of the CN2 profile and the cross-correlation maps
computed from one minute series of the on-sky measurements by the correlation
SLODAR. The field is Fd38. In the profile, the ground layer has more than 80% of the total turbulence power and the weak turbulences are at around 12 km.
Only one central peak is seen on the cross-correlation maps, which corresponds to
the ground layer. Fig. 6.3 shows another example of a SLODAR result with CN2
profile with a significant component at high altitudes. The field is Fd406. The
turbulences spread over a wide range of altitudes and 50 % of the total turbulence
power occurs at the altitudes higher than 10 km. On the cross-correlation maps,
there are strong correlations along the base lines of the NGS pairs.
In most of the time during the on-sky observations, we detect a strong groundlayer as expected [16]. Fig. 6.4 shows a histogram of the ground layer fraction
measured by the covariance SLODAR during the on-sly observations. The groundlayer accounts for the 50 % of the total turbulence power at 75 %-ile. This strong
ground layer suggests the large contribution from the dome seeing.

Chapter 6. On-Sky Evaluation with RAVEN

80

85

2
0.9

NGS 1	

20
0

−20
−40

NGS 3	NGS 2	

−60

0.8
1.5

Fd38, 2015/6/23 19:48:00

0.7
0.6

r0 =0.206 m

0.5

1

0.4
0.3
0.5

0.2

fraction of CN2 (h)dh

Y [arcseconds]40

CN2 (h)dh × 10−13 [m1/3]

60

0.1

−80
−80 −60 −40 −20 0
20 40
X [arcseconds]

0
60

80

WFS 1-2	

0

5

10

15

20

0

Altitude [km]WFS 2-3	

WFS 3-1	

Figure. 6.2: The example of the SLODAR result dominated with a grund layer. The
upper left panel shows the NGS asterism and the dashed circle indicates theangular separation of 60 arcseconds from the center of the FoR. The upper2 profile estimated by the correlation SLODAR andright panel shows the CN
the bottom panel shows the corresponding correlation maps. The red arrows
indicate the baseline of the guide stars, and the signals of the correlation are
expected to apper along these arrows. In this case, most of the turbulence isconcentrated in the ground layer and there are weak turbulences at around
12 km.

Chapter 6. On-Sky Evaluation with RAVEN

60

3.5

NGS 1	Y [arcseconds]

40
20

NGS 3	0

−20

NGS 2	

−40
−60
−80
−80 −60 −40 −20 0
20 40
X [arcseconds]

380

WFS 1-2	

0.35

r0 =0.086 m

0.3

2.5

0.25
2

0.2

1.5

0.15

1

0.1

0.5

0.05

0
60

0.4

Fd406, 2015/6/25 22:02:00

0

5

10

15

20

fraction of CN2 (h)dh

4CN2 (h)dh × 10−13 [m1/3]

80

86

0

Altitude [km]WFS 2-3	

WFS 3-1	

Figure. 6.3: Another example of the SLODAR result with a complex profile. The turbulences occur over wide range of altitudes and the turbulences at altitude
higher than 10 km have significant contribution to the total turbulence
power. On the correlation maps, there are strong correlations along thebase lines of the NGS pairs

Chapter 6. On-Sky Evaluation with RAVEN

6.1.3

87

Wind Estimation

Fig. 6.5 and Fig. 6.6 show the wind speeds, wind directions, and the temporalcorrelation maps with diﬀerent time delays, δt, at each altitudes for the cases ofFig. 6.2 and Fig. 6.3, respectively. In Fig. 6.5, which is the result of the profile
dominated with a ground layer, a peak on each temporal-correlation map seems
not to move from the center even with increasing the time delay, δt, up to 100
ms. This suggests that there is slow or not moving component at all altitude.
As a result, the measured wind speeds are 1 m s−1 at all altitudes. The measured
directions are diﬀerent between the ground layer and high altitude layers.
On the other hand, in Fig. 6.6, which is the case that the turbulences spreadover a wide range of altitudes, the peaks on the temporal correlation maps slowly
move to left with the increase of δt, except for the ground layer. The measured
wind speeds are almost 0 m s−1 at 0 km and around 5 m s−1 at higher altitudes.
In most of the time during the on-sky observations, we detect temporal-correlation
peak moving slowly or not moving with the time delay δt. These peaks are seenat all altitudes as shown in Fig. 6.5 and Fig. 6.6. As a result, the average of thedetected wind speeds are very small, which are equal to or less than 1 m s−1 , forall altitudes.
Although, at low altitudes, this slow wind speeds is not strange because theaverage wind speed at low altitudes are generally considered to be slow comparedwith wind speeds at high altitudes [24], there would be some contributions fromthe dome seeing. The turbulence caused by the dome seeing may be boilingrather than moving with the wind. As explained in Chapter 2, we subtract thetemporal average over a minute from the reconstructed slopes before computing
the temporal-correlations in order to remove the static pattern in one minute. If
the temporal evolution of the dome seeing is faster than 1 minute, it makes thecorrelation peak on the temporal-correlation map of low altitudes, which does not
move with δt. As a result, the estimated wind speeds get small due to the dome
seeing.
For the high altitudes, these slow wind speeds and the temporal correlation
peak without moving are not considered to reflect the real wind speeds, since theaverage wind speeds is faster than 10 m s−1 [24]. It is possible that the temporalcorrelation peaks due to moving atmospheric turbulence layers are hidden by thestrong peak associated with the additional slow or not moving component. We will
discuss this slow wind speeds and the temporal correlation peak without moving
in the oﬀ-line analysis.

120

25%ile75%ile140

8850%ileChapter 6. On-Sky Evaluation with RAVEN

Number

100
80
60
40
20
0
0

0.2

0.4

0.6

0.8

1

Ground Layer Fraction
Figure. 6.4: Histogram of ground layer fraction measured by the covariance SLODAR
during on-sky observations from June 23th to July 2nd in 2015 (HST).

6.1.4

On-sky Performance of the Multi Time-step Reconstruction

We evaluate the performance of the multi time-step reconstruction with the PSF
images taken by IRCS during the on-sky observations. All PSF images which
we use are observed by H-band. The exposure time for each image is roughly
30 seconds. The pixel scale of IRCS corresponds to 0.023 arcseconds, which ismeasured from the image of the binary star taken during the on-sky observations.
We use FWHM, SR, and EE in a 140 mas box for the performance evaluation. These metrics are aﬀected strongly on the background subtraction from theimages. Therefore, the same criteria to subtract the background should be used
for all images. We define the count of the background at radius of 10 times of its
FWHM. When the FWHM is larger than 0.3 arcsecond, the radius of 10 times of
its FWHM is out of the FoV because the size of the FoV of each science channel is
4 arcsecond. In this case, we define the back ground count at the radius as large
as possible.
Fig. 6.7 shows the FWHM, SR, and EE measured from the on-sky PSF images. These metrics are plotted as a function of a r0 measured by the covariance
SLODAR. The red and blue points indicate the results of the single and multitime-step reconstruction, respectively. For the multi time-step reconstruction, we
set ∆t to a time duration corresponding to 20 frames. In most case, the frame

Chapter 6. On-Sky Evaluation with RAVEN

89

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0km

9km

10km

12km

14km

0ms

30
wind speeds

150

25Wind speed [m/s]

100
20

50

15

0

10

−50

wind direction [degree]

wind directions

−100

5

−150
0
0

5

10

15

20

25

Altitude [km]Figure. 6.5: Top: Temporal-correlation maps with diﬀerent time delays, δt, at each altitudes for the cases of Fig. 6.2. The x and y-axes are δi and δj with a range of−9 ≤ δi, δj ≤ 9. One pixel on the temporal-correlation corresponds to one
subaperture size, which is 0.8 m in the RAVEN. The temporal-correlation
maps are normalized by the peak intensity of the map of δt=0 at each altitude. The color scale of the maps is set to [0 0.7]. Bottom: Wind speed(red) and direction profile (blue) as a function of altitude from the telescope.
The wind speed follows the left y-axis, and the wind direction follows theright y-axis.

Chapter 6. On-Sky Evaluation with RAVEN

90

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0km

2km

4km

6km

8km

0ms

30
wind speeds

150

25Wind speed [m/s]

100
20

50

15

0

10

−50

wind direction [degree]

wind directions

−100

5

−150
0
0

5

10

15

20

25

Altitude [km]2 profile shown in Fig. 6.3.
Figure. 6.6: The same figure as Fig. 6.5 with the CN

Chapter 6. On-Sky Evaluation with RAVEN

91

1
FWHM at H-band [arcsecond]Single, Ch1
Single, Ch2

0.8

Multi, Ch1
Multi, Ch2

0.6

0.4

0.2

0
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.3

0.35

0.4

0.3

0.35

0.4

r0 at λ = 500nm [m]Strehl Ratio(a) FWHM
0.14

Single, Ch1

0.12

Single, Ch2
Multi, Ch1

0.1

Multi, Ch2
0.08
0.06
0.04
0.02
0
0

0.05

0.1

0.15

0.2

0.25

r0 at λ = 500nm [m](b) Strehl Ratio

Ensqured Energy in a 140 mas box

0.4
Single, Ch1
0.35

Single, Ch2

0.3

Multi, Ch1

0.25

Multi, Ch2

0.2
0.15
0.1
0.05
0
0

0.05

0.1

0.15

0.2

0.25

r0 at λ = 500nm [m](c) Ensqured Enegy in 140 mas bax

Figure. 6.7: The (a) FWHM, (b) SR, and (c) EE in 140 mas box measured from the
PSF images obtained in the on-sky engineering observations. These metrics
are plotted as a function of a r0 . The red squres indicate the results of thesingle time-step reconstructor and the blue circles are the result of the multitime-step reconstructor. The filled and open symbols indeciate the result of
the science channel 1 and 2, respectively. In the enginnering observation,
we set ∆t to a duretion corresponding to 20 frames. In the most case, theframerate of the control for MOAO system is 150 Hz and ∆t is 133 ms.

Chapter 6. On-Sky Evaluation with RAVEN

92

rate of the MOAO control is 150 Hz and therefore ∆t is 133 ms. We achieve the
FWHM of 0.1 arcseconds and EE more than 0.3 when r0 =0.3 m. On the other
hand, there is the large variance of the EEs at the same r0 =0.3 m. One reason
is the variation of the NGS asterisms. Another reason is the eﬀect of the DM
instability. The initial command of the DM is set to cancel the phase distortiondue to the optical error. This best flat command is determined experimentally
before the observations and used as the reference command of the DM during theobservation. However, we found that the best flat command changes after we use
the DM continuously over a long time. Therefore, especially in the latter half ofthe night, the DMs become unstable and decrease the performance.
In Fig. 6.7, there is no clear improvement due to the multi time-step tomographic reconstruction. The error of the wind estimation may be a main cause
which degrades the performance of the multi time-step reconstruction. Fig. 6.8
shows the PSF images observed with the correction by the single and the multitime-step tomographic reconstruction alternately with two minutes intervals. The
field is Fd38. Since the PSF images of the science channel 1 are aﬀected by the
DM instability, we focus on the result of the science channel 2. The FWHMsof the single time-step tomographic reconstructor are a few mas better than the
FWHM of the multi time-step reconstruction, but this diﬀerence is much smaller
than the pixel scale of 23 mas. It is diﬃcult to conclude that this diﬀerence is dueto the diﬀerence of the reconstructor. The values of SR look to depend on thefried parameter r0 . The EEs in a 140 mas box of the multi time-step tomographic
reconstructor is 0.02 larger than that of the single time-step reconstructor. Fig. 6.9
shows the CN2 profile the wind speeds and directions used for the reconstruction at
this time. The wind speeds at all altitudes are almost 0 m/s, which is aﬀected by
the slow or not moving component which we mentioned in the previous subsection.
In such a case, the improvement with the multi time-step reconstruction can not
be expected. Therefore, the diﬀerence in the EEs may not be due to the diﬀerenceof the reconstructor.
Unfortunately, we can not achieve the clear improvement due to the multitime-step tomographic reconstruction from the on-sky PSF images during engineering observations. In the next section, we improve and evaluate the windestimation method and the multi time-step tomographic reconstruction with theoﬀ-line analysis using the on-sky measurements of SH-WFSs.

Chapter 6. On-Sky Evaluation with RAVEN

r0 [m]	

0.296	

93

0.305	

Single	0.353	

Multi	

0.285	

Single	Multi	

Ch1	

FWHM [mas]	

124	

132	

154	

242	

SR [%]	

0.053	

0.040	

0.035	

0.023	

EE [%]	

0.250	

0.198	

0.180	

0.145	

28

254

Multi	

636

1022

Single	

1404

1790

2172

2554

Multi	

2940

Single	

3322

37

Ch2	
FWHM [mas]	

106	

101	

105	

103	

SR [%]	

0.078	

0.083	

0.084	

0.074	

EE [%]	

0.332	

0.309	

0.333	

0.307	

Figure. 6.8: The on-sky PSF images corrected by the single and multi time-step re21
1292
1950
306
634
964
1623
2609
2936
32
2278
costructor. The size of PSF images is 1.5×1.5 arcsecond.

30
wind speeds (left)

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

150

wind directions (right)

25
Wind speed [m/s]

100
20

50

15

0

10

−50
−100

5

0

5

10

15

Altitude [km]2025wind direction [degree]

0.8

0.8

2 (h)dhfraction of CN

2 (h)dh × 10−13 [m1/3 ]
CN

0.9

−150
0
0

5

10

15

20

25

Altitude [km]2 values (left panel), wind speeds (right panel), and wind
Figure. 6.9: The profiles of CN
directions (right panel) measured on sky. The field is Fd38 and the asterism
2 profile, more than 80 % of the turbulence areis same as Fig. 6.2. In the CN
concentrated in the ground layer and there are weak turbulences at around
10 km. On the left panel, the detected wind speeds are almost zero at all
altitudes.

Chapter 6. On-Sky Evaluation with RAVEN

6.2

Oﬀ-Line Analysis

6.2.1

Improvement in the Wind Estimation

94

First, we improve the wind estimation method. As mentioned in subsection 6.1.3,
the wind speeds are seen to be underestimated due to the additional slow or not
moving component especially at high altitudes. As explained in Chapter 2, wecompute the temporal averaged pattern from the time series of the reconstructedslopes and subtract the computed averaged pattern from the reconstructed slopes
to remove static pattern. In the revised wind estimation method, we perform
this process for every one second time series of the reconstructed slopes before
computing temporal-correlations, and cancel the additional slow or not moving
component from the temporal-correlation maps to detect only moving turbulencelayers. This process corresponds to remove the component changing slower than 1
Hz or moving slower than ∼0.4 m s−1 from the temporal-correlation maps. For the
ground layer, the dome seeing may be removed by this process. Fig. 6.10 and Fig.
6.11 are the revised versions of Fig. 6.5 and Fig. 6.6, respectively. In both of Fig.
6.10 and Fig. 6.11, the only moving peaks are successfully extracted. The measured
wind speeds are larger than the value measured during the on-sky observation in
real-time, and thus, the underestimation is solved by the revised wind estimation
method.
Fig. 6.12 shows example of of the phase distortion pattern averaged over each
one second at diﬀerent altitudes. The averaged pattern at 0 km and 4 km change
slowly with time and some distortion patterns remain over several seconds. Theselong-lived distortion patterns cause the temporal-correlation peak not moving with
δt and the underestimation of the wind speeds. Although the averaged distortion
pattern at 8 km change more rapidly than that of 0 km and 4 km, sometimes
the adjacent distortion maps have marginally same pattern and it is possible thatthese pattern make the temporal-correlation peak not moving.
At the 0 km, the slowly changing distortion pattern can be explained by thedome seeing, but it is diﬃcult to consider that there are additional components
slowly changing at high altitudes, which are diﬀerent from the moving turbulencelayer. In Fig. 6.12, the similar phase distortion patterns are seen on the averaged phase map for all altitudes. These maps are constructed by the tomographic
reconstruction. In other words, it is possible that the reconstructed phase distortion or slope is aﬀected by the degeneracy due to the unoverlapped area, theuncertainty of the turbulence height and strength measured by the SLODAR, andChapter 6. On-Sky Evaluation with RAVEN

95

the measurement noise. Therefore, these similar patterns, which change slowly as
shown in Fig. 6.12, may be caused by the error of the tomographic reconstruction.
In the laboratory experiment, since the magnitude of the NGSs are set to around8 mag, the tomographic reconstruction is performed more accurately comparedwith the on-sky observations, and thus, the slow or not moving correlation peakis not seen in Fig. 5.5. Furthermore, the turbulence layers generated in the laboratory experiment follows perfectly the frozen flow assumption, and therefore, thetemporal-correlation peaks are clearer than the peaks of the on-sky observations.

6.2.2

Atmospheric Turbulence Profile

Here, we summarize the profiles of r0 , CN2 (h), wind speeds and wind directions in
the night on June 23–24th, 2015 measured from the on-sky measurements of the
open-loop SH-WFSs with the oﬀ-line analysis.
Fig. 6.2.2 shows the r0 profiles measured by the correlation SLODAR of RAVEN(the blue plots), the auto-correlation fitting of RAVEN (the red plots), and the
CFHT DIMM (the gray lines) as comparison. In the later half of the night, theRAVEN results are better than the result form the CFHT. This is may be due tothe local diﬀerence since the CFHT sit on a diﬀerent ridge from the Subaru telescope. There is a large diﬀerence between the result from the SLODAR and theauto-correlation fitting in 2:30–3:30. In this time, since the NGS asterism is wide,
the SLODAR can measure the turbulence up to around 9 km and the SLODAR
misses the some turbulence at altitude higher than 9 km. On the other hand, theauto-correlation fitting can measure the total turbulence power including the turbulence missed by the SLODAR. Therefore, there is the large diﬀerence between
the result from the SLODAR and the auto-correlation fitting in this period.
Fig. 6.14 is the CN2 profile. In this day, the atmospheric turbulence is dominatedby the ground layer in all the time. Also, there are weak turbulences at around
10 km. The CN2 profile of this day looks stable throughout the night.
Fig. 6.15 and Fig. 6.16 are the profile of the wind speeds and directions. The
ground layer has the wind of 5 m s−1 and the higher layers have the wind of 10–15 m
s−1 . Similar to the CN2 profile, the wind speed and condition are stable throughout
the night. This stability of the wind speeds and directions is essential for the multitime-step reconstruction.

Chapter 6. On-Sky Evaluation with RAVEN

96

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0km

9km

10km

12km

14km

0ms

30
wind speeds

150

25Wind speed [m/s]

100
20

50

15

0

10

−50
−100

5

−150
0
0

5

10

15

20

25

Altitude [km]Figure. 6.10: The improved version of Fig. 6.5

wind direction [degree]

wind directions

Chapter 6. On-Sky Evaluation with RAVEN

97

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0ms

33ms

66ms

100ms

0km

2km

4km

6km

8km

0ms

30
wind speeds

150

25Wind speed [m/s]

100
20

50

15

0

10

−50
−100

5

−150
0
0

5

10

15

20

25

Altitude [km]Figure. 6.11: The improved version of Fig. 6.6

wind direction [degree]

wind directions

Chapter 6. On-Sky Evaluation with RAVEN

1s

2s

3s

4s

5s

98

6s

7s

8s

9s

10s

0km

4km

8km

Figure. 6.12: Example of the phase distortion pattern averaged over each one second at
diﬀerent altitudes.

0.6
0.5

r0 from auto correlation
SLODAR
CFHT DIMM

r0zenith [m]

0.4
0.3
0.2
0.1
0
20:00 21:00 22:00 23:00 00:00 01:00 02:00 03:00 04:00
Time in HST, 2015/6/23-6/24

Figure. 6.13: r0 profiles in the zenith direction on June 23–24th, 2015. The blue plots
indicate the r0 values measure by the correlation SLODAR of RAVEN.
The red plots are the r0 values measured by the auto-correlation fittingof RAVEN. The gray lines are r0 values measured by the CFHT DIMM
(http://mkwc.ifa.hawaii.edu/).

Chapter 6. On-Sky Evaluation with RAVEN

99

2.5

Altitude [km]

15

2
10

1.5
1

5

0.5
0

2 (h)dh [m1/3 ]10−13 CN

30
20:00 21:00 22:00 23:00 00:00 01:00 02:00 03:00 04:00
Time in HST, 2015/6/23-6/24

2 (h) profile in the zenith direction on June 23–24th, 2015. The y-axis is
Figure. 6.14: CN
altitudes from the telescope.

30

Altitude [km]

20
10

15
10

5

Wind speed [m/s]

25

15

5
0

0
20:00 21:00 22:00 23:00 00:00 01:00 02:00 03:00 04:00
Time in HST, 2015/6/23-6/24

Figure. 6.15: Wind speed profile on June 23–24th, 2015. The y-axis is altitudes from thetelescope.

Chapter 6. On-Sky Evaluation with RAVEN

100

Altitude [km]

15

100
50

10

0
−50

5

−100

Wind direction [degree]

150

−150

0
20:00 21:00 22:00 23:00 00:00 01:00 02:00 03:00 04:00
Time in HST, 2015/6/23-6/24

Figure. 6.16: Wind direction profile on June 23–24th, 2015. The y-axis is altitudes fromthe telescope.

6.2.3

Oﬀ-Line Evaluation of Multi Time-Step Reconstruction

Finally, we evaluate the performance of the multi time-step reconstruction with the
improved wind estimates. We saved the on-sky measurements both of the OL and
the CL-WFSs without an AO correction. The WFE can be computed by comparing the wave-front in the science direction reconstructed from the measurementsof the open-loop SH-WFSs with the wave-front computed from the measurementsof the closed-loop SH-WFSs.
In this subsection, we show the result of the oﬀ-line analysis with two diﬀerent
data set. The first one is shown in Fig. 6.17. Fig. 6.17 shows the asterism, CN2
profile, the wind profile, and the tip-tilt removed residual WFEs. The residual
WFEs are computed from the on-sky measurement of the SH-WFSs with thesingle and multi time-step reconstruction. The multi time-step reconstructor iscomputed with the revised wind speeds and directions and diﬀerent ∆t. The
turbulence spreads over the wide range of altitudes. The wind speeds are 1 m
s−1 at 0 km, and 8 m s−1 at 14 km. In this case, the multi time-step reconstructor
can decrease the WFE compared with the single time-step reconstructior. The
WFE of the multi time-step reconstruction decreases with ∆t up to ∆t=0.16 s,
and this trend is consistent with the results from the analytical model (Fig. 3.4)
and numerical simulation (Fig. 4.5). When ∆t is larger than 0.16 s, the WFE

Chapter 6. On-Sky Evaluation with RAVEN

101

starts to increase. This is because the eﬀect of the turbulence boiling aﬀects on
the multi time-step reconstructor as ∆t > 0.16 s. The expected decay ratio with∆t = 0.16 s for the fastest layer (7.16 m s−1 ) is fdecay =0.76. In other words,
the frozen flow assumption is valid when fdecay ≥ 0.76. Furthermore, there isa gain with the multi time-step reconstruction in this case, even if fdecay ∼ 0.8,
although the boiling of the turbulence aﬀects the performance of the multi timestep reconstructor. Therefore, our assumption in the previous chapters that thefrozen flow holds as fdecay ≥ 0.7 is valid.
Another case is shown in Fig. 6.18. In this case, the ground layer has a large
contribution to the total turbulence power, and it is expected to the gain withthe multi time-step reconstruction would be small compared with the first case,because the multi time-step reconstruction is more eﬀective to the uncovered andunoverlapped areas at high altitudes. In Fig. 6.18, the WFE of the multi timestep reconstructor is larger than that of the single time-step reconstructor. Thisdiﬀerence becomes larger with larger ∆t. The considerable reasons is the error of
the wind estimation.
Here, we discuss the eﬀect of dome seeing in the multi time-step reconstructor
as another reason making the diﬀerence in the performance of the multi time-stepreconstruction between the two cases. Guesalaga et al (2014) [21] founds that thedome seeing has a diﬀerent temporal decay pattern of the cross-correlation peakof the measured slopes from that of the atmospheric ground layer. They measure
the fraction of the dome seeing in the ground layer by using this temporal decay.
We compute this ground layer temporal decay from the temporal cross-correlationof the measured slopes for the two cases, and investigate the eﬀect of the dome
seeing. The intensity of the correlation peak is measured by summing correlationsin the central 3×3 pixel on the temporal cross-correlation maps.
Fig. 6.19 shows the decay ratio of the ground layer for the first case. This canbe fitted by only one linear line, and it suggests that the temporal decay of thecross-correlation can be describe with moving turbulence. The fitted line becomes
0 at the time delay of 1 s. Therefore, the ground layer is completely de-correlated
in 1 s or the peak corresponding to the ground layer goes to out of the central
3×3 box in 1s (i.e, the ground layer moves 1.5 pixel on the correlation map or 1.2
m on the aperture plane). The observed wind speeds at the ground layer, shownin Fig. 6.17, is 1.14 m s−1 , and the travel distance of the ground layer in 1 s is 1.14m, which corresponds to ∼1.5 pixel on the correlation map. Thus, the decay ratio
is consistent with the observed wind speeds.
On the contrast, the decay ratio of the ground layer consists of two slopes inChapter 6. On-Sky Evaluation with RAVEN

102

the second case. This feature is also observed by GEMS [21]. The steeper slope(blue) represents the temporal decay of the moving turbulence, and the moderateslope (green) represents the temporal decay of the dome seeing. If the dome seeing
decays linearly as the time delay is smaller than 0.8 s, the initial fraction of thedome seeing in the ground layer is 64 %. In the multi time-step reconstructor,
these two components moving diﬀerently at the same altitude can not be modeled
and, as a result, it causes the tomographic error. The increase in the WFE withthe multi time-step reconstructor in the second case shown in Fig. 6.18 is possibly
aﬀected by the dome seeing.
Such a dome seeing eﬀect will aﬀects not only the multi time-step reconstructor,
but also predictive controllers [10, 20]. In order to avoid the eﬀect of the dome
seeing, it is required to monitor the fraction of the dome seeing, for example bychecking the temporal decay of the ground layer, and use the multi time-stepreconstruction in the case without the dome seeing.

Chapter 6. On-Sky Evaluation with RAVEN

80

2.5

NGS 3	NGS 1	

Ch 1	

0

Ch 2	

−40

NGS 2	

r0 =0.247 m
1.5

Elevation=33.969 degree

1
0.5

−60
−80
−80 −60 −40 −20 00

20 40 60 80

0

2

4

X [arcseconds]
30

wind speeds
wind directions

25

Wind speed [m/s]

6

8

10

12

14

16

Altitude [km]

150
100

20

50

15

0

10

−50

5

−100

Wind direction [degree]−20

Fd38, 2015/6/23 19:55:45

240

CN2 [10−13 m1/3]

Y [arcseconds]

60

20

103

−150
0

0

2

4

6

8

10

12

14

16

Altitude [km]Tip-Tilt removed WFE [nm]

680
670
660
Ch1, Single

650

Ch1, Multi

640

Ch2, Single

630

Ch2, Multi

620
610
600
0.02 0.04 0.06 0.08

0.1

0.12 0.14 0.16 0.18

0.2

∆t [s]
2 profile (top right), the wind profile (middle), andFigure. 6.17: Asterism (top left), CN
the tip-tilt removed residual WFE (bottom). The residual WFEs are computed from the on-sky measurement of the SH-WFSs with the single andmulti time-step reconstruction. The multi time-step reconstructor is computed with the revised wind speeds and directions and diﬀerent time delay.

Chapter 6. On-Sky Evaluation with RAVEN

80

2.5

NGS 1	

Ch 2	

0

Ch 1	

−20

NGS 3	NGS 2	

−60
−80
−80 −60 −40 −20 0

r0 =0.150 m
1.5

Elevation=62.996 degree

1
0.5
0

20 40 60 80

Fd38, 2015/6/27 01:02:19

2

0

2

4

X [arcseconds]
30

wind speeds
wind directions

25

Wind speed [m/s]

6

8

10

12

Altitude [km]

150
100

20

50

15

0

10

−50

5

−100

Wind direction [degree]20

CN2 dh [10−13 m1/3]

40

−150
0

0

2

4

6

8

10

12

14

16

Altitude [km]480

Tip-Tilt removed WFE [nm]Y [arcseconds]

60

−40

104

460
440

Ch1, Single
Ch1, Multi

420

Ch2, Single
Ch2, Multi

400
380
360
0.02 0.04 0.06 0.08

0.1

0.12 0.14 0.16 0.18

0.2

∆t [s]

Figure. 6.18: The same figure as the Fig. 6.17 with a diﬀerent case.

1416Chapter 6. On-Sky Evaluation with RAVEN

105

Measurement
Moving Turbulence1

Decay Ratio

0.8
0.6
0.4
0.2
0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Time Delay [s]Figure. 6.19: Decay ratio for the ground layer in the first case shown in Fig. 6.17. The redand blue lines indicate the computed decay ratio and a linear line for themoving component determined by the least squares method, respectively.

Measurement
Moving Turbulence
Dome Seeing

1

Decay Ratio

0.8
0.6
0.4
0.2
0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Time Delay [s]Figure. 6.20: Decay ratio for the ground layer in the first case shown in Fig. 6.18. The
red, blue and green lines indicate the computed decay ratio, a linear line forthe moving turbulence, and a linear line for the dome seeing, respectively.
The linear lines are computed by the least squares method with the timedelay range of [0 0.6] second for the moving turbulence and [0.8 1.5] second
for the dome seeing.

Chapter 7
Conclusion
In this thesis, we develop and evaluate a new tomographic reconstruction method
for the wide-field adaptive optics (WFAO), called multi time-step reconstruction,
to solve the lack of information in the tomographic reconstruction and expand
the size of the FoR of WFAO systems. Also, we evaluate a method to estimate
wind speeds and directions at multiple altitudes from the measurements of multiple wave-front sensors. This wind estimation method is required to implement
the multi time-step reconstruction. The basic idea of the multi time-step reconstruction is increasing the information by using the measurements at both of thecurrent and previous time-steps with using wind speeds and directions and thefrozen flow assumption.
First, we demonstrate the influence of the lack of information in the tomographic reconstruction. Then, we present that our tomographic reconstruction
method is eﬀective to reduce the tomographic error due to the lack of information. In the case of 3 NGSs and a telescope of 30 m circular aperture, the multitime-step reconstruction can expand the FoR size from 80 arcseconds diameterto 200 arcseconds diameter without increasing the tomographic error under theassumption that the turbulence follows the frozen flow assumption and we fullyknow the wind speeds and directions.
Second, we show the result of the numerical simulation of the tomographic
reconstruction on a MOAO system on the future ELT. The numerical simulationshows that the multi time-step reconstruction increases Strehl ratio (SR) over a
field of regard of 10 arcminutes diameter by a factor of 1.5–1.8 if we fully know thewind speeds and directions and the turbulence follows the frozen flow assumption.
Third, we present the laboratory of the multi time-step reconstruction method
and the wind estimation method by using RAVEN, which is a MOAO demonstrator. We can successfully measure the wind speeds and directions in the laboratory
experiment. Also, the multi time-step reconstructor can increase an EE in a 140

106Chapter 7. Conclusion

107

mas box by 0.03–0.05. However, absolute value of SR and EE is smaller thanthe value expected by the numerical simulation. Some part of the discrepancy
can be explained by the implementation error and the residual tip-tilt errors. In
order to understand what limits the performance of RAVEN, more experiments
are required.
Finally, we present the on-sky performance of the multi time-step reconstruction and the wind estimation with RAVEN n the Subaru telescope. Unfortunately,
there is no clear improvement with the multi time-step reconstruction because of
the error of the wind estimation. With the oﬀ-line analysis using the on-sky
measurement of the wave-front sensors, we improve the wind estimation method
and re-evaluate the multi time-step reconstruction. We conclude that the multitime-step reconstruction can reduce the WFE compared with the single time-stepreconstructor if fdecay > 0.7 %. Also, we conclude that if there are both of themoving turbulence and the dome seeing, the multi time-step reconstructor does
not work.

Appendix A
Mathematics
A.1

Minimum Variance Reconstructor and Analytical Tomographic Error

We drive the tomographic wave-front reconstructor Eq. (2.7). Substituting Eq.
(2.5), uβk = Eβk sα , and Eq. (2.10)into Eq. (2.6)ϵβk = ⟨||φβk − φ̂βk ||2 ⟩= ⟨||Pβk ϕ − Nβk ûβk (t)||2 ⟩= ⟨||Pβk ϕ − Nβk Eβk sα ||2 ⟩= ⟨||Pβk ϕ − Nβk Eβk (Γα Pα ϕ + η)||2 ⟩= ⟨||(Pβk − Nβk Eβk Γα Pα )ϕ + Nβk Eβk η||2 ⟩(A.1)

A variance of a column vector x is expressed as ||x||2 = Tr[xxT ], where Tr(X)indicates the trace of a matrix X. By using this expression and putting A =(Pβk − Nβk Eβk Γα Pα ) and B = Nβk Eβk , Eq. (A.1) can be transformed into[]
ϵβk = Tr ⟨(Aϕ − Bη)(Aϕ − Bη)T ⟩[][]= Tr A⟨ϕϕ⟩AT + Tr B⟨ηη⟩B T[][]= Tr AΣϕϕ AT + Tr BΣηη B T ,(A.2)

where ⟨ϕϕ⟩ = Σϕϕ and ⟨ηη⟩ = Σηη are covariance matrix for the phase distortiondue to the atmospheric turbulence and a covariance matrix of the measurement
noise, respectively. We assume that the phase distortion is independent from themeasurement noise , this is, ⟨ϕη⟩ = ⟨ηϕ⟩ = 0. We can compute the residual WFE
due to the tomographic error from Eq. (A.2) analytically by matrix calculation.
The reconstructor Eβk for βk direction is derived by minimizing ϵβk . This is108Appendix A. Mathematics

109

equal to computing Eβk , which meets ∂ϵβk /∂Eβk = 0. The partial diﬀerential of
the first term is computed as∂(AΣϕϕ AT )∂=−Pβ Σϕϕ PαT ΓTα EβTk NβTk
∂Eβk
∂Eβk k∂−Nβk Eβk Γα Pα Σϕϕ PβTk
∂Eβk∂+
Nβk Eβk Γα Pα Σϕϕ PαT ΓTα EβTk NβTk
∂Eβk
= − 2NβTk Pβk Σϕϕ PαT ΓTα + 2NβTk Nβk Eβk Γα Pα Σϕϕ PαT ΓTα=2NβTk (Nβk Eβk Γα Pα Σϕϕ PαT ΓTα − Pβk Σϕϕ PαT ΓTα ).

(A.3)

The partial diﬀerential of the second term is computed as∂(BΣηη B T )∂=Nβk Eβk Σηη EβTk
∂Eβk
∂Eβk
=2NβTk Nβk Eβk Σηη .

(A.4)

Then, combining Eq. (A.3) and Eq. (A.4) leads to[]∂ϵβk= 2NβTk Nβk Eβk (Γα Pα Σϕϕ PαT ΓTα − Σηη ) − Pβk Σϕϕ PαT ΓTα .
∂Eβk

(A.5)

The reconstructor meeting ∂ϵβk /∂Eβk = 0 is given as
Eβk =(NβTk Nβk )−1 NβTk Pβk Σϕϕ PαT ΓTα (Γα Pα Σϕϕ PαT ΓTα − Σηη )−1=Fβk Eϕ .

(A.6)
(A.7)

This is the minimum variance reconstructor for MOAO, which compute the DMcommands to compensate the phase distortion in βk direction. The first matrix
Fβk is a fitting matrix of a DM for βk direction, and the second matrix Eϕ
is a estimation matrix to reconstructs the phase distortion of the atmosphericturbulence layers. The estimation matrix can be transformed as−1 −1 T T −1
Eϕ = (PαT ΓTα Σ−1ηη Γα Pα − Σϕϕ ) Pα Γα Σηη .

A.2

(A.8)

The Scaling of Discrete Laplacian Matrix

In this section, we drive the scaling of discrete laplacian introduced in 2.1.3 basedon [20]. The phase covariance matrix of the turbulence layer indexed by i, Σϕϕ,i ,Appendix A. Mathematics

110

1	
 -4	
 1	
-4	
 2	

2	

2	
1	
1	
 -4	
 1	
1	
Figure. A.1: The

is approximated by using the discrete laplacian Li as Eq. (2.14). In the approximation, each row of the discrete laplacian indicates the local curvature of eachphase point on the turbulence layers, which is defined by surrounding four phase
points as the following stencil.

0 ω1 0
stencil(L̄i ) = ω2 ω3 ω4  .
0 ω5 0(A.9)

For the phase points which is not on the boundary, the weights are ω3 = −4 andω1 = ω2 = ω4 = ω5 = 1. These weights are diﬀerent for the boundary phase points
as shown in Fig. A.1.
The scaling constant ci for the discrete Laplacian is determined to satisfyT⟨ϕi Σϕϕ,i ϕi ⟩ = ⟨ϕTi LTi Li ϕi ⟩ = c2i ⟨ϕTi L̄Ti L̄i ϕi ⟩, where ϕi is a phase distortion
vector for ith turbulence layer. If the number of phase points on ith turbulencelayer is denoted by nilayer , the right term can be computed as ⟨ϕTi Σϕϕ,i ϕi ⟩ = nilayer .
Therefore, the scale constant is given byc2i =nilayer.
⟨ϕTi L̄Ti L̄i ϕi ⟩(A.10)

Denoting kth row vector of L̄i by {l̄i }k , we have
nilayer⟨ϕTi L̄Ti L̄i ϕi ⟩ =

∑ ⟨(
)2 ⟩
T
ϕ {l̄i }k(A.11)

k=1

Since each row vector of the discrete Laplacian
contains
no more than 5 non-zero⟨(⟩)
2
element, which defined by the stencil, ϕT {l̄i }kfor non-boundary phase point

Appendix A. Mathematics

111

is computed as⟨(

ϕT {l̄i }k

)2 ⟩

=

⟨( 5
∑

)2 ⟩
ϕ j ωj

j=1

⟨
⟩
1 ∑
ωj ωj ′ (ϕj − ϕj ′ )2
2 j,j ′ =1
5

=−

1 ∑
ωj ωj ′ Dϕ (ξj,j ′ )=−2 j,j ′ =1
5

(A.12)

where Dϕ is a phase structure function from Eq. (1.7) for the Kolmogorov modelor Eq. (1.14) for the von Karman⟨( model) and⟩ xii,j ′ is a spatial separation between
2′Tphase point j and j . Then, ϕ {l̄i }kfor the boundary phase points can becomputed similarly. Finally, the scale constant ci can be obtained by combining
Eq. (A.10), Eq. (A.11), and Eq. (A.12).

Reference[1] M. Akiyama, S. Oya, Y. H. Ono, H. Takami, S. Ozaki, Y. Hayano, I. Iwata,
K. Hane, T. Wu, T. Yamamuro, and Y. Ikeda, “Tmt-age: wide field of regard
multi-object adaptive optics for tmt,” Proc. SPIE 9148, 914814–914814–14(2014).
[2] S. M. Ammons, L. Poyneer, D. T. Gavel, R. Kupke, C. E. Max, and L. Johnson, “Evidence that wind prediction with multiple guide stars reduces tomographic errors and expands moao field of regard,” Proc. SPIE 8447, 84471U(2012).
[3] D. R. Andersen, S. S. Eikenberry, M. Fletcher, W. Gardhouse, B. Leckie,
J.-P. VÃľran, D. Gavel, R. Clare, R. Guzman, L. Jolissaint, R. Julian, andW. Rambold, “The moao system of the irmos near-infrared multi-object spectrograph for tmt,” Proc. SPIE 6269, 62694K–62694K–12 (2006).
[4] D. R. Andersen, K. J. Jackson, C. Blain, C. Bradley, C. Correia, M. Ito,
O. LardiÃĺre, and J.-P. VÃľran, “Performance Modeling for the RAVEN
Multi-Object Adaptive Optics Demonstrator,” Publications of the Astronomical Society of the Pacific 124, pp. 469–484 (2012).
[5] P. S. Argall, O. N. Vassiliev, R. J. Sica, and M. M. Mwangi, “Lidar measurements taken with a large-aperture liquid mirror. 2. sodium resonancefluorescence system,” Appl. Opt. 39, 2393–2400 (2000).
[6] J. Beckers, “Increasing the size of the anisoplanatic patch with multiconjugate
adaptive optics,” Proc. ESO conference 4007, 693–703 (1988).
[7] R. J. Bouwens, G. D. Illingworth, P. A. Oesch, M. Trenti, I. LabbÃľ,L. Bradley, M. Carollo, P. G. van Dokkum, V. Gonzalez, B. Holwerda,
M. Franx, L. Spitler, R. Smit, and D. Magee, “Uv luminosity functions at
redshifts z âĹĳ 4 to z âĹĳ 10: 10,000 galaxies from hst legacy fields,” The
Astrophysical Journal 803, 34 (2015).
[8] T. Butterley, R. W. Wilson, and M. Sarazin, “Determination of the profile
of atmospheric optical turbulence strength from SLODAR data,” Monthly
Notices of the Royal Astronomical Society 369, 835–845 (2006).
[9] A. Consortini and L. Ronchi, “Choice of the model of atmospheric turbulence,”
Appl. Opt. 11, 1205–1211 (1972).

112

REFERENCE

113[10] C. Correia, K. Jackson, J.-P. Véran, D. Andersen, O. Lardière, and
C. Bradley, “Static and predictive tomographic reconstruction for wide-field
multi-object adaptive optics systems,” J. Opt. Soc. Am. A 31, 101–113 (2014).
[11] C. M. Correia, K. Jackson, J.-P. Véran, D. Andersen, O. Lardière, and
C. Bradley, “Spatio-angular minimum-variance tomographic controller for
multi-object adaptive-optics systems,” Appl. Opt. 54, 5281–5290 (2015).
[12] A. CortÃľs, B. Neichel, A. Guesalaga, J. Osborn, F. Rigaut, and D. Guzman,“Atmospheric turbulence profiling using multiple laser star wavefront sensors,”
Monthly Notices of the Royal Astronomical Society 427, 2089–2099 (2012).
[13] G. P. David J. E. Floyd, Jo Thomas-Osip, “Seeing, wind, and outer scaleeﬀects on image quality at the magellan telescopes,” Publications of the Astronomical Society of the Pacific 122, 731–742 (2010).
[14] E. Diolaiti, “MAORY: A Multi-conjugate Adaptive Optics RelaY for the EELT,” The Messenger 140, 28–29 (2010).
[15] B. L. Ellerbroek, “Eﬃcient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A 19, 1803–
1816 (2002).
[16] S. G. Els, T. Travouillon, M. SchÃűck, R. Riddle, W. Skidmore, J. Seguel,E. Bustos, and D. Walker, “Thirty Meter Telescope Site Testing VI: Turbulence Profiles,” Publications of the Astronomical Society of the Pacific 121,
pp. 527–543 (2009).
[17] D. L. Fried, “Statistics of a geometric representation of wavefront distortion,”
J. Opt. Soc. Am. 55, 1427–1435 (1965).
[18] Gendron, E., Vidal, F., Brangier, M., Morris, T., Hubert, Z., Basden, A.,
Rousset, G., Myers, R., Chemla, F., Longmore, A., Butterley, T., Dipper,N., Dunlop, C., Geng, D., Gratadour, D., Henry, D., Laporte, P., Looker,
N., Perret, D., Sevin, A., Talbot, G., and Younger, E., “Moao first on-sky
demonstration with canary,” A&A 529, L2 (2011).
[19] L. Gilles and B. Ellerbroek, “Shack-hartmann wavefront sensing with elongated sodium laser beacons: centroiding versus matched filtering,” Appl. Opt.
45, 6568–6576 (2006).
[20] L. Gilles, P. Massioni, C. Kulcsár, H.-F. Raynaud, and B. Ellerbroek, “Distributed Kalman filtering compared to Fourier domain preconditioned conjugate gradient for laser guide star tomography on extremely large telescopes,”
J. Opt. Soc. Am. A 30, 898–909 (2013).
[21] A. Guesalaga, B. Neichel, A. CortÃľs, C. BÃľchet, and D. GuzmÃąn, “Using
the Cn2 and wind profiler method with wide-field laser-guide-stars adaptive

REFERENCE

114

optics to quantify the frozen-flow decay,” Monthly Notices of the Royal Astronomical Society 440, 1925–1933 (2014).
[22] W. GÃďssler, S. Rabien, S. Esposito, M. Lloyd-Hart, L. Barl, U. Beckmann, T. Bluemchen, M. Bonaglia, J. L. Borelli, G. Brusa, J. Brynnel,
P. Buschkamp, L. Busoni, L. Carbonaro, C. Connot, R. Davies, M. Deysenroth, O. Durney, R. Green, H. Gemperlein, V. Gasho, M. Haug, P. Hubbard,
S. Ihle, M. Kulas, R. Lederer, J. Lewis, C. Loose, M. Lehmitz, J. Noenickx,
E. Nussbaum, G. Orban de Xivry, D. Peter, A. Quirrenbach, M. Rademacher,
W. Raab, J. Storm, C. Schwab, V. Vaitheeswaran, and J. Ziegleder, “Statusof the argos ground layer adaptive optics system,” Proc. SPIE 8447, 844702–
844702–10 (2012).
[23] F. Hammer, F. SayÃĺde, E. Gendron, T. Fusco, D. Burgarella, V. Cayatte, J. Conan, F. Courbin, H. Flores, I. Guinouard, L. Jocou, A. LanÃğon,G. Monnet, M. Mouhcine, F. Rigaud, D. Rouan, G. Rousset, V. Buat, andF. Zamkotsian, “The falcon concept: multi-object spectroscopy combinedwith mcao in near-ir,” Proc. ESO workshop p. 139 (2002).
[24] J. W. Hardy, Adaptive Optics for Astronomical Telescopes (Oxford University
Press, 1998).
[25] Y. Hayano, H. Takami, S. Oya, M. Hattori, Y. Saito, M. Watanabe, O. Guyon,
Y. Minowa, S. E. Egner, M. Ito, V. Garrel, S. Colley, T. Golota, and M. Iye,
“Commissioning status of Subaru laser guide star adaptive optics system,”
Proc. SPIE 7736, 77360N (2010).
[26] T. M. Herbst, R. Ragazzoni, A. Eckart, and G. Weigelt, “Linc-nirvana: achieving 10 mas imagery on the large binocular telescope,” Proc. SPIE 7014,
70141A–70141A–6 (2008).
[27] G. H. Herriot, D. R. Andersen, J. Atwood, C. Correia, P. Byrnes, C. Boyer,
K. Caputa, J. Dunn, B. Ellerbroek, J. Fitzsimmons, L. Gilles, P. Hickson,
A. Hill, J. Pazder, V. Reshetov, S. Roberts, M. Smith, J.-P. Veran, L. Wang,and I. Wevers, “Nfiraos; tmt early light adaptive optics system,” Imaging andApplied Optics p. AWA4 (2011).
[28] K. Jackson, C. Correia, O. LardiÃĺre, D. Andersen, C. Bradley, L. Pham,
C. Blain, R. Nash, D. Gamroth, and J.-P. VÃľran, “Tomography and calibration for raven: from simulations to laboratory results,” Proc. SPIE 9148,
91482K–91482K–13 (2014).
[29] N. Kobayashi, A. T. Tokunaga, H. Terada, M. Goto, M. Weber, R. Potter,
P. M. Onaka, G. K. Ching, T. T. Young, K. Fletcher, D. Neil, L. Robertson, D. Cook, M. Imanishi, and D. W. Warren, “IRCS: infrared camera andREFERENCE

115

spectrograph for the Subaru Telescope,” Proc. SPIE 4008, 1056–1066 (2000).
[30] A. N. Kolmogorov, “Dissipation of energy in the locally isotropic turbulence,”
Akademiia Nauk SSSR 32, 16 (1941).
[31] M. Lamb, D. R. Andersen, J.-P. VÃľran, C. Correia, G. Herriot, M. Rosensteiner, and J. Fiege, “Non-common path aberration corrections for current
and future ao systems,” Proc. SPIE 9148, 914857–914857–13 (2014).
[32] O. LardiÃĺre, D. Andersen, C. Blain, C. Bradley, D. Gamroth, K. Jackson,
P. Lach, R. Nash, K. Venn, J.-P. VÃľran, C. Correia, S. Oya, Y. Hayano,
H. Terada, Y. Ono, and M. Akiyama, “Multi-object adaptive optics on-sky
results with Raven,” Proc. SPIE 9148, 91481G (2014).
[33] J.-F. Lavigne, F. Lamontagne, G. Anctil, M. Wang, M. Tremblay,
O. LardiÃĺre, R. Nash, D. Andersen, M. Savard, P. CÃťtÃľ, C. H. Bradley,and F. ChÃćteauneuf, “Design and test results of the calibration unit for the
MOAO demonstrator RAVEN,” Proc. SPIE 8447, 844754 (2012).
[34] C. Lawson and R. Hanson, Solving Least Squares Problems (Prentice-Hall,1974).
[35] L. H. Lee, “Sparse-matrix regularization for minimum-variancereconstructionof pseudo-kolmogorov turbulence,” in “Adaptive Optics: Analysis andMethods/Computational Optical Sensing and Imaging/Information Photonics/Signal Recovery and Synthesis Topical Meetings on CD-ROM,” (OpticalSociety of America, 2007), p. JTuA2.
[36] L. Lessard, M. West, D. MacMynowski, and S. Lall, “Warm-started wavefront
reconstruction for adaptive optics,” J. Opt. Soc. Am. A 25, 1147–1155 (2008).
[37] J. V. M. Azouit, “Optical turbulence profiling with balloons relevant to astronomy and atmospheric physics,” Publications of the Astronomical Societyof the Pacific 117, 536–543 (2005).
[38] P. Massioni, C. Kulcsár, H.-F. Raynaud, and J.-M. Conan, “Fast computationof an optimal controller for large-scale adaptive optics,” J. Opt. Soc. Am. A
28, 2298–2309 (2011).
[39] A. Miyashita, N. Takato, T. Usuda, F. Uraguchi, and R. Ogasawara, “Statistics of the weather data, environment data, and the seeing of the Subaru
Telescope,” Proc. SPIE 5489, 207–217 (2004).
[40] B. Neichel, F. Rigaut, F. Vidal, M. A. van Dam, V. Garrel, E. R. Carrasco,
P. Pessev, C. Winge, M. Boccas, C. d’Orgeville, G. Arriagada, A. Serio,
V. Fesquet, W. N. Rambold, J. LÃĳhrs, C. Moreno, G. Gausachs, R. L.
Galvez, V. Montes, T. B. Vucina, E. Marin, C. Urrutia, A. Lopez, S. J.

REFERENCE

116

Diggs, C. Marchant, A. W. Ebbers, C. Trujillo, M. Bec, G. Trancho, P. McGregor, P. J. Young, F. Colazo, and M. L. Edwards, “Gemini multiconjugate
adaptive optics system review âĂŞ ii. commissioning, operation and overall performance,” Monthly Notices of the Royal Astronomical Society 440,
1002–1019 (2014).
[41] A. M. Obukhov, “Structure of the temperature field in a turbulent flow,”
Akademiia Nauk SSSR 13, 58 (1949).
[42] L. Poyneer, M. van Dam, and J.-P. Véran, “Experimental verification of thefrozen flow atmospheric turbulence assumption with use of astronomical adaptive optics telemetry,” J. Opt. Soc. Am. A 26, 833–846 (2009).
[43] F. J. Rigaut, B. L. Ellerbroek, and R. Flicker, “Principles, limitations, andperformance of multiconjugate adaptive optics,” Proc. SPIE 4007, 1022–1031(2000).
[44] F. Rigaut, B. Neichel, M. Boccas, C. d’Orgeville, F. Vidal, M. A. van Dam,
G. Arriagada, V. Fesquet, R. L. Galvez, G. Gausachs, C. Cavedoni, A. W.
Ebbers, S. Karewicz, E. James, J. LÃĳhrs, V. Montes, G. Perez, W. N. Rambold, R. Rojas, S. Walker, M. Bec, G. Trancho, M. Sheehan, B. Irarrazaval,
C. Boyer, B. L. Ellerbroek, R. Flicker, D. Gratadour, A. Garcia-Rissmann,and F. Daruich, “Gemini multiconjugate adaptive optics system review âĂŞ i.
design, trade-oﬀs and integration,” Monthly Notices of the Royal Astronomical Society 437, 2361–2375 (2014).
[45] G. Rousset, T. Fusco, F. Assemat, E. Gendron, T. Morris, C. Robert, R. Myers, M. Cohen, N. Dipper, C. Evans, D. Gratadour, P. Jagourel, P. Laporte,
D. Le Mignant, M. Puech, H. Schnetler, W. Taylor, F. Vidal, J.-G. Cuby,
M. Lehnert, S. Morris, and P. Parr-Burman, “Eagle moao system conceptual design and related technologies,” Proc. SPIE 7736, 77360S–77360S–11(2010).
[46] M. Schöck and E. J. Spillar, “Method for a quantitative investigation of thefrozen flow hypothesis,” J. Opt. Soc. Am. A 17, 1650–1658 (2000).
[47] V. I. Tatarski, Wave Propagaton in a Turbulent Medium (McGraw-Hill, 1961).
[48] A. Tokovinin, “From diﬀerential image motion to seeing,” Publications of the
Astronomical Society of the Pacific 114, pp. 1156–1166 (2002).
[49] C. R. Vogel and Q. Yang, “Fast optimal wavefront reconstruction for multiconjugate adaptive optics using the fourier domain preconditioned conjugate
gradient algorithm,” Opt. Express 14, 7487–7498 (2006).
[50] R. W. Wilson, “SLODAR: measuring optical turbulence altitude with a

REFERENCE

117

ShackâĂŞHartmann wavefront sensor,” Monthly Notices of the Royal Astronomical Society 337, 103–108 (2002).
[51] S. A. Wright, E. J. Barton, J. E. Larkin, A. M. Moore, D. Crampton, and
L. Simard, “The infrared imaging spectrograph (iris) for tmt: sensitivities andsimulations,” Proc. SPIE 7735, 77357P–77357P–13 (2010).