Beam shaping for laser-based adaptive optics in astronomy

The availability and performance of laser-based adaptive optics (AO) systems are strongly dependent on the power and quality of the laser beam before being projected to the sky. Frequent and time-consuming alignment procedures are usually required in the laser systems with free-space optics to optimize the beam. Despite these procedures, significant distortions of the laser beam have been observed during the first two years of operation of the Gemini South multi-conjugate adaptive optics system (GeMS). A beam shaping concept with two deformable mirrors is investigated in order to provide automated optimization of the laser quality for astronomical AO. This study aims at demonstrating the correction of quasi-static aberrations of the laser, in both amplitude and phase, testing a prototype of this two-deformable mirror concept on GeMS. The paper presents the results of the preparatory study before the experimental phase. An algorithm to control amplitude and phase correction, based on phase retrieval techniques, is presented with a novel unwrapping method. Its performance is assessed via numerical simulations, using aberrations measured at GeMS as reference. The results predict effective amplitude and phase correction of the laser distortions with about 120 actuators per mirror and a separation of 1.4 m between the mirrors. The spot size is estimated to be reduced by up to 15% thanks to the correction. In terms of AO noise level, this has the same benefit as increasing the photon flux by 40%.


The need for laser beam optimization
The performance of adaptive optics (AO) is strongly dependent on the quality of the wavefront sensing measurements from the guide stars because they directly affect the correction error [1]. For the Shack-Hartmann (SH) wavefront sensors, currently the better adapted sensor for laser guide stars on large telescopes, the standard deviation measurement error (σ meas ) is a function of the spot size and the number of photons received per subaperture and per frame according to [1]: where θ image is the full width at half maximum (FWHM) of the spot image in the focal plane of a SH subaperture, and N ph is the average number of photons received in a SH subaperture each frame. Therefore, the smaller the spot image or the greater the number of received photons, the better the measurement accuracy. This article analyzes the case of the Gemini South multiconjugate AO system (GeMS) [2] which uses a single sodium laser to generate a constellation of 5 laser guide stars on sky. The size of the spot image θ image depends on various parameters, among which the laser spot size in the sodium layer, the elongation of the seen spot due to the sodium layer thickness, the diffraction produced by the SH subaperture size and the atmospheric seeing.
For the first two years of operations of GeMS during the season of low sodium concentration (December-January), the system frame rate had to be frequently reduced compared to its maximum of 800 Hz in order to guarantee the average detection of at least 135 photons per subaperture per frame from each laser guide star [3,4]. This minimum flux is required for adequate closed-loop AO correction by GeMS. From Eq. (1), if the image spot size θ image is reduced by 15%, it is equivalent for a given AO noise level to increase the number of photons by 40%. Laser-based AO systems usually being designed with little margin in terms of expected photons return, it appears all the more important to minimize the spot size.
In laser specifications, the smallest laser spot possible is obtained if the 589 nm beam to be projected has a beam quality factor M 2 = 1 [5]. Any discrepancy from this ideal shape results in M 2 values greater than 1. An estimation of the M 2 of a laser can be computed using a SH to measure the beam irradiance and wavefront at waist [6] and by using the discrete Fourier transform of the resulting complex field, the far-field distribution can also be computed. From each of these fields, the standard deviation of the x and y axis margin distribution of the beam intensity, can be obtained. Using the notations D x and θ x for this quantity along x−direction at the waist and in the far-field respectively, the M 2 is estimated by During the 2011/2012 commissioning campaign, the measured M 2 factors on GeMS were below 1.5 and 1.3, in x and y respectively [4]. A beam quality M 2 close to 1 has not been achieved so far with GeMS laser. The M 2 values even worsened during the first half of 2013, revealing significant degradation of the laser beam quality with time. The laser beam has been measured right at the output of the laser bench using a dedicated Shack-Hartmann sensor. The obtained irradiance maps are shown in Fig. 1 for sets of measurements taken in April and October 2013. Note that these measurements were both made after several days of alignment optimization, i.e. on the best laser shape at each time. From the data of April 2013, M 2 factors of 2.21 and 1.21 along x and y directions respectively were found. In October 2013, the corresponding M 2 factors were 2.23 and 1. 15.
These values highlight the strong distortion of the amplitude along x-axis, also visible in all plots of Fig. 1. The dotted line represents the half of the maximum to quantify the FWHM and it shows that the beam diameter is between 1.5 and 2 times larger along x than along y. The beam quality is always worse in the x-axis with GeMS because it corresponds to the unguided dimension of the amplification stages in the waveguide amplifier module [4]. In addition, a slow increase of the M 2 values over the length of the runs has also been noticed, suggesting that the system suffers from constant misalignment [4].
The impact of this degraded shape of the beam on the laser intrinsic spot size has been previously quantified by including the measured irradiance of April 2013 (top of Fig. 1) in simulations of the uplink propagation of the laser from the launching telescope to the mesosphere [7]. No phase aberrations were considered, but the beam was simulated for amplitude shapes different from that of a perfect Gaussian shape. The same uplink simulation code is used to evaluate the impact of amplitude and phase distortions on the short-exposure spot size at the sodium layer. Previous on-sky measurement of the laser intrinsic spot size of 1.0" has  been reported for GeMS [8]. This number cannot be directly compared to the values obtained with the simulation because there is not a strictly linear dependence of the laser return flux and the short-exposure irradiance focused at the sodium altitude. However, GeMS laser is not expected to induce saturation at the sodium layer, so the variance of the short-exposure irradiance over the SH subaperture field-of-view provides a good estimate of the spot size at a first order [9]. In addition, the vertical distribution of the sodium layer is not taken into account in this computation, but its effect could be considered to globally increase of the spot size by 10% [9]. Simulating a seeing of 0.73 arcseconds at 500 nm (median seeing over 33 GeMS observations nights between December 2012 and June 2013 [3]), on averaged over 100 successive independent atmospheres, the obtained FWHM of the short-exposure spot (arcseconds) is: • 0.61" for a perfect Gaussian beam with FWHM of 25 cm. For a given AO noise level (in Eq. (1)), the reduced size of the short-exposure laser spot with a perfect Gaussian beam compared to the distorted beam would be equivalent to an increase of the received number of photons by 30 to 45%. The beam degradation is understood as the result of quasi-static aberrations in both the laser bench and the beam transfer optics (BTO) path [4]. First, the co-alignment of the optical system for sum frequency generation used in GeMS is highly subject to temperature hysteresis, as well as telescope vibrations. This critical alignment tends to drift along a run. Secondly, the BTO is in a closed environment, with no turbulence, but the outside temperature can vary by 10 to 15 degrees, so aberrations along the BTO are subject to slow evolutions. Frequent and constraining calibrations and alignment procedures (quasi-static aberrations) are thus required in both parts to optimize the beam to launch. These complex and time-consuming alignment procedures currently used will strongly reduce the availability of GeMS [3]. It may also affect the future generations of AO requiring high power laser and using beam transport by mirrors, as planned for the Thirty Meter Telescope [10]. It is worth noting that recent research and development of fiber lasers for sodium guide stars [11] have shown very promising laboratory results to produce lasers with better M 2 quality and high power. There is still few on-sky experience to guarantee that this new technology can solve all the problems for an AO system like GeMS, requiring a 50-Watt laser. Therefore the study of an automatic technique to improve the beam shape of the current laser in GeMS is considered highly desirable.
The concept of a beam shaping system for laser-based astronomical AO was recently proposed [12] and its potential benefit on the GeMS signal-to-noise ratio has been highlighted with preliminary simulations of the uplink propagation of the laser toward the sodium layer [7]. This concept is studied here in greater detail, quantifying with numerical simulations the expected performance of the method with respect to design parameters. These results will be further used to develop an experimental prototype for such beam correction to be tested on GeMS laser.
In Sect. 2, the architecture with two deformable mirrors for amplitude and phase correction is detailed. In Sect. 3, the iterative algorithm used to deduce the commands for both mirrors is described with emphasis given to the iterative unwrapper required to solve the problem in practice. Last in Sect. 4, preliminary simulations of the beam shaping correction method are presented. The results confirm the ability to correct both amplitude and phase with this type of architectures. The influence of parameters like the separation distance between the mirrors, the number of iterations of the algorithm and the number of actuators of the mirrors are emphasized.

2-DM concept for laser beam shaping
In order to optimize the laser beam to be launched, correction of both amplitude and phase distortions of the beam is required. This is referred to as optical field conjugation or amplitudephase correction [13].
The presented concept takes its origins in the beam shaping research for laser radar and active coherent imaging (e.g. [14][15][16][17][18]). In these first studies, the common goal was to produce an optical device which modifies the beam coming from a given source to lead to a desired irradiance distribution in the near-field. In order to avoid loss of energy, only refractive or reflective surfaces must be used to redistribute the energy. Early works were focused on static shaping and lead to the difficult manufacturing of very specific optical components [15][16][17]. Eismann et al. [18] proposed a more versatile solution where they transform a Gaussian beam into a uniformly illuminated field using two computer-generated holographic elements. During the past two decades, several wavefront correction technologies were developed to match this low loss of energy and to provide also adaptive correction for evolving distortions introduced by propagation through turbulent media. Among these, the micro-electro-mechanical systems (MEMS) deformable mirrors (DM) [19] are now compact and robust device for beam shaping at least cost. In our study, we foresee to use two MEMS DMs.  In the fields of defense industry and free space communications, the beam propagation through the atmosphere leads to strong amplitude distortions degrading the system performance. During the last two decades, a significant amount of research was generated for these applications to develop adaptive optics systems correcting for both amplitude and phase distortions of a beam [20][21][22][23][24] .
One of the first approach to tackle phase and amplitude correction simultaneously with two DMs was proposed by Kanev et al. [13]. The first DM corrects for the amplitude and the second one corrects for the phase. The concept has later been revised by Roggemann et al. [20,25,26] in a configuration where the second mirror is in the far-field of the first one. The scheme of such correction system is presented in Fig. 2(a), and it is referred to as mono-static because the receiving and projecting apertures are the same. The phase and amplitude corruption induced by the atmosphere in the optical path from the target above the aperture is supposed to be measured with a suitable device. It could be a SH sensor, from which outputs need to be processed to compute the received amplitude and the phase of the field in the aperture plane of the telescope. The measurements are analyzed to deduce the commands to be sent to the DMs in order to shape the laser beam to the conjugate of the received field, and project it. First, the laser collimated beam is supposed to fall upon DM1, where a phase shape φ 1 is added to the wave. This applied phase aims at modifying the amplitude shape of the wave in the Fraunhofer region, after the Fourier transforming mirror, when the wave is incident on DM2. The second DM, DM2, is optically conjugated to the pupil plane of the beam projecting aperture, so it only modifies the phase distribution of the output beam. A novel modification of the latter approach to near-field correction has been proposed by Barchers et al. [21,[27][28][29] showing the potential interest of conjugating the first DM to a finite range where a turbulent layer could be present, instead of using a Fourier transforming mirror. This kind of dual-DM systems to precompensate atmospheric turbulence effects on beacon has been used during the last two decades in the military and free-space communication fields.
In the adaptive optics systems for astronomy, only few studies tackled full optical field conjugation, because phase correction is by far more important than amplitude correction [30]. The use of multiple DMs in AO for astronomy are common for phase correction only and in tomographic configurations, aiming to correct the phase over a wide field of view like GeMS [2]. Nevertheless, amplitude distortions or scintillation, being a source of loss of contrast in astrophysical imaging, Gonsalves suggested [31] to correct for amplitude distortions thanks to a second DM. Again for higher contrast purposes, phase-induced amplitude apodization [32] is an active are of research using aspheric mirrors for exo-planet imaging. And in the same field, multiple DMs have recently be shown to provide polychromatic correction and allow to null the intensity in High Strehl images [33,34]. With respect to the lasers for astronomical AO, an on-going project studies the potential benefit of only phase correction before projection to the mesosphere [35,36]. More recently another concept of laser beam cleanup with 2 DMs has been presented [37] but again it targets only phase distortions since the second DM is conjugated to the first one and is used only to increase the maximum spatial frequency of the correction and the stroke of the overall system. To the knowledge of the authors, the present study of amplitude and phase correction of the laser beam before projection for astronomical AO is new. The 2DM-system proposed here for astronomical purposes uses a second DM conjugated to a finite range, z (see Fig. 2(b)). DM2 is in the near-field of the DM1 (i.e. in the Fresnel region), at a separation distance z of the order of 1 meter. DM2 remains conjugated to the pupil plane of the beam projecting aperture as in Fig. 2(a)). The near-field configuration should allow matching the space constraints of the laser systems on the telescope, although maintaining the ability to correct for amplitude and phase of the field.
The primary objective of the 2-DM concept is to provide correction of quasi-static amplitude and phase aberrations affecting the beam before its projection. The distortions of GeMS laser such as the ones presented in Fig. 1 have been observed to only evolve over a time scale of hours or days [38]. The developed solution focuses only on these quasi-static distortions but the choice of deformable mirrors would allow in a further extension of the study to tackle optimizations at shorter time scales. For instance, it may be used to finely tune the focus of the output beam when variations of the sodium layer height occurs. The DMs could also precompensate for the dynamical aberrations induced by the atmospheric turbulence in the propagation path of the laser above the launch telescope. This last extension is far from being straightforward however since the turbulence of interest is not directly measured. In astronomical laser-based systems, the configuration would usually be bistatic (cf. Fig. 2(b)), meaning that the laser launching telescope (left) is not the telescope measuring the atmospheric turbulence (right). For on-axis launching of the laser, the turbulence above the secondary mirror, i.e. behind the central obscuration of the laser-based AO system of the telescope, might be extrapolated for low-order modes corrected by the ground-conjugated DM of the AO. The atmospheric wavefront distortion would thus be known thanks to the adaptive optics of the telescope (represented by AO sensing device in Fig. 2(b)), for which the receiving aperture is the primary mirror and does not coincide with the launching telescope. Demonstration of such capability is yet unknown to the authors, so such dynamical extension of the correction is left for further research. In Fig. 2(b), the output of the AO sensing device is connected with a dashed arrow to the processing computer illustrating the fact that this option is not considered further in this paper.
The correction of quasi-static amplitude and phase distortions is based on a given desired beam shape at the projection aperture, U t . The phase patterns to be applied to the DMs are determined by a phase retrieval method [39] like in most of the amplitude-phase correction systems mentioned before. The principle of the method is to iteratively estimate in two planes the phase shapes which match the amplitude and phase constraints in other planes (e.g. image plane, conjugated plane), as far as they are linked by a propagation relation. The first phase retrieval methods were developed by Gerchberg and Saxton [40], Gonsalves [41] and Fienup [42]. Fienup showed that these methods have the same local convergence properties around a fix point as the steepest descent and conjugated gradient algorithms [39]. The intensive research on optical field conjugation for laser defense applications led to eventually demonstrate [21] that algorithms for optical field conjugation, usually referred to as Gerchberg-Saxton algorithms, are applications of what is known as the method of sequential projections onto constraints sets [43]. The algorithm presented in Sect. 3 is a version of Gerchberg-Saxton method adapted to account for the near-field separation of the DMs and to unwrap the phases to be applied on DMs.
Last, to bring the 2DM correction concept to a practical implementation the method needs to provide phase estimates that can be efficiently reproduced by the mirrors. Non-segmented mirrors do not perform well when phase cuts are present [44]. Phase cuts can occur because the phase estimate is wrapped and because of branch points [45]. As shown later with simulations in Sect. 4, the estimated phases in our context are very likely to be wrapped and to generate branch points where the field amplitude vanishes in the aperture. To avoid phase cuts, we developed a regularized iterative unwrapper (see Sect. 3.3) to be applied to the phase estimates after every iteration of the phase-retrieval algorithm. The use of least squares and weighted least squares unwrappers have been reported in the literature related to beam shaping with DMs [29,44]. The least squares approach is not adapted to cases where the field intensity vanishes. A weighted least squares unwrapper allows to deal with degeneracies of the phase estimates where the intensity is zero, and thus avoid phase cuts. However, it is difficult to implement a weighted least squares unwrapper in practice because the weight based on the intensity distribution of the computed propagated field changes at every iteration of the phase retrieval algorithm. So the pseudo-inverse matrix required to apply the weighted least squares cannot be precomputed only once. To overcome these difficulties, a novel iterative unwrapper is detailed in Sect. 3.3. Its weighting can be easily updated on the fly with every new field estimate and no precomputation of pseudo-inverse matrix is required.

Algorithm for amplitude and phase correction
For the 2-DM design presented in Fig. 2(b), the following notations are used. The incident field on DM1 is the one coming out from the laser bench and is noted where x is the coordinate vector in the plane of DM1, and u l (x) and φ l (x) account for the aberrated amplitude and phase in this plane respectively. The field leaving DM1 is noted U 1 and is expressed as where φ 1 is the phase pattern applied to DM1 and m 1 represents the reflection mask for the mirror, in order to account for the finite extension of the device. The physics between DM1 and DM2 is modeled with the Fresnel integral as an approximation of the propagation. This free-space propagation of the electromagnetic wave U 1 , at distance z, is described by the linear unitary transformation T z such that [27,46] where f x is the spatial frequency vector, λ is the laser wavelength and U z is the incident field on DM2. After DM2 phase correction φ 2 , the output field is where m 2 represents the mask of this mirror pupil. DM2 is conjugated to the projecting aperture such that U 2 is assumed to be the transmitted field. In our context, the optical field desired to be transmitted by the launching aperture is assumed known, and we refer to it as the desired output field U t , noted where u t is its amplitude and φ t its phase. This desired output field is a constraint in the algorithm. Its choice is discussed in [7] and in the present paper we only consider it as being known and fixed.
In laser applications where strong scintillation is produced by the turbulence during the propagation above the launching aperture, the transmitted field U t needs to be the conjugate of the field which would be received in a mono-static configuration like in Fig. 2(a) [27]. In the astronomical context, the weak turbulence above the telescope leads to negligible amplitude distortions, such that the sensing device located after the receiving aperture in Fig. 2(b) is aimed at analyzing the phase distortions φ b only.
Assuming U t and U l known, the goal of the 2DM concept is to determine φ 1 and φ 2 to be applied to DM1 and DM2 respectively in order to obtain the field U t on the projecting aperture [20]. For this, a phase retrieval method is used, iteratively providing estimates of φ 2 and φ 1 to satisfy the intensity constraints fixed by U l and U t respectively. The main steps of the algorithm, similar to the one described by Roggemann and Lee [20] are described in Sects. 3.1 and 3.2. In practice, the algorithm must also unwrap the phase estimates in order to guarantee efficient continuous-mirror correction. A specific iterative unwrapper is developed for this purpose in Sect. 3.3.

Forward iterative step for DM2
Given an initial estimate for φ 1 , the incoming field on DM2 is computed Since we want the field after DM2 to be equal to the desired output U t , the required phase to apply on DM2 φ 2 is constrained to satisfy Once the phase φ 2 is obtained from Eq. (9), we compute an additional parameter α, which is the ratio of the incident energy E 2 on DM2 over the energy E t of the desired output field U t . This factor α is used later in Eq. (11) to rescale the field to be propagated backward in the iterative step for determination of φ 1 . The value of α is mathematically defined by

Backward iterative step for DM1
This second step estimates the phase φ 1 to be applied to DM1 [20]. It is obtained by back propagation of the rescaled desired output field. The incident field on DM2 is thus and the input field on DM1 can be written as The last equality is obtained by the property of the adjoint operator T * z of T z . Indeed, for any given complex field U, T * z (U) = T −z (U * ). The estimation of phase φ 1 is obtained applying the constraint that the phase and amplitude after backpropagation beyond DM1 matches the input laser field phase φ l and amplitude u l . Thus, Equations (8)-(10) and Eqs. (11)-(13) represent the core of the 2-DM correction algorithm, being computed alternatively and iteratively, in order to obtain the best amplitude and phase outputs.
It is important to observe that the computation of Eqs. (9) and (13) leads to wrapped phase values. In addition, the field amplitude for which the argument must be obtained can be null or numerically very close to zero at some locations within the DMs aperture, which could lead to undesired branch points [45] appearing in the phase estimations of φ 1 and φ 2 . In order to obtain phase functions φ 1 and φ 2 likely to be efficiently reproduced by the continuous sheet of a DM, the phases must be unwrapped and branch points must be avoided.

Phase unwrapping
In order to overcome these difficulties of wrapped phase and branch points, an iterative and regularized unwrapper estimator is presented. This unwrapper is modified from the weighted least squares unwrapper presented in [12], and is obtained by iteratively solving the linear system where G is the discrete representation of the gradient operator as described by Fried [45]. The PV symbol stands for the principal-value operator which output is an angle that falls in the range ±π (modulus 2π), W g is a diagonal array of weights, C is a discrete representation of the curvature operator and µ is a scalar weighting the curvature regularization term C T C.
A wrapped phase has the gradient of its unwrapped phase plus or minus a multiple of 2π. This is why the unwrapper in Eq. (14) involves the gradient operator on both sides of the equation and the principal value operator on the right hand side where the wrapped phase is. In addition, φ 2 and φ 1 are computed from Eqs. (9) and (13) only where the field amplitude u in the argument is numerically non zero. This leads in both cases to wrapped phase maps φ , and the wrapped phase samples are arbitrarily set to zero where the amplitude is also zero, i.e. where degeneracy occurs.
The matrix W g is used as a diagonal weighting related to the level of confidence of the gradient value of the wrapped phase, i.e. accounting for the fact that the computed phase and gradients values do not represent the true values where the field amplitude reaches zero. The weighting is defined differently for gradient in x and gradient in y, such that the component in (i, j) associated to gradient along x at location (x i , y j ), can be mathematically described as with ε a small numerical value, e.g. 10 −23 , representing the threshold of the product u(x i+1 , y j ) u(x i , y j ) of two neighbour amplitude samples above which the associated gradient is considered with a strictly positive weight. This ε is chosen as a negligible value compared to the average field amplitude over the aperture. This logarithmic scale has the benefit to spread the weighting values of W g mainly in the range 1 to 10 where non negligible amplitude values have been obtained. The weight become significantly lower as soon as one field amplitude involved in the gradient computation is zero, revealing the proximity of a non-illuminated area. The weight is exactly zero where the gradient computation involves two field samples of null amplitude. With such approach, the weighting matrix W g is updated on the fly, everytime the unwrapper is called with a new wrapped phase and amplitude field. This update is easily handled since there is no requirement to precompute a matrix inverse to apply the iterative unwrapper.
Another interesting feature of the unwrapper is the regularizing term C T C (in Eq. (14)), which favors solutions ofφ with low curvature. The relative importance of this penalty term can be tuned thanks to the scalar µ in Eq. (14). In practice, we have found that µ of the order of 10 −3 is apropriate to influence the solution only in areas of negligible amplitude.

Convergence criteria
The convergence of the beam shaping iterative method is assessed thanks to two complementary criteria. The first one ε 2 ampl measures the relative distance between the currently achieved output amplitude shape and the desired one [20]. After i iterations of the method, where the notation U (i) 2 stands for the achieved output field after DM2, when the i-th estimations of φ 1 and φ 2 are applied to DM1 and DM2 respectively. The scalar ε 2 ampl thus represents a relative amplitude error with respect to the initial amplitude error.
The second criterion is denoted as S and it is sensitive to the phase matching between the reached output field and the desired output one. It is written where index (i) again refers to the i−th iteration of the procedure. This criterion is inspired from the far-field intensity criterion used in the literature [27,47]. These two criteria are used in the simulations of Sect. 4 to enhance the convergence of the algorithm. They both evolve in the range [0, 1]. ε 2 ampl starts at 1 by definition and decreases toward 0 when the correction improves. On the contrary, S is expected to increase with iterations toward unity.

Simulations
In this section, we present numerical simulations of the optical field conjugation method proposed in the previous section. In a first step, we check that the optimization of the amplitude and phase outputs is effective, following the algorithm presented in Sect. 3. Next, the influence on performance of system parameters is studied: the separation distance z between the mirrors; the number of iterations; and the grid sampling step of the propagated fields and phase estimates.

Amplitude and phase correction
As a first step, we consider an arbitrary amplitude distortion of the incoming laser [12]. No phase aberrations are considered, but a particular output phase shape, a defocus, is specified. The laser wavelength is 589 nm and the distance between the DMs is chosen to be z = 1.8 m.
Without any beam shaping correction, i.e. φ 1 = φ 2 = 0, the input distortions lead to output amplitude very similar to the input one, slightly modified by the propagation in the near-field (z = 1.8 m). This output amplitude is represented in Fig. 3(a), over the 1 cm diameter of the aperture of DM2. Cuts along x and y central axes are represented with solid curves in Fig. 3(b). It can be observed that the amplitude of the beam is far from being Gaussian.
The desired output amplitude is set to a Gaussian shape, represented by dotted slices along x and y central axes in Fig. 3(b). The desired output phase has been arbitrarily chosen to be a defocus (see dotted curve in Fig. 4(b)). After 400 iterations of the phase retrieval method presented in Sect. 3, the output amplitude and phase after DM2 almost perfectly match the desired ones. The two-dimensional representation of the output amplitude is shown in Fig. 3(c), and the one of the output phase in Fig. 4(a). The central slices along x and y are plotted with solid curves, in Fig. 3(d) for the amplitude and in Fig. 4(b) for the phase. The convergence of the algorithm along the 400 iterations is illustrated in Fig. 5. Both corrections of amplitude and phase work perfectly. It is worth noting that output phase estimated in Fig. 4 differs from the desired phase shape near the aperture edges. This is due to the fact that the dotted curve represents a given phase function, the desired output. On the contrary, the solid curve represents an estimate of the resulting phase, computed from the simulation of the field propagation. Therefore, in the same way as explained in Sect. 3.3 for φ 1 and φ 2 computations, the estimation of the resulting phase is initially wrapped and only at non zero amplitude locations. The field amplitude (shown in Fig. 3(c)) for radii greater than 4 mm is so small that the phase cannot be estimated with accuracy at these locations. The solid curve in Fig. 4 is thus the result of the iterative unwrapper applied to the output phase estimate. In the outer ring area, the regularization of the iterative phase unwrapper detailed in Sect. 3.3 mainly drives the estimation of the phase. However, whatever the phase estimate in this area, it has no effect on the beam shaping performance because the field amplitude is almost zero. Note also that the criterion S in Fig. 5 is not notably affected by the phase discrepancy in this area because of the amplitude weighting in Eq. (17).
From Fig. 5, it can be concluded that the beam shaping has worked to perfection after 150 iterations. The correction of the output phase of the beam is almost instantaneous, as demonstrated by the value of S already jumping to 85% after a single iteration. Amplitude correction drives the global convergence. The influence of the DMs separation distance on the convergence speed is further studied in Sect. 4.2.
In this first simulation, the grid sampling size of φ 1 and φ 2 is identical to the resolution δ of the introduced aberrations U l and of the desired output field U t , with δ 161 µm. The desired output amplitude is a perfect Gaussian with FWHM of 1.7 mm, and the DMs pupils are considered to be 1 cm in diameter. The phase distortions φ 1 and φ 2 which should be applied on DM1 and DM2 in order to obtain such correction are presented in Fig. 6.
These smooth shapes are otbained thanks to the regularized unwrapper described in Sect. 3.3, applied to the wrapped phases computed from iterative procedures of Sects. 3.1 and 3.2.
The action of the phase unwrapper is illustrated in Fig. 7, where the central cuts along y for  φ 1 is shown in dotted line for the wrapped solution (between −π and π) and in solid line once unwrapped. The limit of the mirrors pupil diameter are represented by the vertical dashed lines. Again it can clearly be noticed that in the center of the slice (where the amplitude of the field is larger), the unwrapper accurately unwraps the phase estimates. On the pupil edges however, where the amplitude vanishes, the regularizing term of the unwrapper becomes significant and the unwrapped estimate is extrapolated, penalizing local curvature of the phase.

Influence of DM separation
In this section, the influence of the separation distance between the 2 DMs is investigated. The idea of the near-field position of DM2 comes from possible space constraints on the telescope. The simulations presented in this section demonstrate that such near-field propagation is able to produce satisfactory amplitude corrections for the laser beam.
The simulated system is identical to the one presented in Sect. 4.1, but a flat output phase is sought, in order to focus the attention on the amplitude correction abilities. The final correction performance is evaluated in terms of the M 2 factor of the output beam. This criterion is applicable in this case because the desired output amplitude is a Gaussian shape. The results of the correction performance after 400 iterations are presented in Fig. 8.
The output M 2 x and M 2 y are shown in solid and dashed lines as a function of the separation distance z between the DMs. The dotted horizontal lines stand for the M 2 x and M 2 y values when no beam shaping correction system exists. These values are computed assuming that the output beam is U l . When the beam shaping system is considered, the M 2 factor is computed after the DM2, which means that even if no DM correction is applied, the output beam is slightly affected in the near-field. Note that for z ≤ 1 m, M 2 y factor achieved after 400 iterations is worse than when no correction is applied. This highlights the difficult amplitude correction at short DMs separation. For distances greater than 2 m, it can be observed that the beam correction is very efficient, with resulting M 2 factor values as low as 1.02. Such values mean very good quality Gaussian beams. In addition to the evolution of the correction quality with the DMs separation z, its influence on the convergence speed is also studied. In Fig. 9 (left), the relative amplitude error of the output beam ε 2 ampl is represented as a function of the iteration number. The slowest convergence is obtained for the shortest distance z = 1 m, where convergence is reached after more than 350 iterations. The convergence is accelerated progressively when the separation z is increased. For z = 5 m, less than 20 iterations suffice. The convergence acceleration is due to the fact that a phase deformation of minor amplitude is required when the separation distance z increases. This is illustrated on the right of Fig. 9, where the highest curve is the phase deformation φ 1 required for separation z = 1 m. The required deformation is reduced by approximately a factor equal to p when the separation is changed from z = 1 m to z = p, with p in the range of 1 to 5 meters. It is thus an important parameter to take into account in the design of the beam shaping system, depending on the mechanical properties of the DMs to be used.

Influence of number of degrees of freedom
In practice, expensive costs of the DMs with large number of actuators may force to consider the application of phase deformations with a lower spatial resolution than the one considered in Sects. 4.1 and 4.2. So a study of the impact of the phase grid sampling for the estimates φ 1 and φ 2 on the performance of the beam shaping is required. We specify this constraint in terms of the number of estimated phase samples across the DMs diameter, n, such that a DM with ∼ 3/4n 2 degrees of freedom over the pupil area is expected to be able to reproduce it satisfactorily.
We choose the laser amplitude distortions u l as the ones measured from GeMS in April and October 2013 (Fig. 1). Again, we do not consider any phase aberrations φ l for the input. The simulation grid sampling has been chosen equal to the resolution of the laser wavefront sensor measurement at GeMS (Fig. 1), i.e. δ = 198 µm. The DMs aperture diameter is asumed to be 4.9 mm. The desired output amplitude is set to a perfect Gaussian with FWHM of 1.7 mm and a flat output phase.
The correction performance assessed by ε 2 ampl as a function of z is represented in the plots of Fig. 10 for various values of n. Input amplitude distortions measured in April 2013 are used for the plot on the left and the ones measured in October 2013 are used on the right plot. With a large number of degrees of freedom for the correction, it is possible to reduce the error to a few percentage points of its original value placing the DMs at short distances (z below 1 m). At such short separations, the error reduction significantly degrades when the number of phase samples n across the pupil is too small.
If the mirrors have only 11 actuators across the aperture for instance, the remaining error is about 10% of its original value, and it requires greater separation distance z, i.e. beyond 1.2 m. Note that the distortions measured in October 2013 appear more difficult to correct than the ones observed in April 2013. For separations z beyond 2 meters, there is no remarkable benefit in using 31 × 31 phase samples compared to using 11 × 11.    Fig. 11. Same as in Fig. 10 but the error is presented as an absolute value, not relatively to ε 2 ampl (0). The dashed lines represent ε 2 ampl (0) as a function of z.
It is worth noting that a minimum for the relative error exists within the DM separation range studied. For instance, if n = 31, the relative error ε 2 ampl with z = 2 m, is degraded by about 10 percentage points compared to the achieved error with z = 80 cm. However, this is mainly due to the normalization value in the denominator of ε 2 ampl (Eq. (16)), which decreases with z. In order to clarify this, the same simulations results are presented in Fig. 11 but showing the amplitude error ε 2 ampl before its normalization by the error without correction in Eq. (16). The amplitude error without any beam shaping correction (dashed line) progressively decreases with increasing separation distance z between the DMs. This effect simply results from the propagation of the beam in the near-field. This absolute representation reveals that good correction quality is obtained if z is greater than 1.4 m for all three cases (n = 11, 17 and 31). If n is greater than 11, shorter separations z can be chosen for approximately the same good correction quality.

Benefit of correcting GeMS amplitude and phase distortions
According to the results of Sects. 4.2 and 4.3, we are now interested in quantifying how much the 2-DM correction system could reduce the size of the laser spot image. For this, the measurements made at GeMS in April and October 2013 are considered with both their amplitude and wavefront aberrations. The quality of the correction is evaluated by simulating the uplink propagation of the resulting beam at the launching aperture to the sodium layer as mentioned in Sect. 1. The spot size at 90 km and its intensity are compared to the ones evaluated without correction. The reduction of the laser spot size at the sodium layer because of the correction is represented as an equivalent increase of the number of photons (in Eq. (1)) for a given AO noise level. These factors of increase are shown in Fig. 12 as a function of the simulated atmospheric seeing. Left graph stands for the aberrations measured in April 2013, and right graph is for the ones of October 2013. The reduction of the spot size, in terms of photons increase, produced by a perfect Gaussian beam is also plotted.  First, one observes that the reduction of the spot size is more significant in good seeing conditions. The laser beam aberrations thus degrade the AO performance more severely when the turbulence conditions are good.
Second, because the 2-DM correction does not lead to a perfect Gaussian beam the three curves are not superimposed with the upper one (Gaussian). Still, with a correction system using 11 × 11 mirrors, the equivalent increase of photons is in the range 20 to 40%. A higher-order correction with 17 actuators across the aperture improves the photons increase to the range of 25 to 45%. Obviously the correction of GeMS laser aberrations will not allow to double the AO frame rate (which would require an equivalent increase by 100%), but still for median seeing of observing conditions for GeMS an increase of the number of photons by 20 to 35% would already be achievable. In the most favorable conditions for a correction with 11 actuators across the diameter, the frame rate could even be increased by a factor 1.4, which could mean a jump from 350 Hz to almost 500 Hz.

Conclusion
A two-deformable-mirror concept is developed to optimize the beam shape of the laser used as beacons for adaptive optics in astronomy. The optimization system aims at correcting quasistatic beam distortions in amplitude and phase thanks to a phase retrieval technique. Nearfield algorithms previously developed in the literature for defense and space communications applications are adapted to the projection of lasers in the astronomical context. In particular, an iterative and regularized unwrapping method is developed to deduce commands for the mirrors. The method is computationally efficient and robust to handle any input amplitude distortions and any desired output amplitude, thanks to its iterative feature and on-the-fly update of the weighting matrix used for the unwrapping.
Preliminary simulations are presented to illustrate the effective amplitude and phase correction achieved by the method. This correction reduces the laser spot size thus improving the AO measurement accuracy. The convergence is shown to be accelerated when the separation distance between the 2 DMs increases. Regarding the number of degrees of freedom required for the mirrors, a reasonable order of 11 actuators across the aperture is shown to provide good beam shaping correction when applied to cases of amplitude and phase distortions measured at GeMS during 2013. Additional degrees of freedom on the deformable mirrors would improve the correction quality. Since the algorithm converges in a few seconds, it can definitely be applied to the the quasi-static aberrations observed at GeMS, which evolve at the time-scale of hours or days.
The simulated cases highlight the fact that for reduced number of degrees of freedom the considered examples of distortions are better corrected with DMs separations greater than 1.4 meter approximately. The simulation of uplink propagation toward the sodium layer show that the correction of aberrations measured on GeMS laser will not allow to double the AO frequency. Still in median seeing conditions for GeMS observations a 2-DM correction system with 11 actuators across the aperture may allow a reduction of the spot size by 10 to 15% compared to the spot generated by the laser beam with aberrations. For a given AO noise level, this is equivalent to an increase of the number of received photons by up to 40%. Space constraints within the telescope and costs of DMs will drive the trade-off between DM separation and number of degrees of freedom when designing the prototype to be tested on GeMS.
Finally, the measured quasi-static aberrations of the laser shown in this paper are the ones observed after alignment and manual optimization procedures on GeMS laser beam. This corresponds to a time-consuming process usually involving 1 or 2 laser engineers for at least 1 week before every run of GeMS. The presented 2-DM correction can relax the need for such complex alignment optimizations. Actually the laser optimization procedure at GeMS includes a trade-off between beam quality and laser power. So relaxing in a first step the pressure on laser beam quality could lead to improved laser power [38].

Aknowledgments
This work was supported by the Chilean Research Council (CONICYT) through grant Fondecyt 1120626. Benoit Neichel aknowledges the French ANR program WASABI -ANR-13-PDOC-0006-01. Particular thoughts and aknowledgements go to our co-author Vincent Fesquet who recently passed away.