11institutetext: Max-Planck Institut für Sonnensystemforschung, Justus-von-
Liebig-Weg 3, 37077 Göttingen, Germany 22institutetext: Georg-August-
Universität, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany 33institutetext:
School of Space Research, Kyung Hee University, Yongin, Gyenoggi-Do 446-701,
Republic of Korea

# Correction of atmospheric stray light in restored slit spectra

S. Saranathan 1122 M. van Noort 11 S. K. Solanki 1133

###### Abstract

Context. A long-standing issue in solar ground-based observations has been the
contamination of data due to stray light, which is particularly relevant in
inversions of spectropolarimetric data.

Aims. We aim to build on a statistical method of correcting stray-light
contamination due to residual high-order aberrations and apply it to ground-
based slit spectra.

Methods. The observations were obtained at the Swedish Solar Telescope, and
restored using the multi-frame blind deconvolution restoration procedure.
Using the statistical properties of seeing, we created artificially degraded
synthetic images generated from magneto-hydrodynamic simulations. We then
compared the synthetic data with the observations to derive estimates of the
amount of the residual stray light in the observations. In the final step, the
slit spectra were deconvolved with a stray-light point spread function to
remove the residual stray light from the observations.

Results. The RMS granulation contrasts of the deconvolved spectra were found
to increase to approximately 12.5%, from 9%. Spectral lines, on average, were
found to become deeper in the granules and shallower in the inter-granular
lanes, indicating systematic changes to gradients in temperature. The
deconvolution was also found to increase the redshifts and blueshifts of
spectral lines, suggesting that the velocities of granulation in the solar
photosphere are higher than had previously been observed.

###### Key Words.:

techniques: image processing – instrumentation: spectrographs –
instrumentation: adaptive optics – atmospheric effects

##  1 Introduction

The degradation of observations due to stray light has been a long-standing
issue in solar physics. The term stray light, in recent literature, has come
to refer to extraneous undesirable contributions to pixel intensities in
images and spectra caused by one or more independent physical processes. Stray
light is currently understood to be the sum of multiple components with
possibly differentiable characteristics (see Martínez Pillet (1992); Beck et
al. (2011) for a comprehensive review). For ground-based observations, these
include the refraction of light due to seeing, as well as large-scale
scattering due to the presence of aerosols, dust, molecules, etc. The former
effect varies rapidly with time and is the dominant effect over distances
typically up to a few arcseconds, while the latter exhibits slower modulation
and brings contributions from sources that are much farther away. For example,
the inexplicable brightness of sunspot umbrae (Zwaan, 1965) was found to be
due to atmospheric scattering of light from the full solar disk. Stray light
is also generated within telescopes and instruments by the diffraction of
light, as well as scattering due to the surface roughness of optical elements.
It is even present in space-based observations and may be anisotropic and
inhomogeneous over the field of view (FOV).

Stray-light contamination manifests itself differently in images and spectra.
Images suffer a loss of spatial resolution, as well as a reduction in
intensity contrasts - an effect verified by observations obtained through both
space-based and ground-based instruments (see Uitenbroek et al., 2007; Sánchez
Cuberes et al., 2000; Danilovic et al., 2008). While this contrast reduction
has been satisfactorily modelled in space-based observations (e.g. Mathew et
al., 2009; Wedemeyer-Böhm & Rouppe van der Voort, 2009), a discrepancy still
exists in ground-based observations (Scharmer et al., 2010). The effects of
stray light on spectra are revealed by inversions of spectropolarimetric data.
The mixing of intensities alters the depth, broadening, and symmetry of
spectral lines, introducing different errors in different spectral lines. The
quality of the fit in inversions is therefore generally poorer, and the
reliability of the retrieved physical parameters is reduced. Several
prescriptions to tackle spectral contamination can be found in the literature.

Orozco-Suárez et al. (2007) first proposed such a method by constructing a
local stray-light profile from a 1″region surrounding every pixel. The
inversion code was then allowed to fit the weight of this profile in the pixel
being inverted. Although a significant step forwards, the method is formally
incorrect as the stray-light contribution in a pixel cannot be directly
estimated from the contaminated data. An improvement over this method through
a self-consistent coupled inversion scheme was proposed by van Noort (2012).
However, the method requires prior knowledge of the transfer functions that do
not vary within the FOV. These constraints are likely not satisfied by ground-
based observations (especially from spectrographs), and the preferred
approach, therefore, is to restore the data as perfectly as possible before
inversions.

Restoration of degraded observations is only possible with an estimate of the
transfer function of the optical system. Systematic efforts have been
undertaken along these lines. The earliest such attempts were made by Wanders
(1934), Stumpff (1961), and David & Elste (1962). They defined an analytical
point spread function (PSF) that quantified the extraneous contribution to
intensity in a given pixel, as a function of distance. Subsequent theoretical
and empirical refinements to the model were made by Zwaan (1965), Staveland
(1970), Mattig (1971), and later Martínez Pillet (1992). Although instructive,
early models of the spread function were limited in their ability to improve
the quality of ground-based observations due to the use of coarse assumptions.
A significant boost came from parallel developments in data acquisition,
adaptive optics (AO), and image restoration techniques, such as speckle
(Keller & von der Luehe, 1992), phase-diversity (Paxman et al., 1992; Löfdahl
& Scharmer, 1994), and multi-frame blind deconvolution (MFBD; van Noort et
al.2005), which attempt to model the true PSF as it changes with time. In the
current paradigm, AO systems and image restoration methods are routinely
employed to circumvent degradation and/or reverse it to an acceptable level.

Despite the progress, persistent signs of degradation in ground-based
observations have prompted a careful search for the sources of residual stray
light. Scharmer et al. (2010) found that high-order wave-front aberrations -
corresponding to the wings of the PSF, which were systematically
underestimated by MFBD - could be partly responsible for residual stray light
in ground-based observations. They used a statistical method based on
Kolmogorov statistics and knowledge of the Fried parameter $r_{0}$ (Fried,
1966) to approximate corrections to the transfer functions returned by MFBD
and subsequently reduce the amount of residual stray light in observations. We
propose a refinement to this method by using the values of the Fried parameter
to generate fully turbulence-degraded synthetic observations. Restoring the
synthetic observations along with the actual observations in parallel, using
MFBD, allows us to self-consistently quantify the efficiency with which the AO
reversed seeing-induced degradation in the actual observations. Given that AO
systems are not perfect, these efficiencies can be used to arrive at estimates
of the residual degradation present in the AO-corrected, MFBD-restored
observations. In this paper, we discuss the theory of the method and its
implementation on spectropolarimetric data obtained from the Swedish Solar
Telescope (SST, Scharmer et al., 2003).

##  2 Residual degradation in ground-based observations

When observing astrophysical objects through a telescope of finite aperture,
diffraction sets a fundamental limit to the spatial resolution of the final
image(s). In the presence of atmospheric seeing, or due to the instruments
themselves, additional degradation may be introduced. For a given optical
system, quantification of image degradation in the observations is made with
its transfer function. If the optical system includes the atmosphere, such as
is the case for ground-based observations, the transfer functions are also
time-dependent. We write the image formation model in the fourier domain:

$\hat{I}(\@vec{k},t)=\hat{T}_{atm}(\@vec{k},t)\hat{T}_{inst}(\@vec{k})\hat{O}(\@vec{k})+\hat{N}(\@vec{k}),$
(1)

where $\@vec{k}$ and $t$ denote spatial frequency and time, respectively,
$\hat{T}_{atm}(\@vec{k},t)$ is the transfer function of the atmosphere,
$\hat{T}_{inst}(\@vec{k})$ is the transfer function of the telescope and
relevant instruments, and $\hat{O}(\@vec{k})$ represents the true solar scene.
$\hat{N}(\@vec{k})$ represents additive noise, and $\hat{I}(\@vec{k},t)$ is
the flat-fielded, dark-current subtracted image, recorded by the detector. In
ground-based observations, the explicit dependence of $\hat{O}(\@vec{k})$ on
time is dropped because the timescales of solar evolution are much larger than
the fluctuation timescales of the Earth’s atmosphere.

To model $\hat{T}_{atm}(\@vec{k},t)$, the approximation that it is the result
of a pure phase error in the wave front is commonly made. Wave-front phase
errors are introduced either by the optical elements themselves (static) or by
the atmosphere (dynamic). The $\hat{T}_{atm}(\@vec{k},t)$ is then computed
using:

$\displaystyle\zeta(\@vec{x},t)$ $\displaystyle=\exp{(i\phi(\@vec{x},t))},$
(2) $\displaystyle\hat{T}_{atm}(\@vec{k},t)$
$\displaystyle=R(\zeta,\zeta^{*}),$ (3)

where $\phi(\@vec{x},t)$, $\zeta(\@vec{x},t)$, and $R(\zeta,\zeta^{*})$ are
the wave fronts, the pupil function, and the auto-correlation function,
respectively. A perfectly flat wave front corresponds to an ideal transfer
function that drops nearly linearly with $\@vec{k}$. In the presence of
aberrations, the decline with $\@vec{k}$ is much steeper, forming images with
a reduced signal to noise ratio (S/N) at high spatial frequencies. To prevent
this, an AO system can be deployed combining a wave-front sensor (WFS) with a
deformable mirror (DM) to measure and compensate wave-front aberrations in
real time. The efficacy of such systems is a function of the number of
actuators driving the DM, sampling of the wave front by the WFS, and the
bandwidth of the AO. A review on the design and performance of AO systems is
found in Rimmele & Marino (2011).

Aberrations are often quantified by projecting the corresponding spatial
distribution of the wave-front phase error onto an ordered set of basis modes
(Zernike polynomials, Karhunen-Loève functions, etc.):

$\phi(\@vec{x},t)=\sum_{l=0}^{\infty}{C_{l}\leavevmode\nobreak\
Z_{l}(\@vec{x},t)},$ (4)

where $l$ denotes the mode number, $\phi$ is the wave-front phase error,
$Z_{l}$ are the basis modes, and $C_{l}$ are mode amplitudes. In closed-loop
operation, AO systems act as high-pass filters (Rimmele & Marino, 2011).
Specifically, the amplitudes of low-order aberrations are suppressed better
than the high orders. The meaning of ”low” and ”high” in this context is
somewhat arbitrary as no such sharp boundary exists. Nonetheless, it is
instructive to rely on such a classification for the purposes of motivating
what follows. The action of AO on the wave-front error can then be written as:

$\phi^{AO}=\sum_{l=0}^{\infty}(1-\epsilon^{AO}_{l})\leavevmode\nobreak\
C_{l}\leavevmode\nobreak\ Z_{l},$ (5)

where explicit dependence on $\@vec{x}$ and $t$ is dropped for brevity.
$\phi^{AO}$ is the phase-error after compensation by the DM, and
$\epsilon^{AO}_{l}$ quantifies the suppression of mode $l$. That is,
$\epsilon^{AO}_{l}$ characterises the efficiency of the AO in compensating for
phase errors with a spatial distribution according to mode $l$. For perfect
compensation, $\epsilon^{AO}_{l}=1$, and the residual wave-front error would
be identically zero over the entire aperture. Realistically,
$\epsilon^{AO}_{l}$ would depend on $l$ and typically vary between zero and
one. Given that AO systems act as high-pass filters, $\epsilon^{AO}_{l}$ would
be close to one at the lower orders and close to zero at the higher orders.
One would also expect the AO to become completely ineffective beyond a certain
point, that is:

$\epsilon^{AO}_{l}=0,\hskip 10.0pt\forall\hskip 10.0ptl>L.$ (6)

However, during the application of our method to actual observations, several
scenarios were observed where $\epsilon^{AO}_{l}$ was non-zero, even when
$l>200$. Further interpretations of this characteristic are provided in Sec.
(7) of this paper. Additional complications concerning the nature of
$\epsilon^{AO}_{l}$ are that they may potentially vary with seeing, and within
the FOV. The dependence on seeing is due, in part, to the way in which the WFS
detects the shape of incident wave fronts. For instance, a Shack-Hartmann WFS,
such as the one used at SST, simultaneously images a local region in the FOV
through an array of sub-apertures. Cross-correlations of the sub-images are
then used to compute their relative shifts, and subsequently the local tip-
tilt of wave fronts over the sub-apertures. The reliability of these shifts is
particularly impacted when the seeing is very poor. Correspondingly, the
$\epsilon^{AO}_{l}$ would be low. The second complication - dependence of
$\epsilon^{AO}_{l}$ on the location - is due to anisoplanatism. When using an
AO with a single mirror and a WFS, only the seeing that is common to the whole
FOV can be corrected over the whole FOV. The more distant the atmospheric
turbulence is, the smaller the area of the FOV that shares the same wave-front
aberrations (isoplanatic patch), and the smaller, therefore, the area that can
be fully corrected. However, any corrections that are specific to one
isoplanatic patch are unavoidably applied to the whole FOV, and therefore
added to the aberrations everywhere else, thus causing a statistical increase
in the wave-front errors there, rather than a decrease. To avoid this, a
ground-layer WFS can be used that relies on a larger area than one isoplanatic
patch to calculate the wave-front slope, so that the cross-correlations are
less sensitive to wave-front errors that are not shared across that area.

In practice, this method reduces, but does not completely eliminate, the
sensitivity to anisoplanatic wave-front aberrations, so that the correction of
the wave-front phase errors becomes increasingly ineffective with distance
from the region to which the AO is locked onto (hotspot). This leads to
statistically poorer corrections at the edges of the FOV. A combined effort of
evaluating how the corrections get worse with distance, and how the high
altitude seeing intrinsically varies within the FOV is required in order to
characterise this effect. In the context of this section, the most important
realization concerning $\epsilon^{AO}_{l}$ is that they are generally unknown.
This leaves the profile of residual wave-front phase errors unknown, and
therefore the corresponding residual transfer functions unknown. Revising the
image formation model to allow AO-compensations, and denoting the residual
transfer functions corresponding to $\phi^{AO}$ as $\hat{T}^{AO}$, we write:

$\hat{I}(\@vec{k},t)=\hat{T}(\@vec{k},t)\hat{O}(\@vec{k})+\hat{N}(\@vec{k}).$
(7)

It is implicit that $\hat{T}(\@vec{k},t)$ corresponds to the combined optical
system where the time-dependence is solely caused by the atmosphere. The
subscripts have, therefore, been dropped. Further improvements to the quality
of observations are availed with the use of post-facto processing methods such
as MFBD(van Noort et al., 2005; van Noort, 2017) that do not require prior
knowledge of the degrading transfer functions. In MFBD image restoration,
estimates of the transfer functions are made directly from the degraded
observations. This is achieved by expanding the wave-front phase errors onto
an orthogonal basis and fitting for the amplitude of every mode. However, the
basis used by MFBD is ultimately finite and the expansion is truncated to a
mode limit $M$, which causes a systematic underestimation of the wings of the
Point Spread Functions (PSFs). Therefore, while MFBD corrects for degradation
caused by lower-order modes, the effects of the higher orders (stray light)
are left largely untreated. Since image restoration with MFBD is normally
performed on observations that are already AO-compensated, one may assume
that:

$\displaystyle\epsilon^{AO+MFBD}_{l}=1,\hskip 10.0pt\forall\hskip 10.0ptl<M,$
(8)

that is, residual degradation in AO-compensated MFBD-restored observations are
exclusively due to the higher-order modes. A peculiar state of compensation is
realised when the AO corrects more modes than MFBD ($L>M$). In this case,
modes with $l$ ¡ M are perfectly compensated by the AO and MFBD, modes with M
¡ $l$ ¡ L have been compensated by the AO (but with unknown efficiency), and
modes with $l$ ¿ L have not been compensated at all. Alternatively, one may
try to perform image restoration with $M$ increased to ~$L$ so that the issue
of the unknown correction efficiency of the AO is circumvented entirely.
However, given that the efficiency of wave-front corrections made by the AO
vary within the FOV, mitigation of this issue is attempted during image
restoration - by segmenting the FOV into a grid of patches and restoring each
patch individually. In doing so, the assumption that wave-front phase errors
over each patch are isoplanatic is made. This imposes a limitation on the
admissible physical size of the PSFs, and indicates that restoring with
arbitrarily high mode limits is not beneficial. The preferable approach is
therefore to develop a self-consistent method of removing residual stray light
from the restored data directly.

We clarify that the term ’residual stray light’ here refers exclusively to the
component of stray light arising from high-order aberrations. Other sources of
stray light, such as true scattering in the atmosphere and within the
instrument, are beyond the scope of this paper. Correspondingly, we refer to
the PSF describing stray-light contamination as a stray-light PSF.

##  3 Residual degradation in slit spectra

In the previous section, we reviewed the limitations of AO and MFBD and
qualitatively addressed the nature of residual degradation in ground-based
observations. We now describe the characteristics of residual stray-light
contamination in spectropolarimetric data obtained with Spectrographs. A
detailed discussion on how such data are obtained and restored is found in van
Noort (2017). Here, we only briefly review the relevant details. In the
interest of keeping the description concise, we ignore the polarimetric aspect
of such data. We assume that the conclusions drawn in this section apply to
all the Stokes components equally.

Slit-spectra are obtained with a detector that collects light that has passed
through a slit, and that has been dispersed by a dispersive element, usually a
grating. Another detector (slit-jaw camera) is used to image the slit and the
surrounding FOV through a broadband filter spanning the spectral region
covered by the spectra. Both the detectors are operated synchronously so that
the recorded images/spectra can be assumed to be degraded by the same seeing,
and the frame-rates of the detectors are kept high enough to prevent seeing-
induced distortions, and seeing-induced crosstalk. The images recorded by the
slit-jaw camera are restored in bursts with MFBD to yield estimates of the
degrading PSFs as a function of space and time. However, to restore the slit-
spectra, one cannot directly deconvolve them as the intensities perpendicular
to the slit have not been recorded co-temporally. Instead, the problem is
reformulated into a set of linear equations that connect the observed
intensities with the undegraded solar spectra through the PSFs. That is, Eq.1
can be re-written into a matrix equation of the form:

$\hat{M}\textbf{x}=\textbf{y},$ (9)

where x represents the coveted undegraded spectra, y represents the seeing-
degraded intensities recorded by the spectral camera, and $\hat{M}$ is, in
principle, a non-sparse non-square matrix that represents an instantaneous
PSF. The values of such PSFs are obtained from the MFBD restoration of the
slit-jaw images, and therefore the matrix $\hat{M}$ is completely known. Eq. 9
can be constructed for every wavelength point, and for each instant of time,
yielding an over-determined system of coupled equations. Since $\hat{M}$ is
not a square matrix, one could instead solve the pseudo-inverse problem:

$\displaystyle\hat{A}\@vec{x}=\@vec{b},$ (10)

where $\hat{A}=\hat{M^{T}}\hat{M}$ is a square matrix, and
$\@vec{b}=\hat{M^{T}}\@vec{y}$. However, since a large number of equations and
unknowns are typically involved, it is difficult to directly invert the matrix
$\hat{A}$. The solution therefore, is obtained using a variation of the
iterative Lucy-Richardson scheme. Since the matrix $\hat{A}$ is constructed
using PSFs whose wings are underestimated, the solution thus obtained is
likely to still contain some residual stray light. If we were to estimate a
correction $\delta\hat{A}$ to $\hat{A}$ that addressed this issue, the most
self-consistent way forwards would be to construct a new matrix
$\hat{A^{\prime}}=\hat{A}+\delta\hat{A}$ and restore the raw data again.
However, due to the rapid scaling of computational workload with the size of
the domain, using $\hat{A^{\prime}}$ becomes prohibitively expensive. We
therefore propose an alternate method of removing stray light directly from
the restored spectra.

First, the Stokes signals are converted to modulated intensities using an
arbitrary modulation scheme for each wavelength, yielding four modulated
intensities corresponding to each Stokes parameter. This is followed by
degrading the modulated intensities with a diffraction limited PSF
corresponding to the aperture of the telescope. The degraded intensities
corresponding to each Stokes parameter are then deconvolved using the same set
of PSFs that describe residual stray light exclusively across all wavelengths.
Although the deconvolution is bound to amplify the noise in the signals, this
amplification is uniform across the Stokes parameters and does not, therefore,
introduce spurious polarization signals(polarimetric precision is preserved).
After the deconvolution is performed, the same modulation scheme used in the
forward modulation is used to invert the deconvolved intensities back into
Stokes spectra. The method described above requires that a set of suitable
stray-light PSFs be found. Two key aspects of such a set of stray-light PSFs
have to be noted.

Firstly, at any given wavelength, slit-spectra are not truly images because
the intensities along the scan-direction have been recorded at different
times. It then follows that residual stray light generally varies along the
scan-direction as the seeing evolves with time. Therefore, to remove residual
stray light from a single ’image’ at a particular wavelength, multiple stray-
light PSFs spanning the duration of the scan have to be computed. This
complication is not shared by observations made with filtergraphs, as
intensities there are obtained strictly co-temporally. Secondly, consider a
fixed point on the Sun. Let the slit sample this point at time $t_{0}$.
According to Eq. (10), the true solar intensity of this point (solution) is
computed from the intensities in the slit recorded between $t_{0}-\Delta t$
and $t_{0}+\Delta t$, where $2\Delta t$ is the time it takes for the slit to
scan the spatial extent of the PSF. Naturally, intensities recorded at $t_{0}$
carry the highest weights towards the solution while those recorded at
$t_{0}\pm\Delta t$ carry the lowest. The correct amount of residual stray
light is therefore described not just by one stray-light PSF computed at
$t_{0}$, but through a weighted average of many stray-light PSFs computed
between $t_{0}-\Delta t$ and $t_{0}+\Delta t$, with the weights specified in
the matrix $\hat{A}$.

##  4 Estimating the stray-light PSF

Figure 1: Kolmogorov wave fronts simulated for different values of the Fried
parameter $r_{0}$, for an aperture size of 1 m. The mean values of all the
wave fronts have been subtracted, and a uniform grey scale has been used in
the illustration of the wave fronts.

Scharmer et al. (2010) followed a statistical method, where values of the
Fried parameter ($r_{0}$) were used to simulate a set of wave fronts obeying
Kolmogorov statistics. $r_{0}$ values were calculated using a setup employing
a wide-field wave-front sensor installed before the tip-tilt mirrors and the
DM. At the time, the DM was driven by a 37-electrode configuration,
corresponding to roughly 37 KL-modes. However, the DM at SST has since been
upgraded to a configuration using ~85 electrodes, or roughly 100 KL-modes. As
mentioned earlier, image-restoration is typically performed using MFBD with a
mode limit of ~40, so that the amplitude of KL-modes between $l=40$ and
$l=100$ in the residual wave fronts is unknown. We propose that the average
correction efficiency of the AO can be computed from the following procedure:
(i) the $r_{0}$ values are first used to simulate wave fronts obeying
Kolmogorov statistics; (ii) the corresponding transfer functions computed
through Eq. (2) & Eq. (3) are used to degrade an image synthesised from MHD
simulations; (iii) both the sets of observations(synthetic and actual) are
restored using MFBD (which returns the mode amplitudes); (iv) the ratio of the
standard deviations of the mode amplitudes (Eq. 11) in both the restorations
yields the efficiency of the AO as a function of the mode. This efficiency is
mathematically expressed as:

$\displaystyle\epsilon^{AO}_{l}=1-(\sigma_{l}^{obs}/\sigma_{l}^{syn}),$ (11)

where $l$ denotes the mode number, and $\sigma^{syn}_{l}$ and
$\sigma^{obs}_{l}$ are standard deviations of the mode amplitudes returned by
MFBD restoration of the synthetic and original recorded observations,
respectively. Interpretation of Eq. (11) is straightforward for the trivial
cases. For example, if the recorded observations are found to not have any
degradation caused by mode $l$, the standard deviation for mode $l$ computed
from image restoration would ideally be zero. This then implies that the AO
perfectly compensated mode $l$ in the incident wave fronts
($\epsilon^{AO}_{l}=1$).

A potential source of error in the calculation of $\epsilon^{AO}_{l}$ is the
assumption that MFBD works identically on the input observations, irrespective
of the state of their degradation. To address this issue, the procedure
outlined above can be followed iteratively - after the first round of
$\epsilon^{AO}_{l}$ are found, these can be used to subtract modes from the
incident wave fronts. These partially compensated wave fronts can then be used
to generate another set of degraded observations followed by another round of
image-restoration, and the cycle can be followed until the computed
efficiencies are found to converge. Once a converged set of
$\epsilon^{AO}_{l}$ are found, residual stray light in the observations is
estimated by first perfectly subtracting all modes with $l<M$, followed by
subtracting modes with $l>M$ with efficiency $\epsilon^{AO}_{l}$. The residual
wave fronts thus generated are largely dominated by higher-order aberrations,
and they approximate the amount of residual stray light still present in the
observations.

##  5 Kolmogorov phase screens

A frequent approach to characterizing AO-systems is through the use of
Kolmogorov phase screens, which are simulations of the spatial distribution of
phase errors (wave fronts) incident on the telescope. These simulations are
statistical in nature since modelling the exact propagation of a wave front
requires knowledge of the refractive index at all spatial locations within the
path of the wave front. The statistical properties of these Kolmogorov phase
screens are well established in the works of Noll (1976), McGlamery (1976),
Lane et al. (1992), and Glindemann et al. (1993). They conform to the
Kolmogorov-Obukhov power law, and can be numerically simulated (Nagaraju &
Feller, 2012) as a function of $r_{0}$:

$\displaystyle\Gamma(\@vec{k})=0.023(D/r_{0})^{-5/3}|\@vec{k}|^{-\beta},$ (12)

where $\Gamma(\@vec{k})$ is the power at spatial frequency $\@vec{k}$, $D$ is
the diameter of the phase screen, and $\beta=11/3$. Fig. 1 shows examples of
such phase screens simulated over a circular aperture of $1$ m diameter. The
mean value of all wave fronts has been subtracted, but no other modes have
been removed. The phase screens are then used to compute the PSFs with which
the synthetic MHD-generated images are degraded, and an appropriate amount of
noise, by comparison with slit-jaw observations for example, is added to the
synthetic data. We make a note that the validity of the MHD simulations is
irrelevant as we are only interested in the synthetic images. Restoring the
synthetic data and the original recorded observations with MFBD facilitates a
comparison of their respective mode amplitudes.

##  6 Observations and data reduction

Figure 2: Panels showing the FOV and its location with respect to the full
solar disk. Left: Continuum intensity image from the restored observations.
The red arrow points towards the disk centre. Right: Line-of-sight magnetogram
from HMI taken at approximately the same time. The red box highlights the FOV
of the observations.  Figure 3:  Top to bottom: Stokes $Q$, $U$, and $V$
profiles averaged over the FOV. The dotted lines and solid lines show the
spectra before and after crosstalk correction, respectively. The grey shaded
areas indicate the wavelengths over which the dot products are computed. The
units of intensity are arbitrary, and the profiles are not normalised with
respect to Stokes $I$.

The dataset used in this paper was obtained on 28 June 2016 at the SST
(Scharmer et al., 2003), using the TRIPPEL spectrograph, a detailed
description of which has been provided by Kiselman et al. (2011). A fast
polarimeter designed for slit-spectra - based on the concept of FSP (Iglesias
et al., 2016) \- was installed and used during the same campaign (TRIPPEL-SP,
Doerr et. al 2019, in prep). The scan was made under moderate seeing
conditions with the help of an AO consisting of an 85-electrode bi-morph
mirror, a Shack-Hartmann wave-front-Sensor (WFS), a correlation tracker, and a
separate mirror for tip-tilt adjustments. The AO reported an average Fried
parameter value of $13$ cm over the duration of the scan.

Two Jai SP20000 CMOS cameras were used for recording the slit-jaw images,
whereas the spectra were recorded using two Jai SP12000 CMOS cameras. The
latter allowed 4-Mpixel frames to be recorded at a rate of up to $400$ Hz, the
designed modulation rate of the FSP modulator. All cameras were synchronised
using an external trigger signal. Although frames were recorded with a spatial
sampling of $0.044$″, the sampling of the restored slit-spectra was scaled
down to $0.06$″- corresponding to the critical sampling of the SST at $6300$
Å.

###  6.1 Image restoration

The slit-jaw frames were restored in bursts of 2000 frames (5 s) with a
temporal overlap of 3.75 s. As indicated in earlier sections, the restoration
was performed over a grid of patches that were 128 $\times$ 128 pixel2 in
size, with an overlap of $118\times 118$ pixel2. The PSFs returned by MFBD for
each patch were interpolated across the length of the slit. Karhunen-Loève
modes were used as the basis functions, and the mode limit for the
restorations was set to $M=45$. As outlined in Sec. (3), the values of the
PSFs were used to construct the matrix $\hat{A}$, and subsequently restore the
data recorded by the spectral camera.

The restored slit-spectra cover a 53″$\times$ 53″FOV around a mostly uni-polar
active region Plage, which is located towards the east limb at a heliocentric
angle of 23 degrees ($\mu=0.92$). The left panel in Fig. 2 shows a Stokes I
image from the restored dataset, extracted at a continuum wavelength near
$6296.0$ Å. The spectral window of the observations ranges from 6288.0 Å to
$6304.0$ Å, with a sampling of 8.6 mÅ. The wavelengths were calibrated by
fitting the flat-field (obtained at disk centre) using a least-squares method
to the standard FTS spectrum from Delbouille et al. (1973). A total of 16
usable solar spectral lines were identified in the spectral window, including
the iron pair at $6301.5$ Å and $6302.5$ Å. On the panel in the right of Fig.
2, a LOS magnetogram from HMI (Scherrer et al., 2012) is displayed for
context.

###  6.2 Crosstalk and grey stray light

Figure 4: Panels showing degraded images synthesised from an MHD simulation at
a continuum wavelength near $6300\AA{}$ for different values of $r_{0}$
(mentioned in the panels). The grey scales of the degraded images are
identical. The grey scale range of the undegraded image has been enlarged to
avoid extreme contrasts. The sampling in all images is $0.044\arcsec$. Top
left: Undegraded image (no added noise) corresponding to $r_{0}$ = $\infty$.

The restored images indicated the presence of crosstalk from
$I\xrightarrow{}Q,U,V$ and from $V\xrightarrow{}Q,U$, a consequence of
neglecting the polarizing effects of the telescope optics so far. Although
models of the Mueller matrix of SST exist, we found it convenient to apply ad
hoc corrections to the spectra. The crosstalk from $I\xrightarrow{}Q,U,V$ was
corrected by subtracting the average offsets (averaged over the FOV at
continuum wavelengths) of Stokes $Q,U$, and $V$ (normalised with respect to
Stokes $I$). To remove crosstalk from Stokes $V$, we parameterised the Stokes
vectors on the Poincaré sphere, and rotated the Stokes components such that
the FOV-averaged Stokes $Q,U$, and $V$ were symmetric, symmetric, and anti-
symmetric in wavelength with respect to the line-centres, respectively. In
doing so, we assume that the core of the plage makes the dominant contribution
to the Stokes $V$ spectra, and that in these regions, gradients in the line-
of-sight velocities and magnetic fields are not too strong to produce
significant intrinsic asymmetry. The shaded region in Fig. 3 indicate
weighting of the profiles to include only solar lines in the
(anti)symmetrization procedure. Since most of the crosstalk is from
$V\xrightarrow{}Q,U$, the linear polarization signals look significantly
different after crosstalk removal, while Stokes $V$ remains largely unchanged
(Fig. 3).

To remove the wavelength-independent component of the scattered light possibly
caused by the extended wings of the spectral PSF (spectral veil, (Borrero, J.
M. et al., 2016)), we make the assumption that the spectral PSF of the
instrument is approximated to first order by a gaussian kernel with an initial
guess value for its Full-Width Half Maxima (FWHM). The flat-field spectra
obtained at disk centre is then expressed as the sum of the standard FTS
spectra convolved with this gaussian kernel, and the additive spectral veil
(expressed as a fraction of the continuum intensity). We then retrieved, both
the FWHM of the gaussian kernel and the additive component, by iteratively
fitting the flat-field spectra using a $\chi^{2}$ minimization method. The
FWHM was found to be 41mÅ and the fraction of spectral veil was found to be 5%
of the continuum intensity.

###  6.3 Fried parameter

Figure 5: Plots of quantities obtained from MFBD restorations. Left: Standard
deviation of mode amplitudes plotted as a function of mode number. Mode
amplitudes were obtained from MFBD restoration of degraded images. Only up to
100 modes are shown here. Right: Average efficiency of the AO obtained from
Eq. (11) plotted for different different values of $r_{0}$ Figure 6:  Stray-
light PSFs for different $r_{0}$, computed from the residual phase screens.
The Strehl ratios of the PSFs are indicated in top left of each panel. At
$r_{0}$ = 20 cm the radius encircling 90% energy is around 1.1″.

In closed-loop operation of the AO at the SST, an $r_{0}$ value for every
second is computed from the voltages output to the DM-electrodes, as well as
the variance of residual shifts of the $24\times 24$ pixel2 sub-images in four
select sub-apertures, similar to a Differential Image Motion Monitor (Sarazin
& Roddier, 1990). Further, another concurrent set of $r_{0}$ values are also
reported using the same method, but with smaller $8\times 8$ pixel2 sub-
images. Since cross-correlations between larger images are more sensitive to
shifts that are global, the values of $r_{0}$ computed from the $24\times 24$
pixel2 sub-images are biased more towards ground-layer seeing that is shared
by the full FOV. To avoid this bias, we use the values of $r_{0}$ computed
from the $8\times 8$ pixel2 sub-images instead, and consider these to be more
representative of seeing originating in both the ground and upper layers.

###  6.4 Calculating the stray-light PSF

An ensemble of $2000$ Kolmogorov phase screens were created for each value of
$r_{0}$ between $r_{0}=6$ cm and $r_{0}=21$ cm in steps of $1$ cm. A snapshot
from an MHD simulation of a quiet Sun, performed with MURaM (Vögler et al.,
2005), was then used to synthesise an image at a continuum wavelength near
$6300$Å with the SPINOR radiative transfer code (Frutiger et al., 2000). The
synthesised image was subsequently degraded with the transfer functions
corresponding to each of the $2000$ Kolmogorov phase screens simulated for
each $r_{0}$, followed by the addition of Poisson noise to each of the
degraded images. The amount of added noise was calibrated from the power-
spectra of the slit-jaw images in the recorded observations. Fig. 4 shows a
panel of these degraded images for the binned values of $r_{0}$.

We then setup parallel image restorations of the synthetic degraded data
generated above, and the recorded observations. The mode limit used in the
restorations was set to $M=250$. Since MFBD is applied on AO-compensated data
in practice, we avoided restoring the synthetic data for $r_{0}<10$ cm.
Further, due to the increased mode limit, the restoration of the recorded
observations was constrained to the patch at the centre of the FOV. The
implications of this restriction are that the stray-light PSFs that we compute
would likely underestimate the residual stray-light contamination in the edges
of the FOV.

MFBD returns the fitted mode amplitudes for each KL-mode, and for each
degraded image. This allowed us to compute their standard deviations, which
are plotted in the left panel of Fig. 5 for $r_{0}=5,10,20$ cm. The dotted and
solid standard-deviation curves correspond to the synthetic degraded data, and
the original recorded observations, respectively. Only the first $M=100$ modes
are shown to magnify the behaviour of these curves at the lower orders, but
they continue to decrease all the way upto mode $M=250$. It is immediately
seen that the standard deviations of the mode amplitudes in the recorded
observations are lower than those corresponding to synthetic data, as it
should be since the recorded observations are AO-compensated while the
synthetic data are not AO-compensated. Another point to note from these curves
is that the standard deviations of the synthetic degraded data decrease more
rapidly than the recorded observations. This suggests that the AO-system at
the SST works to keep the level of residual degradation in the observations
fairly stable, irrespective of the seeing conditions.

From the standard-deviation curves, we computed the efficiency of the AO-
system as described by Eq. (11), for every second of the scan. The panel on
the right of Fig. 5 shows the computed efficiency of the AO
($\epsilon^{AO}_{l}$) as a function of the basis mode, averaged over the
realizations for $r_{0}=5,10,15$ cm. We now discuss the trends and features of
these efficiency curves.

The first is the presence of a periodic saw-tooth pattern, observable also in
the standard-deviation curves although not as strongly. These peaks correspond
to the radial KL-modes, and may be caused by the ordering of the KL-basis
modes. Since the radial KL-modes correspond to fluctuations at a lower spatial
frequency compared to the preceding and succeeding azimuthal modes, they would
always have more power according to the Kolmogorov-Obukhov power law. However,
these peaks should partly cancel each other out in the efficiency curves as
efficiency is obtained as a ratio of standard deviations. As to why they are
very prominent is not fully clear. A possible cause may be that the MFBD
restoration is somewhat sensitive to the level of residual degradation in the
data. An approach to circumvent this issue would be to use an iterative
scheme, where AO-efficiencies at a given iteration are used to simulate AO-
compensation on the fully uncompensated phase screens, and thereby generate
another dataset. Upon restoring this dataset with MFBD again, and comparing
the resulting standard deviations with the observations, one can arrive at a
correction to the AO-efficiencies. The procedure can be repeated iteratively
until the AO-efficiencies are found to converge. In our application, we noted
that AO-efficiencies were usually overestimated at the very first iteration,
and therefore the stray-light PSFs corresponding to this iteration likely
underestimated the amount of residual degradation. However, upon subsequently
deconvolving the data, we noted the amplification of artefacts in the spectra.
The more aggressive the stray-light PSF, the larger the artefacts were
amplified beyond the noise level. The properties of these artefacts are
further discussed in Sec. (7). In order to avoid these artefacts, we therefore
considered the stray-light PSFs obtained at the very first iteration to be
final. We note that owing to this choice, we may slightly underestimate the
amount of residual stray light in our slit spectra.

The second is the indication that the AO continues to compensate very high-
order modes consistently, particularly when the seeing is poor. While it is
possible that certain higher-order modes are able to be compensated by the AO,
the indication that modes upto $L=250$ are corrected with an efficiency
averaging $30\%$ when $r_{0}=5$ cm could also be due to theoretical errors in
the variances of mode amplitudes in the simulated phase screens. For example,
modifications to the exponent $\beta$ in Eq. (12) have been suggested in
literature (Nicholls et al., 1995; Rao et al., 2000) with $2<\beta<4$,
implying that the stray-light correction that we perform could be both an
underestimation ($\beta=2$) or an overestimation ($\beta=4$). Another factor
that affects the standard-deviation curves are errors in $r_{0}$. Ideally, for
a completely unbiased estimate of $r_{0}$, the size of the sub-images in the
WFS should be smaller, still, than $8\times 8$ pixel2. However, a lower bound
exists as the AO cannot lock onto regions without sufficient intensity
contrast. Consequently, the $r_{0}$ that is computed from the cross-
correlation of these sub-images may deviate from the true value at any given
time. The third feature of the AO-efficiency curves is that they also appear
to take on negative values. Going back to the standard-deviation curves, this
implies that the recorded observations have more degradation in the higher-
orders than that caused by completely uncompensated Kolmogorov wave fronts!
One possible explanation for this effect is that the AO-mirror adds spurious
high-order aberrations to the incident wave fronts. This may be because the
mirror cannot exactly reproduce the Zernike or KL-basis modes, leading to the
amplification of high-order modes in the residual wave fronts.

Despite some of the surprising results, the method outlines a self-consistent
way of computing the performance of the AO-system from the recorded
observations themselves. Since the slit spectra are AO-compensated and have
been restored with MFBD, the amount of residual degradation in the slit
spectra were estimated in the following way - first, modes until $M=45$ were
subtracted completely from the fully uncompensated phase screens. This step
simulates the joint correction of AO and MFBD to the first $45$ KL-modes,
which is assumed to be perfect. In the next step, modes from $M=45$ to $M=250$
were removed using the AO-efficiencies:

$\displaystyle\phi_{R}=\phi_{S}-\sum_{l=1}^{45}{C_{l}\leavevmode\nobreak\
Z_{l}}-\sum_{l=46}^{250}{\epsilon^{AO}_{l}\leavevmode\nobreak\
C_{l}\leavevmode\nobreak\ Z_{l}}.$ (13)

Here, $\phi_{R}$ is the residual wave front, $\phi_{S}$ is the simulated wave
front, and $C_{i}$ is the amplitude of mode $l$. The PSFs corresponding to the
residual wave fronts were averaged over the realizations to yield a smooth
stray-light PSF for every second of the scan. Following this, the stray-light
PSFs were also running-averaged in time according to the weights specified in
the matrix $\hat{A}$, and reduced to one stray-light PSF per row of the FOV
(row is the scan-direction). These stray-light PSFs are displayed in Fig. 6
for different values of $r_{0}$ in log scale. The Strehl ratios of the PSFs
are indicated in the top left of each panel, and vary from 50% in poor seeing
to approximately 80% in good seeing.

##  7 Results and analysis

Figure 7:  Comparison of the slit spectra at a continuum wavelength near 6300
Å, before and after stray-light deconvolution. The grey scale has been kept
uniform across both images, and the units of intensity shown here are
arbitrary. The images shown here are at the centre of the larger FOV.

Following the procedure outlined in Sec. (3), the Stokes profiles were first
converted to modulated intensities using a balanced modulation scheme. The
four modulated intensities were then subsequently deconvolved with an airy-
disk PSF corresponding to the diffraction limit of the aperture of the SST, at
each wavelength. To perform the deconvolution, the Lucy-Richardson scheme was
used. The deconvolution was non-regularised and terminated after $20$
iterations as tests with synthetic data indicated that convergence is achieved
quickly. We now discuss the effects of the stray-light removal procedure.

###  7.1 Effects on RMS granulation contrast

After deconvolution with the stray-light PSF, RMS granulation contrasts were
found to increase to approximately to $12.3\%$ in the quiet Sun, from $9.0\%$
for the MFBD restored data. In Fig. 7, two continuum intensity images of the
centre of the FOV are shown before (left) and after (right) the deconvolution.
The grey scale has been kept equal so that the increase in contrast is readily
visible in the deconvolved spectra. Since the method we followed works best at
the centre of the FOV, the increase in contrast at the edges of the FOV is
lower than it theoretically should be.

In order to understand the reduction of stray light in our data, a comparison
with synthetic images is required. We therefore synthesised an image at a
continuum wavelength near $6300$ Å (same region as our spectra) from an MHD
simulation, and degraded it to $86.5\%$ of the diffraction limit of the
telescope. This percentage was extracted from the power spectrum of
deconvolved observations, and locating the spatial frequency above which noise
dominates. Comparing the RMS contrasts, we found that a deficiency of
$2.5-3\%$ still exists in our data. Since we are exclusively removing stray
light due to higher-order aberrations, contributions from instrumental
sources, and from atmospheric scattering still exist. However, as to how much
they contribute to an increase in RMS contrast is still unclear.

###  7.2 Effects on spectra

Figure 8: A: Comparison of the average line profile in granules before and
after the deconvolution. Only the Fe I spectral line at 6301.5 Å is shown here
to magnify the differences. B: Comparison reproduced for the inter-granular
lanes. The vertical scale for both plots have been kept the same. Figure 9:
Line profiles of Fe I at $6297$ Å shown before (dotted) and after (solid) the
deconvolution in certain pixels. The grey shaded region highlights the
location of the spike in the blue wing of the spectral line.

The deconvolution was found to affect the shape of spectral lines
significantly. Fig. 8 shows a comparison of the spatially averaged, normalised
line profiles before and after the deconvolution in granules (top) and inter-
granular lanes (bottom) for the Fe I line at $6301.5$ Å. The selection of
granules and inter-granular lanes was made on the basis of an intensity
threshold. Although the differences are diminished due to the averaging,
certain characteristic changes are visible.

The depth of spectral lines in the granules was found to increase, as well as
become more blue-shifted. This is because the stray-light profile in granules
is composed of contributions from the surrounding inter-granular lanes which
are comparatively shallower, and red-shifted due to the presence of down-
flows. Subtracting these contributions, therefore, causes the spectral lines
to become deeper and more blue-shifted in the granules. The converse holds
true for the inter-granular lanes, where removing the stray-light contribution
from granules causes the spectral lines within inter-granular lanes to become
shallower and more red-shifted. These effects collectively indicate that
temperature gradients, as well as line-of-sight velocities in both granules,
and inter-granular lanes are intrinsically larger than what would be observed
without stray-light correction.

In the plage regions, the deconvolution was found to increase the amplitude of
the Stokes $V$ lobes, while decreasing the depth of spectral lines in Stokes
$I$. Since, Stokes $V$ is roughly proportional to the derivative of Stokes $I$
(under the weak-field approximation), deconvolution increases the deduced
field strengths of the plage. However, several degeneracies in the atmospheric
parameters affecting the shape of the spectral lines exist, and an accurate
quantification of the effects of stray-light removal on spectropolarimetric
data requires spectropolarimetric inversions.

The deconvolution was also found to amplify artefacts in certain pixels. Fig.
9 shows a spectral line at $6297.79$ Å in four different pixels (labelled
A-D). The dotted and solid lines show the spectral line before and after
stray-light deconvolution, respectively. In other words, the dotted line
represents slit spectra that were restored with MFBD, and the solid lines
represent slit-spectra that were restored with MFBD, followed by the stray-
light deconvolution. The grey shaded regions indicate the location in the blue
wing of the spectra line where a conspicuous spike appears. These spikes were
found to be visible in all but the telluric lines, and their severity was
found to vary. In some pixels (labelled B and C), the spike was already found
in MFBD-restored data before the stray-light deconvolution. A likely cause for
this may be the differential aberrations between the slit-jaw and the spectral
cameras - when the slit moves or vibrates within the slit-jaw cameras’ FOV
during the scan. This implies that the PSFs returned by MFBD, based on the
slit-jaw images alone, may not be able to precisely describe the degradation
in the spectra recorded by the spectral camera. In pixels labelled A and D, it
is harder to ascertain whether the spike existed before the stray-light
deconvolution. Nonetheless, the stray-light deconvolution appears to have
amplified it. We found the number of pixels displaying such artefacts above
the noise level to be less than $0.1\%$ of the FOV.

##  8 Discussions and conclusions

In this paper, we extended the method initially developed by Scharmer et al.
(2010) \- using Kolmogorov statistics - to self consistently estimate the
amount of residual stray light directly from the observations. Although we
focused on slit-spectra, our method is also applicable to observations
obtained with filtergraphs. We found that the RMS contrasts of granulation in
our observations increased by roughly three percentage points - from $9\%$ to
$12.5\%$. This is much closer to the values of contrasts found in MHD
simulations. Since the RMS contrasts in MHD simulations themselves show an
intrinsic variability of one percent, the discrepancy between our observations
and simulations may be of the order of one to three percent.

In our method, we made several deliberate choices that theoretically lead to a
slight underestimation of the residual stray light in our data for reasons of
expediency, but it is possible to relax some of these constraints. For
example, the estimation of the stray-light PSFs can also be performed as a
function of position along the slit. One would therefore have a set of fully
space variant stray-light PSFs. This requires more computational resources as
multiple independent MFBD restorations need to be performed by segmenting the
FOV into patches, but it would conceptually be similar to the method we have
implemented here.

We note that our method only removes the component of stray light due to high-
order aberrations. Stray-light contamination from instrumental sources, and
optical scattering remain unknown in our observations. Progress, however, is
being made in accounting for these components. Löfdahl & Scharmer (2012) used
a setup with an artificial target - situated at the primary of the focus of
the telescope - to ascertain the amount of stray light originating in the
post-focus optics at SST. They reported that contributions to stray light due
to ghost images and scattering may be up to a few percent, in addition to the
dominant contribution coming from high-order wave-front aberrations
originating in the AO bimorph mirror.

Due to the lack of a comprehensive solution that accounts for all the
components of stray light, our observations are still likely to contain some
residual degradation. However, given that our method accounts for, and
restores a good chunk of the missing contrast, we consider our observations to
be suitable for spectropolarimetric inversions.

###### Acknowledgements.

The authors would like to thank Peter Sütterlin, and Hans-Peter Dörr for their
technical support and useful discussions. This project has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation program (grant agreement No. 695075) and has been
supported by the BK21 plus program through the National Research Foundation
(NRF) funded by the Ministry of Education of Korea. The participation of S.
Saranathan was funded by the International Max Planck Research School for
Solar System Science.

## References

  * Beck et al. (2011) Beck, C., Rezaei, R., & Fabbian, D. 2011, A&A, 535, A129
  * Borrero, J. M. et al. (2016) Borrero, J. M., Asensio Ramos, A., Collados, M., et al. 2016, A&A, 596, A2
  * Danilovic et al. (2008) Danilovic, S., Gandorfer, A., Lagg, A., et al. 2008, A&A, 484, L17
  * David & Elste (1962) David, K. H. & Elste, G. 1962, ZAp, 54, 12
  * Delbouille et al. (1973) Delbouille, L., Roland, G., & Neven, L. 1973, Atlas photometrique du spectre solaire de [lambda] 3000 a [lambda] 10000
  * Fried (1966) Fried, D. L. 1966, J. Opt. Soc. Am., 56, 1372
  * Frutiger et al. (2000) Frutiger, C., Solanki, S. K., Fligge, M., & Bruls, J. H. M. J. 2000, A&A, 358, 1109
  * Glindemann et al. (1993) Glindemann, A., Lane, R. G., & Dainty, J. C. 1993, Journal of Modern Optics, 40, 2381
  * Iglesias et al. (2016) Iglesias, F. A., Feller, A., Nagaraju, K., & Solanki, S. K. 2016, A&A, 590, A89
  * Keller & von der Luehe (1992) Keller, C. U. & von der Luehe, O. 1992, A&A, 261, 321
  * Kiselman et al. (2011) Kiselman, D., Pereira, T., Gustafsson, B., et al. 2011, A&A, 535, A14
  * Lane et al. (1992) Lane, R. G., Glindemann, A., & Dainty, J. C. 1992, Waves in Random Media, 2, 209
  * Löfdahl & Scharmer (1994) Löfdahl, M. G. & Scharmer, G. B. 1994, A&AS, 107, 243
  * Löfdahl & Scharmer (2012) Löfdahl, M. G. & Scharmer, G. B. 2012, A&A, 537, A80
  * Martínez Pillet (1992) Martínez Pillet, V. 1992, Sol. Phys., 140, 207
  * Mathew et al. (2009) Mathew, S. K., Zakharov, V., & Solanki, S. K. 2009, A&A, 501, L19
  * Mattig (1971) Mattig, W. 1971, Sol. Phys., 18, 434
  * McGlamery (1976) McGlamery, B. L. 1976, in Proc. SPIE, Vol. 74, Image processing, ed. J. C. Urbach, 225–233
  * Nagaraju & Feller (2012) Nagaraju, K. & Feller, A. 2012, arXiv e-prints [arXiv:1211.0460]
  * Nicholls et al. (1995) Nicholls, T. W., Boreman, G. D., & Dainty, J. C. 1995, Optics Letters, 20, 2460
  * Noll (1976) Noll, R. J. 1976, Journal of the Optical Society of America (1917-1983), 66, 207
  * Orozco-Suárez et al. (2007) Orozco-Suárez, D., Bellot-Rubio, L. R., & del Toro Iniesta, J. C. 2007, The Astrophysical Journal Letters, 662, L31
  * Paxman et al. (1992) Paxman, R. G., Schulz, T. J., & Fienup, J. R. 1992, J. Opt. Soc. Am. A, 9, 1072
  * Rao et al. (2000) Rao, C., Jiang, W., & Ling, N. 2000, Journal of Modern Optics, 47, 1111
  * Rimmele & Marino (2011) Rimmele, T. R. & Marino, J. 2011, Living Reviews in Solar Physics, 8, 2
  * Sánchez Cuberes et al. (2000) Sánchez Cuberes, M., Bonet, J. A., Vázquez, M., & Wittmann, A. D. 2000, ApJ, 538, 940
  * Sarazin & Roddier (1990) Sarazin, M. & Roddier, F. 1990, A&A, 227, 294
  * Scharmer et al. (2003) Scharmer, G. B., Bjelksjo, K., Korhonen, T. K., Lindberg, B., & Petterson, B. 2003, in Proc. SPIE, Vol. 4853, Innovative Telescopes and Instrumentation for Solar Astrophysics, ed. S. L. Keil & S. V. Avakyan, 341–350
  * Scharmer et al. (2010) Scharmer, G. B., Löfdahl, M. G., van Werkhoven, T. I. M., & de la Cruz Rodríguez, J. 2010, A&A, 521, A68
  * Scherrer et al. (2012) Scherrer, P. H., Schou, J., Bush, R. I., et al. 2012, Sol. Phys., 275, 207
  * Staveland (1970) Staveland, L. 1970, Sol. Phys., 12, 328
  * Stumpff (1961) Stumpff, P. 1961, Veroeffentlichungen der Universitaets-Sternwarte zu Goettingen, 7, 305
  * Uitenbroek et al. (2007) Uitenbroek, H., Tritschler, A., & Rimmele, T. 2007, The Astrophysical Journal, 668, 586
  * van Noort (2012) van Noort, M. 2012, A&A, 548, A5
  * van Noort (2017) van Noort, M. 2017, A&A, 608, A76
  * van Noort et al. (2005) van Noort, M., Rouppe van der Voort, L., & Löfdahl, M. G. 2005, Sol. Phys., 228, 191
  * Vögler et al. (2005) Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335
  * Wanders (1934) Wanders, A. J. M. 1934, ZAp, 8, 108
  * Wedemeyer-Böhm & Rouppe van der Voort (2009) Wedemeyer-Böhm, S. & Rouppe van der Voort, L. 2009, A&A, 503, 225
  * Zwaan (1965) Zwaan, C. 1965, in IAU Symposium, Vol. 22, Stellar and Solar Magnetic Fields, ed. R. Lust, 277