A Frame Decomposition of the Atmospheric Tomography Operator

We consider the problem of atmospheric tomography, as it appears for example in adaptive optics systems for extremely large telescopes. We derive a frame decomposition, i.e., a decomposition in terms of a frame, of the underlying atmospheric tomography operator, extending the singular-value-type decomposition results of Neubauer and Ramlau (2017) by allowing a mixture of both natural and laser guide stars, as well as arbitrary aperture shapes. Based on both analytical considerations as well as numerical illustrations, we provide insight into the properties of the derived frame decomposition and its building blocks.


Introduction
The imaging quality of earthbound astronomical telescopes like the Extremely Large Telescope (ELT) [11] of the European Southern Observatory (ESO), currently under construction in the Atacama desert in Chile, suffers from aberrations due to turbulences in the atmosphere, which result in blurred images. This is commonly counteracted by the use of Adaptive Optics (AO) systems, which use the measurements of one or more Wavefront Sensors (WFS) to suitably adjust Deformable Mirrors (DMs) in such a way that the incoming wavefronts are corrected (flattened) after reflection on the mirrors; see Figure 1. 1 (left). Since the atmosphere is constantly changing, this has to be done in real time. For more details on adaptive optics we refer to [7,29,30].
There are a number of different types of AO systems, the simplest one being Single Conjugate Adaptive Optics (SCAO). Thereby, the light of a so-called Natural Guide  [2]) and sketch of an SCAO system (right, image taken from [6]).
Star (NGS), a bright star in the vicinity of an object of interest, is used to adjust the DM to obtain a corrected wavefront and thus a sharp image of the NGS and the nearby object. See Figure 1.1 (right) for a schematic drawing of an SCAO system. In case that there is no bright enough NGS available in the vicinity of the object of interest, or if one wants to achieve a good correction over either a larger or multiple fields of view, one needs to resort to different, more complex AO systems. A good unidirectional correction in the absence of a NGS in the vicinity of an object of interest is for example achieved by Laser Tomography Adaptive Optics (LTAO), while a correction over a large or multiple fields of view is achieved by Multiconjugate Adaptive Optics (MCAO) and Multiobject Adaptive Optics (MOAO), respectively [1,5,19,24,28]. These are schematically depicted in Figure 1.2. In common with all of those different AO systems is the use of so-called Laser Guide Stars (LGS), artificial stars created by powerful laser beams in the sodium layer of the atmosphere, which are used to increase the number of guide stars available for correction and thus enhance the imaging quality.
Since an SCAO system succeeds at enhancing the imaging quality in one direction of view only, measurements from a single WFS are enough to compute suitable correction shapes of the DM, because wavefront aberrations from two objects close to each other, in this case from the reference NGS and the considered object of interest, are approximately the same. However, this is not the case for LTAO, MCAO, and MOAO, where the NGS and LGS are far away from the object of interest, or a good correction has to be achieved over a large or multiple fields of view. Hence, one has to use multiple WFSs (typically one for each guide star) and DMs, whose correcting shapes have to be determined from the turbulence profile of the atmosphere, which in turn has to be calculated from the WFS measurements. This gives rise to the problem of the atmospheric tomography. LGS, respectively. The light blue areas correspond to those directions of view, for which the AO systems aim at achieving a correction. Image taken from [2].
Unfortunately, and in particular since the separation of the NGSs and LGSs is low (e.g., 1 arcmin for MCAO and 3.5 arcmin for MOAO), the problem of atmospheric tomography falls into the category of limited-angle tomography, which is known to be a severely ill-posed problem [4,21]. In addition, the number of available guide stars is relatively small as well (e.g., 6 LGSs for the ELT), which in combination with the severe ill-posedness makes the reconstruction of the full atmospheric turbulence above the telescope a hopeless endeavour. Hence, one works with the commonly accepted assumption that the atmosphere contains only a limited number of turbulent layers, which are infinitely thin and located at predetermined heights. The problem of atmospheric tomography then becomes the task of reconstructing the turbulences on only a finite number of turbulence layers from the available WFS measurements. For a schematic depiction with three layers, (natural) guide stars, and WFSs see Figure 1.

(left).
A number of numerical reconstruction approaches have been proposed and developed for the atmospheric tomography problem, among them a minimum mean square error method [13], a back-projection algorithm [14], conjugate gradient type iterative reconstruction methods with suitable preconditioning [9,15,18,36,37], the Fractal Iterative Method (FrIM) [33][34][35], the Finite Element Wavelet Hybrid Algorithm (FE-WHA) [20,[38][39][40], and Kaczmarz iteration [27,31]; see also [8,16,17,23,25,26,32] and the references therein. All of these methods work comparatively well, each with its own peculiar advantages and drawbacks, and the resulting reconstructions have been successfully used to enhance the overall imaging quality of the corresponding AO systems. However, these numerical approaches themselves do not provide any deeper insight into the atmospheric tomography problem itself.
Hence, the authors of [22] set out to provide a mathematical analysis of the atmospheric tomography operator (defined below) underlying the problem, which is derived from the Radon transform using the layered and limited-angle structure of the problem [7,13]. In particular, they derived a singular-value-type decomposition of the operator, which not only provides the basis for efficient numerical reconstruction methods Figure 1.3: Schematic depiction of an atmospheric tomography setup with three turbulence layers, (natural) guide stars, and wavefront sensors (left). The coloured areas are those parts of the turbulence layers which are "seen" by each of the different wavefront sensors. Illustration of the cone effect (right) for a single LGS. Images taken from [38].
but also allowed to gain insight into the ill-posedness of the problem itself. In contrast to the already known singular-value decomposition of the limited-angle tomography operator [4,21], the singular values of this decomposition could be computed explicitly.
However, the singular-value-type decomposition derived in [22] is only valid under a very restrictive set of assumptions. In particular, leaving aside some technicalities, it is only valid for square aperture shapes and tomography settings with only NGSs and no LGSs. This is obviously problematic for two reasons: firstly, the aperture shapes of telescopes is usually not square, and secondly, as we have already seen above, most of the AO systems which rely on atmospheric tomography include LGSs as an integral part of their design. Hence, many of the interesting theoretical results derived in [22] no longer hold for those (practically) important cases. Furthermore, while numerical routines based on this singular-value-type decomposition can in principle be adapted via measurement extension to (partly) circumvent the restriction of square aperture shapes, an adaption to include also LGSs is not possible in any straightforward way. This is mainly due to the so-called cone effect: since a LGS is created by a powerful laser beam inside the sodium layer of the atmosphere, it is, in contrast to an "infinitely far away" NGS, located at a finite height. Thus, light travelling from the LGS to the telescope pupil passes through larger areas in lower atmospheric layers than in higher ones; see Figure 1.3 (right) for an illustration. Mathematically, this results in the addition of a layer and guide star dependent scaling parameter in the atmospheric tomography operator (see below), which causes a number of complications.
The aim of this paper is to overcome the restrictions of the singular-value-type decomposition in [22] of the atmospheric tomography operator noted above. In particular, we want to find a decomposition which allows both LGSs and arbitrary aperture shapes. This is done in two steps: First, we consider the case of a tomography setup with only LGSs and no NGSs. Setting aside for the moment considerations on the practicality of such a setup (e.g., the tip/tilt problem), this is at the same time both a completion and a natural extension of the results of [22] , and an ideal starting point for deriving the main results of this paper. Then, we provide a decomposition of the atmospheric tomography operator for general problem settings including both a mixture of NGSs and LGSs, as well as arbitrary aperture shapes. This decomposition is done in terms of a set of functions which form a frame, with important implications for both theoretical as well as numerical aspects of the tomography problem.
The outline of this paper is as follows. In Section 2 we describe the precise mathematical setting of atmospheric tomography considered in this paper, and in Section 3 we review some necessary material on frames in Hilbert spaces. In Section 4 we derive the decompositions for the atmospheric tomography operator mentioned above, first in the case of square domains and LGSs only, and then for a mixture of NGSs and LGSs as well as arbitrary aperture shapes. Section 5 presents some numerical results based on the obtained analytical derivations and is followed by a short conclusion in Section 6.

Mathematical Setting
In this section, we describe the precise mathematical setting of the atmospheric tomography problem considered in this paper. For the sake of consistency, we mainly use the same notations as in [22] .
The atmospheric tomography problem is a limited-angle tomography problem with only finitely many directions of view α g , g = 1, . . . , G, where G denotes the total number of NGSs and LGSs. The directions α g ∈ R 2 are such that if (α x g , α y g , 1) ∈ R 3 points from the center of the telescope to the guide star g, then α g = (α x g , α y g ) ∈ R 2 . We denote the telescope aperture with Ω A ⊂ R 2 and assume that the atmosphere contains L layers, where each layer is a plane at height h l , l = 1, . . . , L parallel to Ω A .
Since we do not only consider NGSs but also LGSs in this paper, we need to, in particular, take the cone effect into account (see above). For this, we need to define the parameters c l,g . Assuming that G NGS and G LGS denote the number of NGSs and LGSs, respectively, such that G = G NGS + G LGS , we set LGS denotes the height of the LGSs. Since h l < h LGS for all l = 1, . . . , L, we have that c l,g ≤ 1. Furthermore, since we assume that the h l are in ascending order and that h L < h LGS , we have that c l,g ≥ 1 − h L h LGS > 0. For every layer l, we can now define the domains where The domains Ω l are exactly those parts of the layers which are "seen" by the wavefront sensors (compare with Figure 1.3) and are therefore those parts of the atmosphere on which one can expect to reconstruct the atmospheric turbulences. Denoting by φ = (φ l ) l=1,...,L the turbulence layers and by ϕ = (ϕ g ) g=1,...,G the incoming wavefronts, the atmospheric tomography operator can now by defined as follows: On the definition and image spaces of A we define the canonic inner products Completely analogous to [22,Theorem 3.1], we have that the operator A with respect to the above scalar products is not compact, and hence, a singular system does not necessarily need to exist for A. However, the authors of [22] managed to derive a singular-value-type decomposition of what they called the periodic atmospheric tomography operatorÃ, given bỹ where Ω T := [−T, T ] 2 with T sufficiently large, such that which form an orthonormal basis of L 2 (Ω T ). However, the derived decomposition introduces artefacts due to the periodicity assumptions and in particular does not cover the case of LGSs and mixtures of NGSs and LGSs. Furthermore, in practice, wavefront measurements are not given on the extended domain Ω T but on Ω A only. Thus, for applying the decomposition derived in [22] , these measurements have to be extrapolated to Ω T in some way, which is also not desirable. Hence, in Section 4, we extend the singular-value-type decomposition from [22] by deriving a frame decomposition of the original operator A, based on the ideas of [22] but avoiding some of their shortcomings. For this, we first review some necessary materials on frames in Hilbert spaces below.

Frames in Hilbert Spaces
For the upcoming analysis, we need to recall some known results on frames in Hilbert spaces, which can for example be found in [3]. First, recall the definition of a frame: Now, for a given frame {e k } k∈K , one can consider the so-called frame (analysis) operator F and its adjoint F * , or synthesis-operator, defined by Due to (3.1) and the general fact that F = F * , there holds that Furthermore, one can define the operator S := F * F , i.e., and it follows that S is a bounded linear operator with B 1 I ≤ S ≤ B 2 I, where I denotes the identity operator. Furthermore, S is invertible and B −1 It follows that if one definesẽ k := S −1 e k , then the set {ẽ k } k∈K also forms a frame, with bounds It is also known that every f ∈ H can be written in the form It is not always possible to compute the dual frameẽ k explicitly. However, since it holds that [3] where R := I − 2 A+B S, the elements of the dual frame can be approximated by only summing up to a finite index N , i.e., (3.6) The induced error of this approximation is controlled by the frame bounds A, B, i.e., For a numerical implementation, (3.6) can also be written in the following recursive form, which allows for an efficient numerical implementation in practicẽ Although for frames the decomposition of f in terms of the functions e k is not unique, the representation in (3.4) is the most economical, in the sense the following Proof. See for example [3].
Note that if {e k } k∈K is a tight frame with bounds B 1 = B 2 = B, then we have that S −1 = B −1 I,ẽ k = e k /B and therefore Using these results, we now proceed to derive a frame decomposition of the atmospheric tomography operator A below.

Frame Decomposition
In this section, we first derive a singular-value-type decomposition of the operatorÃ in the case of only LGSs, following the ideas of [22] . Afterwards, we derive a frame decomposition of the operator A as defined in (2.1) where, in contrast to [22] , we do not restrict ourselves to only NGSs, but allow a mixture of both NGSs and LGSs, as well as (almost) arbitrary aperture shapes Ω A .

The Pure LGS Case
In this section, we consider the periodic tomography operator (2.3) from [22] , but now for the case that instead of only NGSs, we consider a setting using only LGSs. Since in this case c l,g = c l independent of the guide star direction g, the operatorÃ can now be written in the form where, as before, Ω T = [−T, T ], however now with T sufficiently large, such that We now derive a singular-value-type decomposition of the adapted operatorÃ using the ideas from [22] . First, due to the presence of the constants c l , it makes sense to use, for ever layer l, a different orthonormal basis of L 2 (c l Ω T ), namely the functions Any function φ l ∈ L 2 (c l Ω T ) can be written in the form and for φ = (φ l ) L l=1 ∈ L l=1 L 2 (c l Ω T ) we collect the expansion coefficients φ jk,l in the vectors φ jk := (φ jk,l ) L l=1 ∈ C L . With this, we now get Proposition 4.1. LetÃ be defined as in (4.1) and let the G×L matricesÃ jk be defined by Proof. Using the definition (4.1) ofÃ and the expansion (4.3), we get from which the assertion immediately follows after interchanging the series, which is allowed since the norm ofÃ jk is bounded independently of j, k.
As in [22] we consider the singular system for each of the matricesÃ jk , i.e., the vectors v jk,n ∈ C L , u jk,n ∈ C G and numbers σ jk,n , n = 1, . . . , r jk ≤ min{G, L} such that where the superscript H denotes the Hermitian adjoint (i.e. the complex transpose), r jk is the rank of the matrixÃ jk , and the σ 2 jk,n are the positive eigenvalues of the matrices A H jkÃ jk andÃ jkÃ H jk , respectively. Hence, the decomposition ofÃ in terms of the functions w jk is given by which is completely analogous to [22,Equation (10)], however with a different singular systems (σ jk,n , u jk,n , v jk,n ). Thus, for incoming wavefronts the best-approximate solution of the equatioñ Aφ = ϕ is given by whereÃ † denotes the Moore-Penrose generalized inverse ofÃ (see for example [10,22]), which is well-defined if and only if the following Picard-condition holds

The General Case
In this section, we derive a frame decomposition for the general atmospheric tomography operator A as defined in (2.1), i.e., for the case of arbitrary aperture domains Ω A and both NGSs and LGSs.
The main idea of [22] , where some of the ideas of the upcoming analysis are taken from, was to use the special properties of the functions w jk (2.4), in particular that they form an orthogonal basis of L 2 (Ω A ). Unfortunately, for general domains Ω A , the functions w jk do not necessarily form a basis of L 2 (Ω A ). However, they do form a (tight) frame in the sense of Definition 3.1, as we see in the following Lemma 4.2. Let w jk be defined as in (2.4). If T is large enough such that Ω A ⊂ Ω T , then the system {w jk } j,k∈Z forms a tight frame over L 2 (Ω A ) with frame bound 1, i.e., Proof. We start by definingψ : where we have used that since the w jk form an orthonormal basis of L 2 (Ω T ), they are also a tight frame with frame bounds B 1 = B 2 = 1. This proves the assertion. Now, since {w jk } jk∈Z forms a tight frame with bound 1, it is also its own dual frame. Hence, (3.4) implies that any function ψ ∈ L 2 (Ω A ) can be written in the form ψ(x, y) = j,k∈Z ψ, w jk L 2 (Ω A ) w jk (x, y) . (4.7) In particular, we have that the incoming wavefronts ϕ can be written in the form We also want to expand functions on L 2 (Ω l ) in terms of frames. It is not difficult to find frames for L 2 (Ω l ); for example, for large enough T , the sets {w jk } jk∈Z or {w jk,l } jk∈Z already form frames over L 2 (Ω l ). However, for the upcoming analysis, those frames do not satisfy the specific needs of the problem under consideration. Hence, we use a different, problem tailored frame, which we build from the functions w jk,lg (x, y) := c −1 l,g w jk ((x, y)/c l,g ) I c l,g Ω A +αgh l (x, y) , which we assume to hold from now on. That these functions can indeed form frames can be seen from the following Proposition 4.3. Let w jk,lg be defined as in (4.9) and let T be large enough such that (4.10) holds. Then, for fixed l, the set {w jk,lg } jk∈Z, g=1,...,G forms a frame over L 2 (Ω l ) with frame bounds 1 ≤ B 1 ≤ B 2 ≤ G.
Proof. Let ψ ∈ L 2 (Ω l ) be arbitrary but fixed and start by looking at where we have used that c l,g Ω A + α g h l ⊂ Ω l ⊂ c l,g Ω T . Now since for each (fixed) g the functions c −1 l,g w jk ((x, y)/c l,g ) form an orthonormal basis of L 2 (c l,g Ω T ) (and thus a tight frame with bound 1), we get that G g=1 jk∈Z which, together with (4.11), yields (c l,g Ω A + α g h l ), we get that Combing this together with (4.12), we get that which proves the assertion. Now, denoting for fixed l the dual frame of {w jk,lg } jk∈Z, g=1,...,G by {w jk,lg } jk∈Z, g=1,...,G , the above proposition together with (3.4) implies that any function φ l ∈ L 2 (Ω l ) can be written in the form φ l (x, y) = G g=1 jk∈Z φ jk,lgwjk,lg (x, y) , φ jk,lg = φ l , w jk,lg L 2 (Ω l ) . (4.13) Next, we want to find an expansion of Aφ in terms of the frame {w jk } jk∈Z . Due to (4.7) and since Aφ ∈ L 2 (Ω A ), such an expansion is given by w jk , (4.14) As we already know, even though this might not be the only possible expansion, it is the most economical one in the sense of Proposition 3.1. Furthermore, due to (3.9), this expansion allows us to express Aφ − ϕ in terms of the expansion coefficients of Aφ and ϕ, which is important for determining an (approximate) solution of Aφ = ϕ.
We now derive an explicit expression for the coefficients (Aφ) g , w jk L 2 (Ω A ) in terms of the coefficients φ jk,lg in the following Proof. Due to the definition of A, and setting r = (x, y) we have Using substitution in the integral yields c l,g Ω A +αgh l φ l (r)w jk (r/c l,g ) dr .
Since c l,g Ω A + α g h l ⊂ Ω l , due to (4.9) and (4.13), we get Combining the above results, we get c −1 l,g w jk (α g h l /c l,g )φ jk,lg .
which directly yields the assertion.
Hence, combining (4.14) and the above proposition, we get that Thus, if we define (with a slight abuse of notation), the vectors 17) and the C G×(L·G) matrices  The above expansion is obviously similar to (4.4) for the pure LGS case, however now with different (and slightly larger) matrices A jk . Hence, we can again consider the singular value decomposition of each of the matrices A jk , i.e., with another small abuse of notation, the vectors v jk,n ∈ C L·G , u jk,n ∈ C G and numbers σ jk,n , n = 1, . . . , r jk ≤ G such that v H jk,m v jk,n = δ mn , u H jk,m u jk,n = δ mn , σ jk,1 ≥ · · · ≥ σ jk,r jk > 0 , (4.20) to get, in complete analogy to (4.6), the following frame decomposition of the atmospheric tomography operator A defined in (2.1): (Aφ)(x, y) = jk∈Z r jk n=1 σ jk,n v H jk,n φ jk u jk,n w jk (x, y) . which we use to find a solution to the equation Aφ = ϕ in Theorem 4.8 below. Before we proceed to that though, we observe that the structure of the matrices A jk allows to compute their singular-value-decomposition explicitly, which leads to the following Proposition 4.5. Let (σ jk,n , u jk,n , v jk,n ) r jk n=1 , for j, k ∈ Z be the singular systems of the matrices A jk as defined in (4.20). Then there holds r jk = G , σ jk,n = L l=1 c −2 l,n , u jk,n = e n , (4.23) where e n denotes the n-th unit vector in C G .
Proof. Due to the structure of A jk , we have that from which the formula for σ jk,n , u jk,n = e n and r jk = G immediately follow. Furthermore, since there holds v jk,n = (1/σ jk,n ) A H jk u jk,n , the expression for (v jk,n ) immediately follows, which concludes the proof.
Using the explicit representation derived above, we immediately get the following Furthermore, again due to Proposition 4.5, we get that u H jk,g ϕ jk = ( e n ) H ϕ jk = (ϕ jk ) g , which, together with (4.25) and the definition of A immediately yields the assertion.
Concerning the well-definedness of A, we can now derive the following Lemma 4.7. Let ϕ ∈ L 2 (Ω A ) G and let A be defined as in (4.22). Then Aϕ is welldefined and there holds Proof. LetF l be the frame operator (compare with (3.2)) corresponding to the dual frame {w jk,lg } jk∈Z ,g=1,...,G of {w jk,lg } jk∈Z ,g=1,...,G and letF * l be its adjoint. Since the dual frame is also a frame but with inverse frame bounds, it follows from Proposition 4.3 together with (3 Hence, for any sequence a l = {a jk,lg } jk∈Z ,g=1,...,G ∈ 2 there holds which we now use for the choice where {ϕ jk } jk∈Z are the expansion coefficients of ϕ defined in (4.8). Since due to (4.24) and the choice of a l there holds (Aϕ) l =F * l a l , (4.30) we now obtain from (4.28) that Now summing over l we get that where the last equality follows from the definition (4.8) of the coefficients ϕ jk together with Lemma 4.2. Due to Proposition 4.5 and since 0 < c l,g ≤ 1 the singular values σ jk,g are independent of j, k and bounded away from 0. Hence, since ϕ ∈ L 2 (Ω A ) G , it follows that Aϕ is well-defined. Furthermore, the above estimate together with the explicit expression for σ jk,g from Proposition 4.5 implies that which, together with 0 < c l,g ≤ 1, immediately yields the assertion.
With this, we are now able to prove the main theorem of this section. Proof. Let ϕ ∈ L 2 (Ω A ) G and let ϕ jk,g = ϕ g , w jk L 2 (Ω A ) be its canonical expansion coefficients, collected in the vectors ϕ jk ; compare with (4.8). Furthermore, let φ ∈ D(A) and let φ jk,lg = φ l , w jk,lg L 2 (Ω l ) be its canonical expansion coefficients, collected in the vectors φ jk ; compare with (4.13) and (4.17). Due to Proposition 4.4 there holds, where the matrices A jk are as in (4.18). Since due to Lemma 4.2 the set {w jk } jk∈Z forms a tight frame over L 2 (Ω A ) with frame bound 1, it follows with (3.9) that Hence, together with G g=1 jk∈Z Thus, φ is a solution of Aφ = ϕ if and only if its expansion coefficients φ jk satisfy The solutions of these matrix-vector systems can be characterized via the singular systems (σ jk,n , u jk,n , v jk,n ) r jk n=1 of the matrices A jk . Since in Proposition 4.5 we showed that r jk = G and u jk,n = e n , it follows that at least one solution of (4.34) exists. Hence, by the properties of the SVD it follows that the vectors are the unique solution of (4.34) with minimal 2 norm. Hence, if φ can be found such that φ jk = A † jk ϕ jk , then for all other solutions ψ of Aφ = ϕ there holds (4.36) We now show that the choice φ := Aϕ is exactly such that Since due to (4.17) there holds φ jk,lg = (φ jk ) l+(g−1)L , this is equivalent to showing For this, note first that as in the proof of Corollary 4.6 we have that Hence, we now have to show that φ jk,lg = (a l ) jk,g , ∀ jk ∈ Z , 1 ≤ l ≤ L , 1 ≤ g ≤ G .
To do so, note that in (4.30) in the proof of Lemma 4.7 we saw that whereF l is the frame operator of the dual frame {w jk,lg } jk∈Z ,g=1,...,G . Hence, we obtain where we used the fact (see [3]) that F lF * l = P l , with P l denoting the orthogonal projector in 2 onto R(F l ) = R(F l ), and that by assumption a l ∈ R(F l ). However, since this implies that φ jk,lg = φ l , w jk,lg L 2 (Ω l ) = (F l φ l ) jk,g = (a l ) jk,g , it now follows that (4.37) holds, and thus Aϕ is a solution of Aφ = ϕ. Together with (4.36), this also implies the coefficient inequality (4.31). Note that due to Lemma 4.7, Aϕ is well-defined since ϕ ∈ L 2 (Ω A ) G . Finally, since φ jk,lg = φ l , w jk,lg L 2 (Ω l ) = (a l ) jk,g , it follows that the expansion 4.22 is the most economical one in the sense of Proposition 3.1, which concludes the proof.
Remark. In the above theorem we saw that Aϕ is a solution of the atmospheric tomography problem given that a l ∈ R(F l ), which can be interpreted as a consistency condition. On the other hand, if instead we only assume that Aφ = ϕ is solvable, then it follows from (4.33) that for any solution φ there holds A jk φ jk = ϕ jk , and thus that Since the complete set of solutions is given by φ † + N (A), where φ † ∈ N (A) ⊥ denotes the minimum-norm solution of Aφ = ϕ, and since by (4.33) for anyφ ∈ N (A) there holdsφ jk ∈ N (A jk ), it follows that Noting that this can be rewritten as and since due to (4.30) and (4.38) there holds (Aϕ) l =F * l A † jk ϕ jk l+(g−1)L jk∈Z ,g=1,...,G , it now follows from (4.39) that which implies that if a l / ∈ R(F l ) one can still consider Aϕ as an approximate solution.
Remark. The explicit representation of A given in (4.24) can be used as the basis of an efficient numerical routine for (approximately) solving the atmospheric tomography problem Aφ = ϕ. Replacing the infinite sum over the indices j, k by a finite sum, all one needs for implementation are the vectors ϕ jk and the evaluation of the functionsw jk,lg (x, y) on a predetermined grid. For this, note that since the functions {w jk,lg } jk∈Z,g=1,...,G form the dual frame of the functions {w jk,lg } jk∈Z,g=1,...,G , they can be efficiently numerically approximated via the approximation formula (3.8). Furthermore, the coefficients ϕ jk can be efficiently computed via (4.8) using the Fourier transform. Thus, since apart from ϕ jk all other quantities can be precomputed independently of the input data ϕ, the computation of Aϕ via formula (4.24) can be efficiently numerically implemented.
Remark. In Lemma 4.7, we have derived that A ≤ 1/L, i.e., that A is bounded. In particular, for any noisy measurement ϕ δ of the incoming wavefronts ϕ this implies which shows that the (approximate) solution of the atmospheric tomography problem via the application of the operator A is stable with respect to noise in the data. At first, this result seems counter-intuitive, since the atmospheric tomography operator is derived from the Radon transform under the assumption of a layered atmosphere. Hence, since inverting the Radon transform is known to be an ill-posed problem, one expects A to become unbounded with an increasing number of layers L. However, note that on the one hand Ay is only an approximate solution for a l / ∈ R(F l ), and on the other hand, for L going to infinity, the atmospheric tomography operator A as defined in (2.1) does not tend towards the (limited-angle) Radon transform.
To see this, note that the sum in the definition of A stems from the discretization of the line integrals of the Radon transform via the simple quadrature rule where the weights (x k − x k−1 ) were dropped. Assuming that (x k − x k−1 ) = 1/L, which would correspond to equidistant atmospheric layers, the above quadrature rule becomes which indicates that unless the atmospheric tomography operator is multiplied by 1/L in this case, it does not converge to the desired limit as the number of layers tends to infinity. Incidentally, if A is replaced by (1/L)A, then A has to be replaced by LA, and it then follows from (4.27) that resulting in an upper bound for the solution operator LA which tends towards infinity for an increasing number of layers L, as one would expected. Note that further adaptations to the tomography operator A would have to be made in order to ensure its convergence to the Radon transform in the limit for L → ∞. These correspond to the integration weights having to be adapted to different guide star directions. However, since this is practically not relevant as L is generally fixed, we do not go into further details here. In summary, we want to emphasize that the upper bound on A does not contradict the ill-posedness of the general tomography problem.
Remark. Equation (4.31) implies that Aϕ can be seen as the minimum-coefficient solution of Aφ = ϕ. Since the sets {w jk,lg } jk∈Z,g=1,...,G do not form tight frames over L 2 (Ω l ), Aϕ is not necessarily also the minimum-norm solution of the equation.
Denoting again by (σ jk,n , u jk,n , v jk,n ) r jk n=1 the singular value decompositions of the matrices A jk , we arrive at the same decomposition of A as in (4.21), but obviously with the new singular values and functions. All results derived above hold analogously for this new decomposition, with the difference that, due to the tightness of the frame, Aϕ is then not only the minimum-coefficient solution of Aφ = ϕ, but also the minimum-norm solution. Furthermore, one can again derive the singular values of A jk explicitly, namely However, note that using this approach, the number of subdomains Ω l,m can increase strongly with the number of guide stars G. Furthermore, note that the number of frames per layer considered in the two approaches presented in this section are directly related to G and M l , respectively. Hence, since the used frames have discontinuities, using the approach based on the frames {w jk,lgm } potentially introduces more discontinuities into the solution in a numerical implementation than the approach based on the frames {w jk,lg }, which is not necessarily desirable in practice.
Remark. In some situations, it is desirable to, instead of the standard L 2 inner product (2.2), work with the weighted inner product where {γ l } L l=1 denotes a normalized, nonzero sequence of weights. Although we do not focus on this issue in our current investigation, it should be noted that adapting the frame decomposition to include this weighted inner product should be straightforward.

Numerical Illustrations
In this section, we present some numerical illustrations showcasing the building blocks of the frame decomposition of the atmospheric tomography operator derived above, namely the frame functions w jk,lg and their corresponding dual frame functionsw jk,lg . While the functions w jk,lg have been defined explicitly and thus can be easily visualized, this is not the case for the functionsw jk,lg , which are only implicitly defined. However, since they are the foundation of the frame decomposition, it is of interest and importance to get some idea about their visual representation, which we aim to do in this section.
It was already noted above that the dual frame functionsw jk,lg can be numerically approximated via formula (3.8), and that this iterative process can be implemented efficiently using the Fourier transform. However, the functionsw jk,lg are problemadapted, i.e., they are dependent on the concrete parameter setting of the considered atmospheric tomography problem, which we thus first need to decide upon. Hence, for our purposes, we take the same setup as in [22] , which in turn is adapted from the MAORY setup [5,12] as planned for the ELT currently being built by ESO. More precisely, we consider a telescope with a diameter of 42 m (the originally planned ELT size), which leads to a circular aperture domain Ω A , as well as 6 NGS positioned uniformly on a circle with a diameter of 1 arcmin, with the corresponding directions α g given in the following table (ρ = 0.000290888, a = ρ · 0.5, b = ρ · 0.866025): For defining the layer profile we choose the ESO Standard Atmosphere, which consists of 9 layers located at the heights h l given in the following table: For the numerical implementation, and in particular for the approximation of the dual frame functionsw jk,lg via formula (3.8), which in our present case reads as jk,lg , w pq,lr L 2 (Ω l ) w pq,lr , (5.1) one can only consider a finite number of indices j, k (and thus also of p, q). Hence, for our numerical illustrations, we restricted ourselves to j, k, p, q ∈ {−63, . . . , 63}, which is motivated by the fact that the wavefront sensors of MAORY cannot reconstruct higher frequencies which might be present in the wavefronts anyway. For computing the dual frame functionsw jk,lg depicted below, we implemented the iterative procedure (5.1) in Matlab (R2019a) on a desktop computer running on a Mac OS with an 8 core processor (Intel Xeon W-2140B CPU@3.20GHz) and 32GB RAM. The iteration itself was stopped after 500 iterations for each of the functions, which is more than sufficient as suggested by the error estimate (3.7) and the fact that already after much fewer iterations the magnitudes of the updates become almost negligible. The computation time for a single dual frame function was more than two hours, although this can be significantly reduced to a couple of seconds by more sophisticated implementations (currently under development). Note again that the computation of the dual frame functions has to be done only once for each AO setting and can be carried out in advance of the actual atmospheric tomography reconstruction. The resulting frame and dual frame functions w jk,lg andw jk,lg , for different values of j, k, l, g are depicted in Figure 5.1 and Figure 5.2, which show absolute value plots of the real part of those functions in linear scale and (for the dual frame functions) also in logarithmic scale. While the frame functions depicted in Figure 5.1 are living on the same layer (l = 2) and differ only in the guide star direction (g = 2 and g = 5), the frame functions depicted in Figure 5.2 vary also in the other parameters. In particular, they are living on atmospheric layers far apart from each other (l = 2 and l = 9). One can clearly (especially from the log-plots) see that even though by definition the frame functions w jk,lg have local support, namely exactly the domains Ω A (α g h l ), the resulting dual frame functionsw jk,lg spread over the whole domain Ω l , which is due to the fact that all the frame functions w jk,lg , also for all other values of g, contribute to their definition. The influence of all the different guide star directions on each of the dual frame functionsw jk,lg is also apparent in the appearance of ring-like structures corresponding to the the domains Ω A (α g h l ). Furthermore, the nicely visible periodicity of the frame functions w jk,lg due to their definition involving an exponential function manifests itself in a periodic pattern of the dual frame functionsw jk,lg .

Conclusion and Outlook
In this paper, we considered the problem of atmospheric tomography by analysing the underlying tomography operator, and derived a frame decomposition for different problem settings. In particular, in contrast to [22] , we considered settings with a mixture of both natural and laser guide stars, and did not place any restrictions on the shape of the aperture shape. The resulting decomposition yields information on the behaviour of the operator and provides a stable reconstruction algorithm for obtaining a minimumcoefficient least-squares solution of the atmospheric tomography problem. Furthermore, we presented numerical illustrations showcasing some of the building blocks of the frame decomposition. In a forthcoming publication currently under preparation, we plan to  and dual frame functionsw jk,lg (middle, right) for j = 3, k = 2, l = 2, g = 2 (top) and j = 10, k = 4, l = 9, g = 4 (bottom). Plotted is the absolute value of the real part in linear (left, middle) and logarithmic scale (right). present numerical reconstruction results from both simulated and experimental AO data. For this, we are in particular developing efficient, accurate, and stable numerical methods for the computation of the dual frame functions.

Support
S. Hubmer and R. Ramlau were (partly) funded by the Austrian Science Fund (FWF): F6805-N36. The authors would like to thank Dr. Stefan Kindermann for valuable discussions on some theoretical questions which arose during the writing of this manuscript.