# Point spread function modelling for astronomical telescopes: a review
focused on weak gravitational lensing studies

Tobías I. Liaudat 1,2,∗, Jean-Luc Starck 1, Martin Kilbinger 1 and Pierre-
Antoine Frugier 1

###### Abstract

## 1

The accurate modelling of the Point Spread Function (PSF) is of paramount
importance in astronomical observations, as it allows for the correction of
distortions and blurring caused by the telescope and atmosphere. PSF modelling
is crucial for accurately measuring celestial objects’ properties. The last
decades brought us a steady increase in the power and complexity of
astronomical telescopes and instruments. Upcoming galaxy surveys like Euclid
and LSST will observe an unprecedented amount and quality of data. Modelling
the PSF for these new facilities and surveys requires novel modelling
techniques that can cope with the ever-tightening error requirements. The
purpose of this review is three-fold. First, we introduce the optical
background required for a more physically-motivated PSF modelling and propose
an observational model that can be reused for future developments. Second, we
provide an overview of the different physical contributors of the PSF,
including the optic- and detector-level contributors and the atmosphere. We
expect that the overview will help better understand the modelled effects.
Third, we discuss the different methods for PSF modelling from the parametric
and non-parametric families for ground- and space-based telescopes, with their
advantages and limitations. Validation methods for PSF models are then
addressed, with several metrics related to weak lensing studies discussed in
detail. Finally, we explore current challenges and future directions in PSF
modelling for astronomical telescopes.

##  2 Keywords:

point spread function, inverse problems, weak gravitational lensing, image
processing, super-resolution.

##  3 Introduction

Any astronomical image is observed through an optical system that introduces
deformations and distortions. Even the most powerful imaging system introduces
distortions to the observed object. How to characterise these distortions is a
subject of study known as PSF modelling. Specific science applications, like
weak gravitational lensing (WL) in cosmology (see Kilbinger (2015); Mandelbaum
(2018) for reviews), require very accurate and precise measurements of galaxy
shapes. A crucial step of any weak lensing mission is to estimate the PSF at
any position of the observed images. If the PSF is not considered when
measuring galaxy shapes, the measurement will be biased, resulting in
unacceptably biased WL studies. Furthermore, the PSF can be the predominant
source of systematic errors and biases in WL studies. This fact makes PSF
modelling a vital task. Forthcoming astronomical telescopes, such as the
Euclid space telescope (Laureijs et al., 2011), the Roman space telescope
(Spergel et al., 2015; Akeson et al., 2019), and the Vera C. Rubin Observatory
(LSST Science Collaboration et al., 2009; Ivezić et al., 2019), raise many
challenges for PSF models as the instruments are getting more complex and the
imposed scientific requirements tighter. These factors have triggered and
continue to trigger developments in the PSF modelling literature.

Figure 1: An illustration of a field of view showing the PSF modelling
problem. First, the PSF model should be estimated from the stars. The model
should then be used to estimate the PSF at the target positions, e.g., galaxy
positions.

PSF modelling is an interdisciplinary problem which requires knowledge of
optics, inverse problems, and the target science application, in our case,
weak gravitational lensing studies. The objective is to estimate the PSF at
target positions, e.g., galaxy positions, from degraded star observations and
complementary sources of information. Figure 1 shows an illustration of the
problem. The PSF modelling problem is challenging as the model should account
for the different variations of the PSF in the field of view, i.e., spatial,
spectral and temporal. This review is related to these three scientific
fields, discusses in detail the PSF and aims to help understand the different
PSF modelling choices. We start by introducing optical concepts required to
analyse optical imaging systems that are required to understand the more
physically-based PSF models in section 4. Then, motivated by the optical
introduction, we describe the adopted general observational model in section
5. section 6 introduces the different contributors to the PSF at the optical
and detector level. section 7 gives an overview of state-of-the-art PSF
modelling techniques and leads to section 8, which includes comments on the
desirable properties of a PSF model. We end the review by describing different
techniques for validating PSF models in section 9 and concluding in section
10. In addition, we include Table 1, which summarizes the notation and the
different coordinates used throughout this article.

Variable | Description
---|---
Coordinates
$(x,y)$ | Pupil plane or output aperture plane coordinates
$(u,v)$ | Image or focal plane coordinates
$(\xi,\eta)$ | Object plane coordinates
$(\bar{u},\bar{v})$ | Pixel coordinates, the discrete counterpart of the image plane
$\mathbf{p}_{i}$ | 3D spatial coordinate
$\lambda$ | Wavelength
$t$ | Time
Notation
$\mathcal{I},\mathcal{H},\ldots$ | Calligraphic uppercase variables are continuous functions
$I,H,\ldots$ | Uppercase variables are matrices
$c_{m},b_{1}^{k},\ldots$ | Lowercase variables are scalars
$I_{\mathrm{img}}(\bar{u},\bar{v};t|u_{i},v_{i})\in\mathbb{R}$ | Pixel value at position $(\bar{u},\bar{v})$ for the image $I_{\mathrm{img}}$ with its
| centroid at position $(u_{i},v_{i})$ observed at time $t$.
$I_{\mathrm{img},(\cdot|u_{i},v_{i})}\in\mathbb{R}^{p\times p}$ | Observed image with its centroid at position $(u_{i},v_{i})$
Table 1: Coordinates and notation used throughout this article.

##  4 Introduction to optics

This section introduces core optical concepts for understanding the PSF, how
to model it, and certain implicit assumptions usually adopted. The optic
formalism follows the nomenclature of Goodman (2005). For a profound and
rigorous description of optical theory, we refer the reader to the seminal
book of Born and Wolf (1999) or the more concise work seen in Gross (2005). We
refer to Schmidt (2010) for more information on practical wave propagation. If
the reader is familiar with the Fourier optic literature, continue to section
5.

###  4.1 Scalar diffraction theory

When studying the PSF, we are examining how an optical system with a specific
instrument contributes to and modifies our observations. To understand how the
optical system interacts with the propagation of light, we need to dig into
the nature of light, an electromagnetic (EM) wave. To make a fundamental
analysis, one would need to use Maxwell’s equations, solve them with the
optical system under study, and obtain the electric and magnetic fields.
Solving a set of coupled partial differential equations is an arduous task.
Several approximations can be made, if some conditions are met, to alleviate
the mathematical burden of solving Maxwell’s equations without introducing
much error into the analysis.

subfigure0  (a) (b)

() (a) Illustration of the Huygens-Fresnel principle and the modification of a
wavefront due to an obstacle. Reproduced from Liaudat (2022). (b) Illustration
of the different diffraction regions behind an aperture.

Diffraction theory provides a fundamental framework for analysing light
propagation through an optical system. This is especially the case when
working with EM waves in the optical range when the optical image is close to
the focus region. The Huygens-Fresnel principle (Huygens, 1690; Fresnel, 1819;
Crew et al., 1900) states that every point of a wavefront may be considered a
secondary disturbance giving rise to spherical wavelets. At any later instant,
the wavefront may be regarded as the envelope of all the disturbances.
Fresnel’s contribution to the principle is that the secondary wavelets
mutually interfere. This principle provides a powerful method of analysis of
luminous wave propagation. In Figure 2(), the propagation of an incident plane
wavefront through an obstacle, a single slit, is shown. The secondary wavelets
constitute the plane wavefront before the obstacle. Then, the wavefront shape
is modified due to the obstacle, following the Huygens-Fresnel principle.

The secondary waves mutually interfere constructively or destructively,
according to their phases. The analysis of the light propagation in a
homogeneous medium is simple as the spherical wavelets interfere without
obstacles, and the total wavefront propagates spherically in the medium.
However, suppose the wave encounters an obstacle. In that case, the secondary
waves in the vicinity of the boundaries of the obstacle will interfere in ways
that are not obvious from the incident wavefront.

Gustav Kirchhoff was a pioneer in providing a solid mathematic foundation to
the Huygens-Fresnel principle using Green’s theorem (Kirchhoff, 1883). He
started by deriving the integral theorem by Kirchhoff, which expresses the
solution of the homogeneous wave equation at an arbitrary point in terms of
the values of the solution and its first derivative at all points on an
arbitrary closed surface surrounding the point (Born and Wolf, 1999, §4.3.1).
Then, he studied the diffraction by a planar screen, a simple setup which
allows making assumptions simplifying the integral theorem. The assumptions
are commonly known as the Kirchhoff boundary conditions (Kirchhoff, 1883).
Although accurate in practice, these boundary conditions were proved
inconsistent and were eliminated Sommerfeld (1896), giving origin to the
Rayleigh-Sommerfeld diffraction theory (Rayleigh, 1897; Sommerfeld, 1964; Born
and Wolf, 1999, §8.11) that we will continue to describe below.

Figure 2: Illustration of the coordinate system for the diffraction equations.
Figure adapted from Liaudat (2022).

Let us consider a diffractive aperture in a plane $(x,y)$ illuminated in the
positive $z$ direction. We will study the diffracted wave in a parallel plane
$(u,v)$ at a normal distance $z$ from the first plane. The $z$-axis is
orthogonal to both planes and intersects them at their origins. Figure 2
illustrates the coordinate system described above. The diffracted wave is
written as

$\mathcal{U}(\mathbf{p}_{1})=\frac{z}{\text{j}\lambda}\iint_{\Sigma}\mathcal{U}(x,y;0)\frac{\exp(\text{j}\,k\,r_{01})}{r_{01}^{2}}\,\mathrm{d}x\mathrm{d}y\;,$
(1)

where j denotes the imaginary unit, $\lambda$ is the wavelength,
$k=2\pi/\lambda$, $\mathbf{p}_{0}=(x_{0},y_{0};0)$,
$\mathbf{p}_{1}=(u_{1},v_{1};z)$,
$r_{01}=\|\mathbf{p}_{1}-\mathbf{p}_{0}\|_{2}$, $\Sigma$ is the aperture in
the $(x,y)$ plane, and $\mathcal{U}$ is the electric field. The incident wave
is $\mathcal{U}(\mathbf{p}_{0})$ and the diffracted wave is
$\mathcal{U}(\mathbf{p}_{1})$.

There are two approximations in the derivation of Equation 1. The first
approximation is that we are considering a scalar theory of diffraction, a
scalar electric and magnetic field, and not the fields in their complete
vectorial form. The scalar theory provides a full description of the EM fields
in a dielectric medium that is linear, isotropic, homogeneous, and non-
dispersive. However, even if the medium verifies these properties, if some
boundary conditions are imposed on a wave, e.g., an aperture, some coupling is
introduced between the EM field components and the scalar theory is no longer
exact. Nevertheless, the EM fields are modified only at the edges of the
aperture, and the effects extend over only a few wavelengths into the
aperture. Therefore, if the aperture is large compared to the wavelength, the
error introduced by the scalar theory is negligible. Refractive optical
elements can also induce polarisation of the EM field. The level of accuracy
desired will determine if the bias introduced can be neglected or needs to be
taken into account.

Although the current formulation is powerful in representing the diffraction
phenomena, it is still challenging to work with the integral from Equation 1.
As a consequence, we will explore further approximations to the Rayleigh-
Sommerfeld diffraction theory that will give origin to the Fresnel diffraction
and Fraunhofer diffraction.

###  4.2 The Fresnel approximation

The Fresnel approximation is based on the binomial expansion of the square
root in the expression $\sqrt{1+b}$ for some $b$111The binomial expansion is
given by $\sqrt{1+b}=1+\frac{1}{2}b-\frac{1}{8}b^{2}+\cdots$.. The distance
$r_{01}$ can be expressed as

$r_{01}=z\sqrt{1+\left(\frac{u_{1}-x_{0}}{z}\right)^{2}+\left(\frac{v_{1}-y_{0}}{z}\right)^{2}}\;,$
(2)

which can be approximated, using the first two terms of the binomial
expansion, as

$r_{01}\approx
z\left(1+\frac{1}{2}\left(\frac{u_{1}-x_{0}}{z}\right)^{2}+\frac{1}{2}\left(\frac{v_{1}-y_{0}}{z}\right)^{2}\right)\;.$
(3)

The $r_{01}$ appearing in the exponential of Equation 1 has much more
influence in the result than the $r_{01}^{2}$ in the divisor. Therefore, we
use Equation 3 to approximate $r_{01}$ in the exponential, and for the divisor
we approximate $r_{01}^{2}\approx z^{2}$. Then, we can express the diffracted
field as

$\mathcal{U}(u_{1},v_{1};z)=\frac{e^{\text{j}kz}}{\text{j}\lambda
z}\iint_{\Sigma}\mathcal{U}(x,y;0)\exp\left(\text{j}\,\frac{k}{2z}\left[\left(u_{1}-x\right)^{2}+\left(v_{1}-y\right)^{2}\right]\right)\,\mathrm{d}x\mathrm{d}y\;,$
(4)

and if we expand the terms in the exponential, we get

$\mathcal{U}(u_{1},v_{1};z)=\frac{e^{\text{j}kz}}{\text{j}\lambda
z}e^{\text{j}\,\frac{k}{2z}\left(u_{1}^{2}+v_{1}^{2}\right)}\iint_{\Sigma}\left\\{\mathcal{U}(x,y;0)e^{\text{j}\,\frac{k}{2z}\left(x^{2}+y^{2}\right)}\right\\}e^{-\text{j}\,\frac{2\pi}{\lambda
z}\left(u_{1}x+yv_{1}\right)}\mathrm{d}x\mathrm{d}y\;.$ (5)

The Fourier transform (FT) expression can be recognised with some
multiplicative factors in Equation 5. The diffracted wave is the FT of the
product of the incident wave and a quadratic phase exponential. In this case,
we have approximated the spherical secondary waves of the Huygens-Fresnel
principle by parabolic wavefronts. The region where the approximation is valid
is known as the region of Fresnel diffraction. In this region, the major
contributions to the integral come from points $(x,y)$ for which $x\approx u$
and $y\approx v$, meaning that the higher-order terms in the expansion that we
are not considering are unimportant. The region of Fresnel diffraction can be
seen as the coordinates $(u,v,z)$ that verify

$z^{3}\gg\frac{\pi}{4\lambda}\left(\left(u-x\right)^{2}+\left(v-y\right)^{2}\right)^{2}\;,\quad\forall(x,y)\in\Sigma.$
(6)

###  4.3 The Fraunhofer approximation

We continue to present a further approximation that, if valid, can
significantly simplify the calculations. The Fraunhofer approximation assumes
that the exponential term with a quadratic dependence of $(x,y)$ is
approximately unity over the aperture. The region where the approximation is
valid is the far field or Fraunhofer region. Figure 2() illustrates the
different diffraction approximations as a function of the aperture’s distance.
The required condition to be in this region reads

$z\gg\frac{k\left(x^{2}+y^{2}\right)}{2}\;,\quad\forall(x,y)\in\Sigma.$ (7)

The Fraunhofer diffraction formula is given by

$\mathcal{U}(u_{1},v_{1};z)=\frac{e^{\text{j}kz}}{\text{j}\lambda
z}e^{\text{j}\,\frac{k}{2z}\left(u_{1}^{2}+v_{1}^{2}\right)}\iint_{\Sigma}\left\\{\mathcal{U}(x,y;0)\right\\}e^{-\text{j}\,\frac{2\pi}{\lambda
z}\left(u_{1}x+yv_{1}\right)}\,\mathrm{d}x\mathrm{d}y\;,$ (8)

where we can reformulate the previous equation using the FT as follows

$\mathcal{U}(u_{1},v_{1};z)=\frac{e^{\text{j}kz}}{\text{j}\lambda
z}e^{\text{j}\,\frac{k}{2z}\left(u_{1}^{2}+v_{1}^{2}\right)}\text{FT}\left\\{i_{\Sigma}(x,y;0)\mathcal{U}(x,y;0)\right\\}\left(\frac{u_{1}}{\lambda
z},\frac{v_{1}}{\lambda z}\right)\;,$ (9)

where $i_{\Sigma}$ is an indicator function over the aperture taking values in
$\\{0,1\\}$. It is also possible to consider image vignetting and multiply the
indicator with a weight function so that the resulting function takes values
in $[0,1]$. Cameras are sensitive to the light’s intensity reaching their
detectors. The instantaneous intensity of an EM wave is equal to its squared
absolute value. Therefore, we can write the intensity of the diffracted wave
as

$\mathcal{I}(u_{1},v_{1};z)=\left|\mathcal{U}(u_{1},v_{1};z)\right|^{2}=\frac{1}{\lambda^{2}z^{2}}\left|\text{FT}\left\\{i_{\Sigma}(x,y;0)\mathcal{U}(x,y;0)\right\\}\left(\frac{u_{1}}{\lambda
z},\frac{v_{1}}{\lambda z}\right)\right|^{2}\;,$ (10)

which is significantly simpler than the original Rayleigh-Sommerfeld
expression from Equation 1.

###  4.4 Analysis of optical imaging systems

Figure 3: Illustration of the generalised black box concept of the imaging
system and its corresponding coordinate systems. The image plane coordinates
are $(u,v)$, the output aperture plane coordinates are $(x,y)$, and the object
plane coordinates are $(\xi,\eta)$. This figure has been adapted from Goodman
(2005).

Let us analyse a general optical imaging system composed of one or many lenses
or mirrors of possibly different characteristics. We treat the optical system
as a black box characterised by the transformations applied to an incident
object scalar wave, $\mathcal{U}_{\rm obj}$, into an output image wave,
$\mathcal{U}_{\rm img}$. Figure 3 illustrates the black box interpretation of
the general optical system. We assume that the effect of the optical system
between the input and output waves is well described by geometrical optics,
which is an affine transformation. All the diffraction effects can be
associated with one of the two pupils, input or output (see (Goodman, 2005, S.
6.1.2) for more discussion). We choose the latter one and consider the
diffraction of the output wave between the output pupil and the image plane.

Let us write the output wave as a function of the input wave using the
superposition integral as follows

$\mathcal{U}_{\rm
img}(u,v)=\iint_{-\infty}^{+\infty}h(u,v;\xi,\eta)\mathcal{U}_{\rm
obj}(\xi,\eta)\mathrm{d}\xi\mathrm{d}\eta\;,$ (11)

where $h$ is the field’s value at image coordinates $(u,v)$ due to a unitary
point-source object at position $(\xi,\eta)$.

Figure 4: Illustration of the isoplanatic region. Two rays from the same
isoplanatic region travel through almost the same turbulence and suffer almost
the same distortions. Figure reproduced from Liaudat (2022).

On the one hand, it is cumbersome to work with non-stationary systems as $h$
in Equation 11. On the other hand, the study of $h$ is the main topic of this
review, which is strongly non-stationary in systems with a large field-of-view
(FOV). However, despite being central to this work, we can consider $h$ to be
stationary in its isoplanatic region. This region consists of close-by points
in the FOV, where the light has travelled similar paths giving small
deviations of $h$. We are assuming a certain regularity in $h$ due to the
physical system under study that allows the deviations to be small. In other
words, we consider $h$ to be locally stationary or stationary in patches.
Figure 4 illustrates the idea of an isoplanatic region. This assumption limits
the size of the objects being imaged, as they should have a certain size range
with respect to the support of $h$ so that all the objects being imaged lie
within the aforementioned region. Then, for a set of positions
$\\{(u,v),(\xi,\eta)\\}$, assuming the stationarity of $h$, we have that
$h(u,v;\xi,\eta)=h(u-\xi,v-\eta)$. A shift in the object space translates as a
shift in the image space without a change of $h$’s morphology.

The effect of the geometrical optics inside the black box on the input
(object) wave writes

$\mathcal{U}_{\rm
g}\left(\tilde{\xi},\tilde{\eta}\right)=\frac{1}{|m|}\mathcal{U}_{\rm
obj}\left(\frac{\tilde{\xi}}{|m|},\frac{\tilde{\eta}}{|m|}\right)\;,\quad\text{and}\quad\tilde{\xi}=m\xi,\quad\tilde{\eta}=m\eta\;,$
(12)

where $m$ is the magnification factor of the optical system. Following the
aforementioned approximations, we can rewrite Equation 11 as a convolution as
follows

$\mathcal{U}_{\rm
img}(u,v)=\iint_{-\infty}^{+\infty}h\left(u-\tilde{\xi},v-\tilde{\eta}\right)\mathcal{U}_{\rm
g}\left(\tilde{\xi},\tilde{\eta}\right)\;\mathrm{d}\tilde{\xi}\mathrm{d}\tilde{\eta}\;,$
(13)

where the impulse response $h$, under the Fraunhofer approximation, writes

$h(u,v)=\frac{a}{\lambda\,z_{i}}\iint_{-\infty}^{+\infty}\mathcal{P}(x,y)e^{-\text{j}\frac{2\pi}{\lambda\,z_{i}}\left(ux+vy\right)}\mathrm{d}x\mathrm{d}y\;,$
(14)

where $\mathcal{P}$ is the pupil function, a binary mask indicating the
obscurations in the exit pupil’s aperture, $a$ is a constant amplitude, and
$z_{i}$ is the distance between the exit pupil plane and the image plane. In
this general case, where the system has no aberrations, we see that the output
image is formed by a geometrical-optics transformation followed by a
convolution with an impulse response that is the Fraunhofer diffraction of the
exit aperture.

####  4.4.1 Introducing optical aberrations

Figure 5: Illustration of the wavefront error in a one-dimensional setting
where the optical system is represented as a single lens. Figure reproduced
from Liaudat et al. (2023).

We considered in the previous development an ideal optical system without any
aberrations which is known to be diffraction-limited. An aberrated optical
system produces the imperfect convergence of rays, which can be expressed
equivalently in wavefront space by deviations from the ideal reference sphere.
The aberrations produce leads and lags in the wavefront with respect to the
ideal sphere, see Figure 5. A complementary interpretation, from Goodman
(2005), is that we start with the previous diffraction-limited system
producing converging spherical wavefronts. Then, we add a phase-shifting plate
representing the system’s aberrations. The plate is located in the aperture
after the exit pupil and deforms the output wave’s phase. To characterise the
aberrations, we will use the generalised pupil function that generalises the
pupil function $\mathcal{P}$ from Equation 14 and writes

$\mathcal{P}_{\rm
gen}(x,y)=\mathcal{P}(x,y)e^{\text{j}\frac{2\pi}{\lambda}\mathcal{W}(x,y)}\;,$
(15)

where $\lambda$ is the central wavelength of the incident wave, $\mathcal{P}$
is the previous pupil function, an obscuration mask, and $\mathcal{W}$
represents the optical path differences (OPD) between a perfectly spherical
and the aberrated wavefront. We also refer to the OPD as the wavefront error
(WFE). Figure 5 illustrates the concept of WFE. It is common to represent the
WFE using a Zernike polynomial decomposition (Noll, 1976) as they are
orthogonal in the unit disk and we generally use circular apertures in
telescopes and optical systems in general. Figure 6() shows the first Zernike
polynomials.

subfigure0  (a) (b)

() (a) Euclid’s pupil function, or obscurations, in the $(x,y)$ plane for a
given position in the $(\xi,\eta)$. (b) Example of the first Zernike
polynomial maps.

The aberrations, $\mathcal{W}$, and the pupil function, $\mathcal{P}$, depend
on the object’s position in the FOV, even though we have not introduced the
dependence in the formulations. To illustrate this fact, we can consider the
obscurations in the aperture as seen in Figure 6(). One can notice a circular
aperture with several obscurations in it, a secondary mirror and three arms
supporting the mirror. The projection of these obscurations into the $(x,y)$
plane depends on the object’s position in the FOV, making the pupil function
$\mathcal{P}$ change with the $(\xi,\eta)$ coordinates.

We have considered aberrations that only depend on the object’s position in
the FOV. However, depending on the optical system under study, there might be
wavelength-dependent aberrations. For example, some refractive components, or
some components implementing complex thin film coatings, may introduce
spurious spectral dependences to the optical system’s response. If this is the
case, we can add a wavelength dependence to the WFE function $\mathcal{W}$ to
account for these effects.

####  4.4.2 Polychromatic illumination: the coherent and the incoherent case

We studied until now a system with ideal monochromatic light. It is time to
shift to polychromatic light as telescopes have filters with finite bandwidths
and hence allow multiple frequencies of light. For a more rigorous analysis of
polychromatic illumination, we refer the reader to the theory of partial
coherence Beran and Parrent (1964, §10); Goodman (1985, §10); Born and Wolf
(1999, §10). Even if we study the system’s behaviour to light with a
particular wavelength, this is practically never the case, as real
illumination is never perfectly chromatic, even for lasers. Therefore, we
consider a narrowband polychromatic illumination centred at a given wavelength
$\lambda$. The narrowband assumption states that the bandwidth occupied is
small with respect to the central wavelength. For polychromatic light, we
follow Goodman (2005) and consider a time-varying phasor of the field,
$\mathcal{U}_{\rm img}(u,v;t)$, where its intensity is given by the time
integration of its instantaneous intensity

$\mathcal{I}_{\rm img}(u,v)=\left\langle\left|\mathcal{U}_{\rm
img}(u,v;t)\right|^{2}\right\rangle_{t}=\frac{1}{T}\int_{-T/2}^{T/2}\left|\mathcal{U}_{\rm
img}(u,v;t)\right|^{2}\mathrm{d}t\;,$ (16)

where $T$ is the detector integration time that is considered much greater
than the optical wave period. We can generalise the field expression from
Equation 13 by considering that light is polychromatic and that the impulse
response $h$ is wavelength independent due to the narrowband assumption. Then,
the field writes

$\mathcal{U}_{\rm
img}(u,v;t)=\iint_{-\infty}^{+\infty}h\left(u-\tilde{\xi},v-\tilde{\eta}\right)\mathcal{U}_{\rm
g}\left(\tilde{\xi},\tilde{\eta};t-\tau\right)\;\mathrm{d}\tilde{\xi}\mathrm{d}\tilde{\eta}\;,$
(17)

where $\tau$ represents the delay of the wave propagation from
$(\tilde{\xi},\tilde{\eta})$ to $(u,v)$. Continuing with the polychromatic
analysis, we rewrite the intensity from Equation 16 as

$\mathcal{I}_{\rm
img}(u,v)=\iint\limits_{-\infty}^{+\infty}\mathrm{d}\tilde{\xi}_{1}\mathrm{d}\tilde{\eta}_{1}\iint\limits_{-\infty}^{+\infty}\mathrm{d}\tilde{\xi}_{2}\mathrm{d}\tilde{\eta}_{2}h\left(u-\tilde{\xi}_{1},v-\tilde{\eta}_{1}\right)h^{*}\left(u-\tilde{\xi}_{2},v-\tilde{\eta}_{2}\right)\mathcal{J}_{\rm
g}\left(\tilde{\xi}_{1},\tilde{\eta}_{1};\tilde{\xi}_{2},\tilde{\eta}_{2}\right),$
(18)

where $h^{*}$ is the conjugate of $h$, $\mathcal{J}_{\rm g}$ is known as the
mutual intensity which describes the spatial coherence of $\mathcal{U}_{\rm
g}$ at two points and writes

$\mathcal{J}_{\rm
g}\left(\tilde{\xi}_{1},\tilde{\eta}_{1};\tilde{\xi}_{2},\tilde{\eta}_{2}\right)=\left\langle\mathcal{U}_{\rm
g}\left(\tilde{\xi}_{1},\tilde{\eta}_{1};t\right)\mathcal{U}_{\rm
g}^{*}\left(\tilde{\xi}_{2},\tilde{\eta}_{2};t\right)\right\rangle\;.$ (19)

We can distinguish two types of illuminations, coherent and incoherent.
Coherent illumination refers to waves whose phases vary in a perfectly
correlated way. This illumination is approximately the case of a laser. In
incoherent illumination, the wave’s phases vary in an uncorrelated fashion.
Most natural light sources can be considered incoherent sources. The mutual
intensity is helpful to represent both types of illumination. In the case of
coherent light, we obtain,

$\mathcal{J}_{\rm
g}^{\text{co}}\left(\tilde{\xi}_{1},\tilde{\eta}_{1};\tilde{\xi}_{2},\tilde{\eta}_{2}\right)=\mathcal{U}_{\rm
g}\left(\tilde{\xi}_{1},\tilde{\eta}_{1}\right)\mathcal{U}_{\rm
g}^{*}\left(\tilde{\xi}_{2},\tilde{\eta}_{2}\right)\;,$ (20)

where $\mathcal{U}_{\rm g}\left(\tilde{\xi}_{1},\tilde{\eta}_{1}\right)$ and
$\mathcal{U}_{\rm g}\left(\tilde{\xi}_{2},\tilde{\eta}_{2}\right)$ are time-
independent phasor amplitudes relative to their time-varying counterpart. As
both time-varying phasors are synchronized, we have taken a reference phasor
and normalized them against their amplitude with respect to a reference point
that can be the origin $(0,0)$. For example,

$\mathcal{U}_{\rm
g}\left(\tilde{\xi}_{1},\tilde{\eta}_{1};t\right)=\mathcal{U}_{\rm
g}\left(\tilde{\xi}_{1},\tilde{\eta}_{1}\right)\frac{\mathcal{U}_{\rm
g}\left(0,0;t\right)}{\left\langle\left|\mathcal{U}_{\rm
g}\left(0,0;t\right)\right|^{2}\right\rangle^{\frac{1}{2}}}\;.$ (21)

Substituting Equation 20 into Equation 18 we obtain

$\displaystyle\mathcal{I}_{\rm img}^{\text{co}}(u,v)=$
$\displaystyle\left|\mathcal{U}_{\rm
img}^{\text{co}}(u,v)\right|^{2}=\left|\iint_{-\infty}^{+\infty}h\left(u-\tilde{\xi},v-\tilde{\eta}\right)\mathcal{U}_{\rm
g}\left(\tilde{\xi},\tilde{\eta}\right)\mathrm{d}\tilde{\xi}\mathrm{d}\tilde{\eta}\right|^{2}\;,$
(22) $\displaystyle\mathcal{I}_{\rm img}^{\text{co}}(u,v)=$
$\displaystyle\left|\left(\mathcal{U}_{\rm g}\star
h\right)(u,v)\right|^{2}\;,$

where we observe the coherent illumination gives a linear system in the
complex amplitude of the field $\mathcal{U}_{\rm g}$. The previous result is
related to the interference of coherent waves. If we now consider incoherent
illumination, the mutual intensity writes

$\mathcal{J}_{\rm
g}^{\text{in}}\left(\tilde{\xi}_{1},\tilde{\eta}_{1};\tilde{\xi}_{2},\tilde{\eta}_{2}\right)=\kappa\,\mathcal{I}_{g}\left(\tilde{\xi}_{1},\tilde{\eta}_{1}\right)\delta\left(\tilde{\xi}_{1}-\tilde{\xi}_{2},\tilde{\eta}_{1}-\tilde{\eta}_{2}\right)\;,$
(23)

where $\kappa$ is a real constant, $\delta$ is Dirac delta distribution, and
$\mathcal{I}_{\rm g}$ is the intensity of the $U_{\rm g}$ field. The constant
$\kappa$ is a result of a simplification from statistical optics giving origin
to Equation 23. The constant depends on the degree of the extension of
coherence when the evanescent-wave phenomenon (Beran and Parrent, 1964) is
taken fully into account. If the coherence extends over a wavelength, $\kappa$
is equal to $\bar{\lambda}^{2}/\pi$, where $\bar{\lambda}$ is the mean
wavelength. See Goodman (1985, §5.5.2) for a deeper discussion on incoherent
illumination and the $\kappa$ constant.

Replacing Equation 23 in Equation 18 the output (image) intensity writes

$\displaystyle\mathcal{I}_{\rm img}^{\text{in}}(u,v)$
$\displaystyle=\left|\mathcal{U}_{\rm
img}^{\text{in}}(u,v)\right|^{2}=\kappa\iint_{-\infty}^{+\infty}\left|h\left(u-\tilde{\xi},v-\tilde{\eta}\right)\right|^{2}\mathcal{I}_{\rm
g}\left(\tilde{\xi},\tilde{\eta}\right)\mathrm{d}\tilde{\xi}\mathrm{d}\tilde{\eta}\;,$
(24) $\displaystyle\mathcal{I}_{\rm img}^{\text{in}}(u,v)$
$\displaystyle=\kappa\left(\mathcal{I}_{\rm
g}\star\left|h\right|^{2}\right)(u,v)\;,$

where $\left|h\right|^{2}$ is the intensity impulse response, also known as
the PSF. In this case, an optical system illuminated with incoherent light is
linear in intensity. Equation 24 shows a commonly used fact; the output
intensity is the convolution of the intensity PSF with ideal image intensity
$\mathcal{I}_{\rm g}$.

####  4.4.3 Usual assumptions adopted in PSF modelling

PSF modelling articles generally implicitly assume specific hypotheses. We
provide some of them on the following list:

  * •

The electric field is assumed to be scalar.

  * •

The incoming light from natural sources is assumed to be ideally incoherent.
Then, the optical system is linear in intensity, as seen in Equation 24, and
the PSF is a convolutional kernel.

  * •

The PSF is considered to be spatially invariant in its isoplanatic region. In
other words, PSF is assumed not to change on objects’ typical length scales.

Although the previous assumptions are standard, certain precision levels
require dropping simplifications. The Euclid mission requirements on the PSF
model accuracy is of $2\times 10^{-4}$ for the root mean square (RMS) error on
each ellipticity component ($\delta e_{i}^{\text{PSF}}$), and $1\times
10^{-3}$ for the relative RMS error on the size ($\delta
R_{\text{PSF}}^{2}/R_{\text{PSF}}^{2}$). The PSF model might need to include
light polarisation to fulfil these extremely tight PSF requirements. Other
assumptions might also be dropped for the precise imaging of widespread
objects. This case might require discarding the spatially invariant assumption
of the PSF in its isoplanatic region.

To conclude, the usual formulation of the PSF, i.e., the intensity of the
impulse response, convolving an image seen in many articles, comes from the
previous assumptions using the results from Equation 24, Equation 14 and
Equation 15. We rewrite this formula as follows

$\mathcal{I}_{\rm
img}(u,v)=\left(\mathcal{H}\star\mathcal{I}_{g}\right)(u,v)\;,$ (25)

where we remind the reader that $(u,v)$ is the image plane, we have dropped
the $\kappa$ term from Equation 24, and $\mathcal{H}$ is the PSF that writes

$\mathcal{H}(u,v)=\frac{a^{2}}{\lambda^{2}\,z_{i}^{2}}\left|\iint_{-\infty}^{+\infty}\mathcal{P}(x,y)\,e^{\text{j}\frac{2\pi}{\lambda}\mathcal{W}(x,y)}\,e^{-\text{j}\frac{2\pi}{\lambda\,z_{i}}\left(ux+vy\right)}\,\mathrm{d}x\mathrm{d}y\right|^{2}\;,$
(26)

where we are studying the PSF for a specific wavelength and FOV position.

##  5 General observational model

We consider the PSF as the intensity impulse response of the imaging system
under study to a point source. The concept of PSF (Born and Wolf, 1999) is
used throughout many imaging applications, including astronomical imaging
(Liaudat et al., 2023; Schmitz, 2019), medical imaging (Dougherty and Kawaf,
2001; Joyce et al., 2018), or microscopy (Soulez et al., 2012; Debarnot et
al., 2021a, b). The central idea behind a PSF is that it represents
transformations done to the imaged object by the imaging system. The PSF is,
in a certain way, a characterisation of the imaging system. Considering
incoherent illumination and the hypotheses from the previous section hold, we
can affirm that the optical system behaves linearly as in Equation 25.
Consequently, the PSF is considered the impulse response of the optical system
and affects the ground truth image through a convolution operation. Focusing
on astronomical imaging, the definition of the imaging system can vary between
the different use cases and telescopes. For example, in a ground-based
telescope, we will consider that the atmosphere belongs to the imaging system
we are modelling. However, naturally, the atmosphere will not be considered in
a space-based telescope. This article will focus on optical systems, which
work with electromagnetic radiation with a wavelength close to the visible
spectrum. For example, Euclid VIS instrument’s theoretical wavelength range is
from $550$nm to $900$nm.

The PSF describes the effects of the imaging system in the imaging process of
the object of interest. The PSF is a convolutional kernel, as we have seen in
subsection 4.4. However, this convolutional kernel varies spatially,
spectrally, and temporally. We give a non-exhaustive list that motivates each
of these variations:

  * •

Spatial variations: The optical system presents a certain optical axis, which
is an imaginary line where the system has some degree of rotational symmetry.
In simpler words, it can be considered as the direction of the light ray that
produces a PSF in the centre of the focal plane for an unaberrated optical
system. The angle of incidence is defined as the angle between an incoming
light ray and the optical axis. The main objective of the optical systems we
study is to make the incoming light rays converge in the focal plane, where
there will be some measurement instruments, e.g., a camera. Depending on the
angle of incidence, the image will form in different positions in the focal
plane. The path of the incoming light will be different for each angle of
incidence, and therefore the system’s response will be different too. In other
words, the PSF will change depending on the angle of incidence or spatial
position in the focal plane where the image is forming. Optical systems with
wide focal planes, generally associated with wide FOVs, present significant
PSF spatial variations.

  * •

Spectral variations: Principally due to the diffraction phenomena and its
well-known wavelength dependence covered in section 4. Refractive222Refraction
refers to the change of direction in the propagation of a wave passing from
one medium to another. Most of the wave energy is transmitted to the new
medium. Reflection refers to the abrupt change of direction of the wave
propagation due to a boundary between mediums. In this last case, most of the
oncoming wave energy remains in the same medium. components of the optical
system under study are also a source of spectral variations (Baron et al.,
2022). Other sources of spectral variations are detector electronic components
(Meyers and Burchat, 2015a) and atmospheric chromatic effects (Meyers and
Burchat, 2015b).

  * •

Temporal variations: The state of the telescope changes with respect to time;
therefore, the imaged object’s transformation also changes. In space-based
telescopes, high-temperature gradients cause mechanical dilations and
contractions that affect the optical system. In ground-based telescopes, the
atmosphere composition changes with time. Consequently, it temporally affects
the response of the optical system, i.e., the PSF.

The PSF convolutional kernel varies with space, time and wavelength. Once we
have set up a specific wavelength and time to analyse our system, we will have
a different convolutional kernel for each position in the field of view. Let
us refer to the PSF field $\mathcal{H}$ as all the PSF representing an optical
system. Then $\mathcal{H}(u,v;\lambda;t|u_{i},v_{i})$ is a specific PSF where
$(u_{i},v_{i})$ represents its centroid, i.e., the first order moments. The
same notation is maintained from Figure 3, where the $(u,v)$ variables
represent the image plane. We can define the PSF field as a varying
convolutional kernel
$\mathcal{H}:\mathbb{R}^{2}\times\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathbb{R}^{2}\to\mathbb{R}$.
This definition would accurately describe how the PSF affects the images
considering the assumptions from section 4 are valid. It is, however,
cumbersome to work with a convolutional kernel that varies at each point in
space. We adopt the approximation that considers the PSF locally invariant in
its isoplanatic region (Born and Wolf, 1999, § 9.5.1), see Figure 4 for an
illustration. This approximation means that in a vicinity of an observed
object, we will consider that the PSF only varies with time and wavelength,
thus facilitating the computation of the convolution. The close vicinity, or
the isoplanatic region, will be defined as the postage stamp to image the
object of interest. The typical galaxies observed for weak lensing have a
comparable size with respect to the PSF size (see Mandelbaum et al. (2018,
Figure 7) for a distribution of relative galaxy to PSF size in the HSC
survey). Consequently, the approximation error is kept low as it is only done
for small patches of the focal plane.

Let us define our object of interest with the subscript ground truth (GT),
$\mathcal{I}_{\mathrm{GT}}(u,v;\lambda;t|u_{i},v_{i})$, that is the $I_{g}$
object from section 4, as a continuous light distribution
$\mathcal{I}_{\mathrm{GT}}:\mathbb{R}^{2}\times\mathbb{R}_{+}\times\mathbb{R}_{+}\times\mathbb{R}^{2}\to\mathbb{R}$.
In this review, we are not considering transient objects, i.e., the time
dependence scale of the object is comparable with the exposure time used to
image it. Therefore, we can ignore the temporal dependency of the GT object,
$\mathcal{I}_{\mathrm{GT}}(u,v;\lambda;t)\neq f(t)$. Let us write our general
observational model that relates our GT object of interest, our PSF and our
observed image as follows

$I_{\mathrm{img}}(\bar{u},\bar{v};t|u_{i},v_{i})=\mathcal{F}_{p}\left\\{\int_{0}^{+\infty}\mathcal{T}(\lambda)\;(\mathcal{I}_{\mathrm{GT}}\star\mathcal{H})(u,v;\lambda;t|u_{i},v_{i})\;\mathrm{d}\lambda\right\\}\circ
N(\bar{u},\bar{v};t|u_{i},v_{i})\;,$ (27)

where $\mathcal{F}_{p}$ is a degradation operator discretising the image to
$\mathbb{R}^{p\times p}$ which includes the image sampling from the
instrument. The variables $(\bar{u},\bar{v})$ denote the discrete (pixelised)
version of the $(u,v)$ variables. Then,
$I_{\mathrm{img}}(\bar{u},\bar{v};t|u_{i},v_{i})\in\mathbb{R}$ corresponds to
the instrument’s measurement at a single pixel $(\bar{u},\bar{v})$, and
$I_{\mathrm{img},(\cdot;t|u_{i},v_{i})}\in\mathbb{R}^{p\times p}$ to the
entire image. The variables $(u_{i},v_{i})$ correspond to the centre location
of the target object $i$. The instrument’s transmission is represented by
$\mathcal{T}:\mathbb{R}_{+}\to\mathbb{R}_{+}$, a function with finite support,
and $N_{(\cdot;t|u_{i},v_{i})}\in\mathbb{R}^{p\times p}$ corresponds to the
noise affecting our observation and possibly a modelling error, where $\circ$
is some composition operator. We have carried out the spectral integration
Hopkins (1957); Eriksen and Hoekstra (2018) on the instrument’s passband
defined in $\tau$. Although Equation 27 provides a general observational
model, it can be unpractical. The continuous functions $\mathcal{H}$,
$\mathcal{T}$, and $\mathcal{I}_{\text{G}T}$ are practically inaccessible. We
make several assumptions to simplify the problem:

  * (a)

the continuous functions $\mathcal{H}$ and $\mathcal{I}_{\text{G}T}$ are well
approximated by piece-wise constant functions over a regular grid in
$\mathbb{R}^{2}$. We assume $\mathcal{H}\approx H$ and
$\mathcal{I}_{\text{G}T}\approx I_{\text{G}T}$, where
$H,I_{\text{G}T}\in\mathbb{R}^{P\times P}$ with $P\geq p$. The resolution of
these two variables has to be greater or equal to the observation resolution,

  * (b)

the noise is additive, i.e. $\circ\equiv+$, although the formulation could be
adapted to consider other types of noise, e.g., Poisson,

  * (c)

the degradation operator is approximated by its discrete counterpart,
$\mathcal{F}_{p}\approx F_{p}$, where $F_{p}:\mathbb{R}^{P\times
P}\to\mathbb{R}^{p\times p}$, that has been discretized in a regular grid. We
assume that the degradation operator is linear, and that includes pixellation,
possibly downsampling, intra-pixel shifts and linear detector effects,

  * (d)

we keep the approximation that the PSF is locally constant within the postage
stamp of $P\times P$ values of the target image,

  * (e)

the integral can be well approximated by a discretised version using
$n_{\lambda}$ bins.

Taking into account the aforementioned assumptions, we can define our
practical observational model as follows

$I_{\mathrm{img}}(\bar{u},\bar{v};t|u_{i},v_{i})={F}_{p}\left\\{\sum_{k=1}^{n_{\lambda}}T(\lambda_{k})\;({I}_{\mathrm{GT}}\star{H})(\bar{u},\bar{v};\lambda_{k};t|u_{i},v_{i})\;\Delta\lambda_{k}\right\\}+N(\bar{u},\bar{v};t|u_{i},v_{i})\;,$
(28)

where $I_{\mathrm{img},(\cdot;t|u_{i},v_{i})}\in\mathbb{R}^{p\times p}$, $T$
is a discretized version of $\mathcal{T}$, and $b^{k}=[b^{k}_{0},b^{k}_{1}]$
is the $k$-th wavelength bin centred in $\lambda_{k}$, with a width of
$\Delta\lambda_{k}=b^{k}_{1}-b^{k}_{0}$.

### Particular case: a star observation

Figure 6: The $3$rd generation set of filters of the MegaCam instrument at
Canada-France-Hawaii Telescope that is currently being used for the Canada-
France Imaging Survey. The transmission filter response includes the full
telescope and $1.25$ airmasses of atmospheric attenuation. The full telescope
includes mirrors, optics, and CCDs.

The case of star observations is of particular interest, as some stars in the
FOV can be considered as a spatial impulse, i.e.,
$\mathcal{I}_{\mathrm{star}}(u,v;\lambda|u_{i},v_{i})=\delta(u_{i},v_{i};\lambda)=f_{(u_{i},v_{i})}(\lambda)$.
Therefore, if we plug the impulse in Equation 28, we obtain a degraded
observation of the PSF field. These observations will be crucial to constrain
the PSF models. Unluckily, we do not always have access to the star’s spectral
variation, $f_{(u_{i},v_{i})}(\lambda)$. However, we dispose of complementary
photometric observations that can be useful to characterise the spectral
variations. These observations provide us with the star’s spectral energy
distribution (SED), which can be defined as the calibrated flux density as a
function of wavelength, usually at low spectral resolution. The photometric
observations are done in several spectral bands. Figure 6 shows the bands from
the MegaCam instrument at the Canada-France-Hawaii Telescope (CFHT)333The
filter curves can be downloaded from the Spanish Virtual Observatory (SVO)
webpage, http://svo2.cab.inta-csic.es/svo/theory/fps/index.php .. We refer the
reader to Hogg (2022) for more information about SEDs and stellar photometry.
The SED is a normalised low-resolution sampling of the star’s spectral
variations. We can write the SED definition we will use as

$\text{SED}_{b^{k}}(\lambda_{k})=\frac{1}{z_{n_{\lambda}(b)}}\int_{b^{k}_{0}}^{b^{k}_{1}}f_{(u_{i},v_{i})}(\lambda)\;\mathrm{d}\lambda\;,$
(29)

where we continued to use the $b^{k}$ bin definition from Equation 28, and
$z_{n_{\lambda}(b)}$ is a constant used so that the SED is normalised to
unity. We have that
$\sum_{k=1}^{n_{\lambda}}\text{SED}_{b^{k}}(\lambda_{k})=1$. We continue by
considering that the GT image in Equation 28 is a star, and we use the
spectral bins from the SED definition to discretise the spectral integration.
Finally, we write the practical star observation model as

$I_{\text{star}}(\bar{u},\bar{v};t|u_{i},v_{i})={F}_{p}\left\\{\sum_{k=1}^{n_{\lambda}}T(\lambda_{k})\text{SED}_{b^{k}}(\lambda_{k})\;{H}(\bar{u},\bar{v};\lambda_{k};t|u_{i},v_{i})\;\Delta\lambda_{k}\right\\}+N(\bar{u},\bar{v};t|u_{i},v_{i})\;,$
(30)

where we consider the star observation
$I_{\text{star},(\cdot;t|u_{i},v_{i})}\in\mathbb{R}^{p\times p}$ as a degraded
version of the PSF field
$\tilde{H}_{(\cdot;t|u_{i},v_{i})}\in\mathbb{R}^{p\times p}$.

##  6 PSF field contributors and related degradations

So far, we have described how the PSF interacts with the images we observe and
how we can model an observation. However, we have not given much information
about the different PSF field contributors and the different degradations
represented by $F_{p}$ in Equation 28 that can occur when modelling
observations. We provide a non-exhaustive list of contributors to the PSF
field, sources of known degradations, and the atmosphere’s effect on our PSF
modelling problem.

###  6.1 Image coaddition

A fundamental contributor to the PSF is the choice of image coaddition scheme.
A coadded image is a composite image created by combining multiple individual
exposures of the same region of the sky in some way. This process can help
increase the signal-to-noise ratio of the observation. Motivated by the
analysis of the LSST data, Mandelbaum et al. (2022) explores different
coaddition schemes and studies how it affects the PSF of the resulting coadded
image. In particular, Mandelbaum et al. (2022) define under which schemes the
coadded image accepts a well-defined PSF, i.e., the observation can be
described by the convolution of an extended object and a uniquely defined
coadded PSF. Bosch et al. (2017) describes the strategy for image and PSF
coaddition in the HSC survey.

###  6.2 Dithering and super-resolution

Dithering consists of taking a series of camera exposures shifted by a
fractional or a few pixel amount. There are several advantages of using a
dithering strategy, which include the removal of cosmic rays and
malfunctioning pixels, improving photometric accuracy, filling the gap between
the detectors, and improving the sampling of the observed scene. Dithering
allows estimating a sampling density of the images that is denser than the
original pixel grid, in other words, to super-resolve the image. Regarding
PSFs, it allows recovering Nyquist sampled PSF from undersampled observations.
Bernstein (2002) studied the effect of dithering and the choice of pixel sizes
in imaging strategies. Naturally, as we will later see, the dithering strategy
is helpful for space-based telescopes thanks to their stability. In ground-
based telescopes, the atmosphere constantly changes the PSF, making a
dithering strategy less effective. However, a dithering strategy can be
helpful if the telescope is equipped with adaptive optics technology, which
will be described in subsection 6.6. An example is the Spectro-Polarimetic
High contrast imager for Exoplanets REsearch (SPHERE) instrument (Beuzit,
J.-L. et al., 2019) built for European Southern Observatory’s (ESO) Very Large
Telescope (VLT) in Chile.

Lauer (1999b) discusses the limiting accuracy effect of undersampled PSFs in
stellar photometry and proposes ways to correct it with dithered data (Lauer,
1999a). Fruchter and Hook (2002) presents the widely-used Drizzle algorithm
that consists of shifting and adding the dithered images onto a finer grid.
Rowe et al. (2011) proposed a linear coaddition method coined Imcom to obtain
a super-resolved image from several undersampled images. Hirata et al. (2023)
later studied the use of Imcom on simulations (Troxel et al., 2023) from the
Roman space telescope, while a companion paper, Yamamoto et al. (2023),
explored its implications for weak lensing analyses. Ngolè et al. (2015)
proposed a super-resolution method coined SPRITE targeting the Euclid mission
based on a sparse regularisation technique. More recent PSF models handle the
undersampling of the observations directly in their algorithm for estimating a
well-sampled PSF field, as we will later see.

###  6.3 Optic-level contributors

These contributors affect the PSF by modifying the wave propagation in the
optical system. In other words, they affect the wavefront’s amplitude and
phase.

subfigure0  (a) (b) (c)

() (a) shows a reduced range of variations to show the surface errors measured
for the Hubble Space Telescope, where the scale has been reduced from $\pm
22$nm to better show details. (b) and (c) show the wavefront error
measurements from the JWST cycle $1$ science operation on the $30$th of July
of $2022$ with the total and the reduced range of variation, respectively.

  * •

Diffraction phenomena and the aperture size: As we have seen in section 4, the
diffraction phenomena happening in the optical system plays an essential role
in the formation of the PSF. The size of the optical system aperture and the
wavelength of the light being studied are of particular interest. Equation 26
shows us that under some approximations, the PSF is the Fourier transform of
the aperture. Therefore, the size of the aperture and the PSF are closely
related. For example, if we consider an ideal circular aperture, its
diffraction pattern is the well-known Airy disk. The relation between the
width of the PSF and the diameter of the aperture is given by

$\theta_{\text{FWHM}}=1.025\frac{\lambda}{d}\;,$ (31)

where $\theta_{\text{FWHM}}$ is the full width at half maximum (FWHM)
expressed in radians, $\lambda$ is the wavelength of the light being studied,
and $d$ is the diameter of the aperture. The width of the PSF is a fundamental
property of an optical system as it defines the resolution of the system. In
other words, the PSF size defines the optical system’s ability to distinguish
small details in the image.

  * •

Optical aberrations: These aberrations are due to imperfections in the optical
elements, for example, a not ideally spherical mirror or a not perfectly
aligning of the optical components. The optical aberrations play a significant
role in the morphology of the PSF and can be modelled using the WFE introduced
in the generalised pupil function from Equation 15. Some aberrations have a
distinctive name, e.g., coma, astigmatism, and defocus, and they represent a
specific Zernike polynomial (Noll, 1976).

  * •

Surface errors or polishing effects: One would ideally like perfectly smooth
surfaces in mirrors and lenses. However, imperfections arise in the optical
surfaces due to imperfect surface polishing. Krist et al. (2011) shows the
measurement of surface errors (SFE) in the Hubble Space Telescope (HST). Gross
et al. (2006, Section 35.2) gives a more in-depth analysis of surface errors
focusing on the tolerancing of SFE. Figure 7() shows the surface errors
measured for the Hubble space telescope (HST). Krist and Burrows (1995)
studied HST’s SFE before and after its iconic repair in $1993$ with parametric
and non-parametric, e.g., Gerchberg and Saxton (1972), phase retrieval
algorithms.

  * •

Obscurations: Complex optical systems have telescope designs where some
elements can obscure some part of the pupil. Obscurations are an essential
contributor to PSF morphology and result from projecting a $3$D structure onto
the $2$D focal plane. The resulting projection depends on the considered
position of the focal plane. Accurate modelling of telescopes with wide-field
imagers, e.g., Euclid, requires the computation of the obscuration’s position
dependence arising from the $3$D projection. The Euclid’s obscurations are
presented in Figure 6(). Fienup et al. (1993); Fienup (1994) studied HST’s
obscuration from phase retrieval algorithms and noticed a misalignment that
caused a pupil shift.

  * •

Stray and scattered light: Optical elements and instruments give rise to light
reaching the detectors. Krist (1995) studied the problem for the HST. Storkey
et al. (2004) developed methods to clean observations with scattered light
from the SuperCOSMOS Sky Survey Hambly et al. (2001). Sandin, Christer (2014)
studied the effect of scattered light on the outer parts of the PSF.

  * •

Material outgassing and ice contamination: Material outgassing leads to
molecular contamination that alters different properties of the imaging
system. Water is the most common contaminant in cryogenic spacecraft, which
then turns into thin ice films. A notable example is the Gaia mission that
suffered from ice contamination (see Gaia Collaboration et al. (2016, Section
4.2.1)) and required several decontamination procedures to slowly remove the
ice from the optical system. Euclid Collaboration et al. (2023) studied the
ice formation and contamination for Euclid. The article also reviews the
lessons learned from other spacecraft on the topic of material outgassing. A
companion paper of Euclid Collaboration et al. (2023) is expected to be
published soon addressing the quantification of iced optics impact on Euclid’s
data.

  * •

Chromatic optical components: These components have a particular wavelength
dependence, excluding the natural chromaticity due to diffraction. They are
usually spectral filters and depend on the optical’s system design. A
particular example is a dichroic filter which ideally serves as an ideal band-
pass filter. The Euclid optical system includes a dichroic filter which allows
using both instruments, VIS and NISP, simultaneously as their passbands are
disjoint. A dichroic filter is made of a stack of thin coatings of specific
materials and thicknesses. Even if these components have a high-quality
manufacturing process, they can induce significant chromatic variations in
reflection affecting the PSF morphology. There is an ongoing effort in the
Euclid consortium to characterise the chromatic variations introduced by
Euclid’s dichroic filter. Baron et al. (2022) proposed a test bench to
characterise Euclid’s dichroic filter and a numerical model of its chromatic
dependence.

  * •

Light polarisation: In section 4, we studied the scalar diffraction theory,
thus neglecting light polarisation. First, the optical system can induce
polarisation even when the incoming light is not polarised. Breckinridge et
al. (2015) studied the effect of polarisation aberrations on the PSF of
astronomical telescopes. The study of polarisation is carried out using Jones
matrices (Jones, 1941). These matrices describe a ray’s polarisation change
when going through an optical system. See McGuire and Chipman (1990, 1991);
Yun et al. (2011) for more information on polarisation aberrations. Second,
there are some regions in space where the incoming light has been polarised by
different sources, e.g., Galactic foreground dust. Lin et al. (2020) studied
the impact of light polarisation on weak lensing systematics for the Roman
space telescope (Spergel et al., 2015). The study found that the systematics
introduced by light polarisation are comparable to Roman’s requirements.

  * •

Thermal variations: The thermal variations in a telescope introduce mechanical
variations in its structure that affect the performance of the optical system.
The origin of the thermal variations is strong temperature gradients due to
the sun’s illumination. It is sometimes referred to as the telescope’s
breathing (Bély et al., 1993) for the periodical pattern consequence of its
orbit. Thermal variations can introduce a small defocusing of the system that
will change the PSF morphology. This phenomenon was first identified in the
HST (Hasan et al., 1993). Nino et al. (2007) studied HST focus variations with
temperature and Lallo et al. (2006) studied HST temporal optical behaviour,
where temperature variations play a principal role. Later works (Makidon et
al., 2006; Sahu et al., 2007; Suchkov and Casertano, 1997) studied the impact
of the thermal variations, and consequently PSF variations, on different
science applications. A Structural-Thermal-Optical Performance (STOP) test
helps predict thermal variations’ impact on the optical system. This effect is
naturally more significant in space-based telescopes as the temperature
gradients in space are considerably more prominent than the ones found on the
ground. Space-based telescopes located at the stable ${\rm L}_{2}$ Lagrange
point, e.g., Euclid and JWST, are less prone to thermal variations compared to
telescopes orbiting the Earth, e.g., HST.

As an example, Figures 7() and 7() show the measured optical contribution for
the James Webb Space Telescope (JWST) PSF. Rigby et al. (2022) presents a
detailed analysis of JWST’s state, from its commissioning, including its PSF.

###  6.4 Detector-level degradations

The detector-level degradations are related to the detectors being used and,
therefore, to the intensity of the PSF. They affect the observed images
through the degradation operator $F_{p}$ from Equation 28, and as we will use
star images, or eventually other observations, to constrain PSF models, it is
necessary to consider their effects. Some of these degradations are non-
convolutional and will not be well-modelled by a convolutional kernel.
Nevertheless, we expect that image preprocessing steps will mainly correct
these effects. However, the correction will not be perfect, and some modelling
errors can propagate to the observations.

Figure 7: Example of two different pixellations on the same high-resolution
image representing an Airy PSF. The difference between the two pixellations is
an intra-pixel shift of $(\Delta x,\Delta y)=(0.35,0.15)$ between them. Figure
reproduced from Liaudat (2022).

  * •

Undersampling and pixellation: The EM wave that arrives at the detectors is a
continuous function. The discrete pixels in the detectors integrate the
functions and measure the intensity of the wave in their respective area. We
name this process pixellation, also known as sampling. Some authors, e.g.,
Anderson and King (2000); Bernstein (2002); Kannawadi et al. (2016), define an
effective PSF as the convolution of the optical PSF, i.e., the flux
distribution at the focal plane from a point source, with the pixel response
of the instrument, e.g., a $2$D top-hat function. High et al. (2007) performed
an early study on the effects of pixellation in WL and the choice of pixel
scale for a WL space-based mission. Krist et al. (2011, Section 3) gives some
insight on pixellation effects for HST. Two aspects of pixellation play a
crucial role in PSF modelling. First, the sampling is done with the same grid,
but it is indispensable to consider that the continuous function is not
necessarily centred on the grid. This difference means that intra-pixel shifts
between the different pixellations will be found. Figure 7 shows how two pixel
representations of the same light profile change due to two different
pixellations. When optimising a PSF model to reproduce some observed stars,
the centroids of both images must be the same. Suppose the image centroids are
the same, and the underlying model represents the observations satisfactorily.
In that case, the residual image between the two pixellated images will be
close to zero. If the centroids are not the same, the residual can be far from
zero even though the model is a good representation of the observation, as
illustrated in the residual image in Figure 7. The second aspect is related to
the Nyquist-Shannon sampling theorem. The theorem states the required number
of samples we need to use to perfectly determine a signal of a given
bandwidth. In the telescopes we study, the bandwidth and number of samples are
related to the aperture’s diameter and pixel size. Depending on the
telescope’s design, the sampling may not verify the Nyquist-Shannon theorem.
If the images are undersampled, meaning that the theorem is not verified, a
super-resolution step is required in the PSF modelling, which is the case of
Euclid. Using an observation strategy with dithering, as described in
subsection 6.2, can significantly mitigate the undesired effects of
undersampling and pixellation. Kannawadi et al. (2021) studies ways of
mitigating the effects of undersampling in WL shear estimations using
metacalibration (Huff and Mandelbaum, 2017; Sheldon and Huff, 2017; Sheldon et
al., 2020), which is a method for measuring WL shear from well-sampled galaxy
images. Finner et al. (2023) studies near-IR weak-lensing (NIRWL) measurements
in the CANDELS fields from HST images. The authors find that the most
significant contributing systematic effect to WL measurements is caused by
undersampling.

  * •

Optical throughput and CCD Quantum Efficiency (QE): The optical throughput of
the system is the combined effect of the different elements composing the
optical system, such as mirrors and optical elements like coatings (Venancio
et al., 2016). The filter being used in the telescope forms part of the
optical throughput, as can be seen in Figure 6 for the MegaCam set of filters.
Figure 6 also includes the CCD QE, which describes the sensibility of the CCD
to detect photos of different wavelengths. Commonly, CCDs do not have a
uniform response to the different wavelengths. Therefore, we must multiply the
CCD QE with the telescope’s optical throughput to compute the total
transmission.

  * •

CCD misalignments: Ideally, we expect that all the CCDs in the detector lie in
a single plane that happens to be the focal plane of the optical system.
However, it is not the case in practice, as there might be small misalignments
between the CCDs introducing small defocuses that change from CCD to CCD. See
Jee and Tyson (2011, Figure 8) for a study of this effect for the Vera C.
Rubin Observatory.

  * •

Guiding errors: Even if space telescopes are expected to be very stable when
doing observations thanks to the attitude and orbit control systems (AOCS),
there will exist a small residual motion that is called pointing jitter. The
effect on the observation is introducing a small blur that can be modelled by
a specific convolutional kernel that depends on the pointing time series.
Fenech Conti (2017, Section 4.8.3) proposes to model the effect for Euclid
with a Gaussian kernel.

  * •

Charge Transfer Inefficiency (CTI): CCD detectors are in charge of converting
incoming photons to electrons and collecting them in a potential well in the
pixel during an exposure. The charge on each pixel is read when the exposure
finishes. The collected electrons are transferred through a chain of pixels to
the edge of the CCD, amplified and then read. High energy radiation above the
Earth’s atmosphere gradually damages CCD detector (Prod’homme et al., 2014b,
a). The silicon damage in the detectors creates traps for the electrons that
are delayed during the reading procedure. This effect is known as CTI,
producing a trailing of bright objects and blurring the image. This effect is
noticeably significant for space telescopes, given the harsh environment. CTI
effects are expected to be corrected in the VIS image preprocessing. Rhodes et
al. (2010) carried out a study on the impact of CTI on WL studies. Massey et
al. (2009) developed a model to correct for CTI for the HST and later improved
it in Massey et al. (2014).

  * •

Brighter-fatter effect (BFE): The assumption that each pixel photon count is
independent of its neighbours does not hold in practice. There is a
photoelectron redistribution in the pixels as a function of the number of
photoelectrons in each pixel. The BFE is due to the accumulation of charge in
the pixels’ potential wells and the build-up of a transverse electric field.
The effect is stronger for bright sources. Antilogus et al. (2014) studied the
effect and observed that the images from the CCDs do not scale linearly with
flux, so bright star sizes appear larger than fainter stars. Guyonnet et al.
(2015) and Coulton et al. (2018) proposed methods to model and correct this
effect. The preprocessing of VIS images is supposed to correct for the BFE,
but there might be some residuals.

  * •

Wavelength dependent sub-pixel response: There exist a charge diffusion
between neighbouring pixels in the CCD. Niemi et al. (2015) studied this
effect for an Euclid’s VIS CCD and modelled the response of the CCD. Niemi et
al. (2015) proposed to model the effect as a Gaussian convolutional kernel
where the standard deviations of the 2D kernel are wavelength dependent,
$\sigma_{x}(\lambda)$ and $\sigma_{y}(\lambda)$. Niemi et al. measured the
proposed model with a reference VIS CCD. Krist (2003) studied the charge
diffusion in HST and proposed spatially varying blur kernels to model the
effect.

  * •

Noise: There are several noise sources in the measurements. Thermal noise
(Nyquist, 1928) refers to the signal measured in the detector due to the
random thermal motion of electrons which is usually modelled as Gaussian.
Readout noise (Basden et al., 2004) refers to the uncertainty in the
photoelectron count due to imperfect electronics in the CCD. Dark-current shot
noise (Baer, 2006) refers to the random generation of electrons in the CCD,
and even though it is related to the temperature, it is not Gaussian. There
are also unresolved and undetected background sources that contribute to the
observation noise. Which are the statistics of the predominant noise will
depend on the imaging setting of the instrument and its properties.

  * •

Tree rings and edge distortions: There exist electric fields in the detector
that are transverse to the surface of the CCD. The origin of these fields
includes doping gradients or physical stresses on the silicon lattice. This
electric field displaces charge, modifying the effective pixel area.
Consequently, it changes the expected astrometric and photometric
measurements. This electric field also generates concentric rings, tree rings,
and bright stripes near the boundaries of the CCD, edge distortions. Given the
close relationship between this effect and the detector, its importance
depends strongly on the instrument being used. This effect is unnoticed in the
MegaCam used in Canada-France Imaging Survey (CFIS) as it depends on the CCD
design. However, it was a major concern in the Dark Energy Camera used in the
Dark Energy Survey (DES), as shown by Plazas et al. (2014). Jarvis et al.
(2020, Figure 9) illustrates the consequence of tree rings in the PSF
modelling.

  * •

Other effects: These effects include detector nonlinearity (Plazas et al.,
2016; Stubbs, 2014), sensor interpixel correlation (Lindstrand, 2019),
interpixel capacitance (McCullough, 2008; Kannawadi et al., 2016; Donlon et
al., 2018), charge-induced pixel shifts (Gruen et al., 2015), persistence
(Smith et al., 2008b, a), reciprocity failure (flux-dependent nonlinearity)
(Biesiadzinski et al., 2011; Bohlin et al., 2005), and detector analogue-to-
digital non-linearity.

###  6.5 The atmosphere

The atmosphere plays a central role in ground-based telescopes’ PSFs. See
Roddier (1981) for an in-depth study of the subject. How the atmosphere
affects our images will strongly depend on the exposure time used to image an
object. The PSF induced by the atmosphere for a very short exposure will look
like a speckle, while a long exposure will produce a PSF that resembles a 2D
Gaussian, or more precisely, a Moffat profile (Moffat, 1969). Figure 9 shows
examples of atmospheric PSFs with different exposure times. As a first
approximation for long exposures, the effect of the atmosphere on the PSF is
that of a spatially varying low-pass filter. Therefore, broadening the PSF and
limiting the telescope’s resolution. Astronomers usually use the term seeing
to refer to the atmospheric conditions of the telescope, and it is measured as
the FWHM of the PSF. The loss of resolution due to the atmosphere is one of
the main motivations for building space telescopes like Euclid and Roman,
where the PSF is close to the diffraction limit and very stable.

The atmosphere is a heterogeneous medium whose composition changes with the
three spatial dimensions and time. The inhomogeneity of the atmosphere affects
the propagation of light waves that arrive at the telescope. Instead of
supposing that the incoming light waves are plane, as emitted by the faraway
source under study, these waves already have some phase lags or leads with
respect to an ideal plane wave. The atmosphere introduces a WFE contribution
to the optical system. These effects can be resumed as an effective phase-
shifting plate, $\Phi_{\text{eff}}(x,y,t)$. However, calculating this
effective plate is cumbersome as it involves having a model of the atmosphere
and integrating the altitude, $z$, so that we have the spatial distribution,
$(x,y)$, of the effective WFEs. The model of the atmosphere is represented by
the continuous $C_{n}^{2}(z)$ (Roddier, 1981) profile, which represents the
variations of the refractive index due to the atmospheric turbulence as a
function of height. However, the $C_{n}^{2}(z)$ is challenging to model and
measure, and even if it is possible, it is computationally expensive to
exploit.

Figure 8: Illustration of six von Kármán phase screen layers at different
altitudes simulated for LSST. The simulations were produced with the GalSim
package (Rowe et al., 2015) using the parameters from Jee and Tyson (2011).

We can discretise the integral over the altitude into $M$ thin-phase screens
of variable strength at different altitudes to simulate the effect of the
atmosphere. Each phase screen will have specific properties and move at
different speeds in different directions. These assumptions are known as the
frozen flow hypothesis. Each phase screen will be characterised by its power
spectrum that can be modelled by a von Kármán model of turbulence (Kármán,
1930). The power spectrum of the atmosphere’s WFE contribution writes

$\Psi(\nu)=0.023\,r_{0}^{-5/3}\left(\nu^{2}+\frac{1}{L_{0}^{2}}\right)^{-11/6}\;,$
(32)

where $\nu$ is a spatial frequency, $r_{0}$ is the Fried parameter, and
$L_{0}$ is the outer scale. Both parameters, $r_{0}$ and $L_{0}$, are
generally expressed in meters. The Fried parameter relates to the turbulence
amplitude, and the outer scale relates to the correlation length. See Figure 8
for an example of atmospherical phase screens. For lengths longer than
$L_{0}$, the power of the turbulence asymptotically flattens. If we take the
limit of $L_{0}$ to infinite, we converge to the Kolmogorov model of
turbulence (Kolmogorov, 1991). See Sasiela (1994) for more information on
electromagnetic wave propagation in turbulence.

Once the phase screens, $\Phi_{m}(x,y|u_{i},v_{i})$, have been simulated
following Equation 32, the temporal variation of the screen has to be taken
into account. The phase screens contribute to the WFE of the PSF, which is why
it depends on the pupil plane variables $(x,y)$. The temporal variation is
usually modelled with the wind’s properties at the phase screen’s reference
altitude. We describe the wind with two components, $v_{u}$ and $v_{v}$, where
we have assumed that $v_{z}=0$. We then obtain the effective phase screen by a
weighted average of the phase screens at the different altitudes as

$\Phi_{\text{eff}}(x,y;t|u_{i},v_{i})=\sum_{m=1}^{M}c_{m}\Phi_{m}(x,y;t|u_{i},v_{i})\;,$
(33)

where $\\{c_{m}\\}$ are some weights. The difficulty of modelling the
atmosphere is that the time scales are comparable with the exposure time.
Therefore, the PSF we estimate for a given time snapshot will change with
respect to another PSF at another snapshot within the same camera exposure.
This change means that we need to integrate the instantaneous PSF over time to
model the PSF physically, which corresponds to

$I_{\text{img}}(\bar{u},\bar{v}|u_{i},v_{i})=\int_{t_{o}}^{t_{o}+T_{\text{exp}}}I_{\text{img}}(\bar{u},\bar{v};t|u_{i},v_{i})\;dt\;,$
(34)

where $I_{\text{img}}(\bar{u},\bar{v},t|u_{i},v_{i})$ is the instantaneous
image of the object affected by the PSF
$\mathcal{H}(u,v;\lambda;t|u_{i},v_{i})$, $t_{o}$ is a random initial time and
$T_{\text{exp}}$ is the exposure time.

Finally, we need to choose the time step size to discretise the integral from
Equation 34. Each instantaneous PSF will look like a speckle. Once we add them
up in the integral, the PSF starts becoming rounder and smoother. Figure 9
shows examples of atmospheric PSFs using different exposure times that were
simulated using $6$ phase screens using the parameters from Jee and Tyson
(2011) that correspond to an LSST-like scenario. It is interesting to see how
the short-exposure PSF looks like a speckle, and then the profile becomes more
and more smooth as the exposure time increases. As a reference, the exposure
time used for the r-band observations in CFIS is $200$s (science
collaboration, 2016). de Vries et al. (2007) studied the PSF ellipticity
change due to atmospheric turbulences as a function of the exposure time. de
Vries et al. observed that the ellipticity of the PSF decreases its amplitude
as the exposure time increases.

Figure 9: Example of atmospheric PSFs with different exposure times. The
simulation was done using the atmospherical parameters from Jee and Tyson
(2011) for an LSST-like scenario.

Another effect that should be considered is atmospheric differential chromatic
refraction. This effect represents the refraction due to the change of medium
from vacuum to the Earth’s atmosphere. The effect varies as a function of the
zenith angle and wavelength. Meyers and Burchat (2015b) performed a study on
the impact of the atmospheric chromatic effect on weak lensing for surveys
like LSST and DES.

Heymans et al. (2012) performed a study on the impact of atmospheric
distortions on weak-lensing measurement with real data from CFHT. Heymans et
al. characterised the ellipticity contribution of the atmosphere to the PSF
for different exposure times. To achieve this, they computed the two-point
correlation function of the residual PSF ellipticity between the observations
and a PSFEx-like PSF model (described in detail in section 7). Salmon et al.
(2009) studied the image quality and the observing environment at CFHT. (Xin
et al., 2018) carried out a study of the PSF and the variation of its width
with time and wavelength for the Sloan Digital Sky Survey (SDSS) (York et al.,
2000; Gunn et al., 1998). Jee and Tyson (2011) carried out a simulation study
to evaluate the impact of atmospheric turbulence on weak-lensing measurement
in LSST. Jee and Tyson (2011) used the atmospherical parameters from
(Ellerbroek, 2002) that were measured in the LSST site in Cerro Pachón, Chile.
There is an ongoing project at LSST DESC444https://github.com/LSSTDESC/psf-
weather-station to leverage atmospheric and weather information at or near the
observation site to produce realistically correlated wind and turbulence
parameters for atmospheric PSF simulations.

Another way to simulate the atmosphere and the PSF is to use a photon Monte
Carlo approach, known as photon shooting. This line of work was carried out in
Peterson et al. (2015, 2019, 2020) with a simulator available coined
PhoSim555https://bitbucket.org/phosim/phosim_release/wiki/Home that aims to
model LSST PSF. The method consists of sampling photons from astronomical
sources and simulating their interactions with their models of the atmosphere,
the optics and the detectors. PhoSim is a powerful simulation tool but is not
intended to be a model for estimating the PSF from observations. The
simulation software GalSim666https://github.com/GalSim-developers/GalSim Rowe
et al. (2015) also incorporates options to simulate atmospheric PSFs from
phase-screen exploiting the methods from (Peterson et al., 2015).

To conclude, we have seen that it is possible to develop a physical model of
the atmosphere based on the optical understanding we have from section 4 and
the studies of atmospheric turbulence of Kármán and Kolmogorov. However, this
approach has two inconveniences. First, the approach requires physical
measurements of the atmosphere at the telescope’s site, which is not always
available. Second, it is computationally expensive, as there is an altitude
and temporal integration, to handle varying atmospherical properties and reach
the exposure time, respectively. In practice, it is required to use long
exposure times to obtain deeper observations, meaning observing fainter
objects that are important for weak-lensing studies. This fact simplifies the
PSF modelling task as the long temporal integration smooths the PSF profile
and the PSF spatial variations over the FOV. Therefore, a data-driven approach
to modelling the PSF can offer a feasible and effective solution in this
scenario.

###  6.6 Adaptive optics

An alternative approach to work with ground-based observations affected by the
atmosphere is to use adaptive optic (AO) systems (Beckers, 1993). This
technology significantly improves the observation resolution in ground-based
telescopes that is severely limited due to the atmosphere, as we have seen in
subsection 6.5. An AO system tries to counteract the effect of the atmosphere
on the incoming wavefront by changing the shape of a deformable mirror. The
key components of the AO system are wavefront sensors (WFS), wavefront
reconstruction techniques and deformable mirrors, which operate together
inside a control loop. The WFS provide information about the incoming
wavefront and usually incorporate a phase-sensitive device. The wavefront
reconstruction has to compute a correction vector for the deformable mirrors
by estimating the incoming wavefront from the WFS information. The control
loop works in real-time sensing and modifying the deformable mirrors so that
the wavefront received by the underlying instrument is free from the optical
path differences introduced by the atmosphere. Davies and Kasper (2012)
provides a detailed review of AO systems for astronomy. The LSST that will
carry out weak lensing studies contains an AO system (Neill et al., 2014;
Thomas et al., 2016; Angeli et al., 2014). Exoplanet imaging studies have
greatly benefited from AO systems and impose strict requirements for these
systems. Some examples are the SPHERE instrument in the VLT (Beuzit, J.-L. et
al., 2019) and the Gemini Planet Imager (Graham et al., 2007; Macintosh et
al., 2014, 2018) in the Gemini South Telescope.

##  7 Current PSF models

Let us now discuss some of the most known PSF models, which can be divided
into two main families, parametric and non-parametric, also known as data-
driven.

###  7.1 Parametric PSF models

This family of PSF models is characterized by trying to build a physical
optical system model that aims to be as close as possible to the telescope.
Once the physical model is defined, a few parameters are estimated using star
observations. Such estimation, also called calibration, is required as some
events, like launch vibrations, ice contaminations and thermal variations,
introduce significant variations in the model. These events prevent a complete
on-ground characterization from being a successful model. Parametric models
are capable of handling chromatic variations of the PSF as well as complex
detector effects. Nevertheless, parametric models have only been developed for
space missions and are custom-made for a specific telescope. The parametric
model is compelling if the proposed PSF model matches the underlying PSF
field. However, if there are mismatches between both models, significant
errors can arise due to the rigidity of the parametric models. The difficulty
of building a physical model for the atmosphere, already discussed in
subsection 6.5, makes them impractical for ground-based telescopes.

The parametric model Tiny-Tim777github.com/spacetelescope/tinytim (Krist,
1993; Krist and Burrows, 1995; Krist et al., 2011) has been used to model the
PSF of the different instruments on board the HST. The Advanced Camera for
Surveys (ACS) in HST was used to image the Cosmic Evolution Survey (COSMOS),
which covers a $2$ $\text{deg}^{2}$ field that was used to create a widely
used space-based weak-lensing catalogue. The first WL shape catalogue used the
Tiny-Tim model (Leauthaud et al., 2007). Rhodes et al. (2007) studied the
stability of HST’s PSF, noticing a temporal change of focus in the images.
Besides the parametric model Tiny-Tim model, Anderson and King (2000)
developed the concept of effective PSF, which is the continuous PSF arriving
to the detectors, i.e., Equation 26, convolved with the pixel-response
function of the detector. Anderson and King (2000) proposed an algorithm to
model the effective PSF iteratively from observed stars. Anderson and King
(2006) continued the work on the effective PSF, adding some improvements and
detailed a model for the HST instruments ACS and Wide Field Camera (WFC).
Hoffmann and Anderson (2017) then carried out a comparison between Tiny-Tim
and the effective PSF approach for ACS/WFC. The study shows that the effective
PSF approach consistently outperforms the Tiny-Tim PSFs, exposing the
limitations of parametric modelling. Anderson (2016) describes the adoption of
the effective PSF approach applied to Wide Field Camera 3 (WFC3/IR)
observations, which are undersampled. The software
Photutils888github.com/astropy/photutils (Bradley et al., 2022) provides an
implementation of the effective PSF approach from Anderson and King (2000)
with the enhancements from Anderson (2016). Schrabback, T. et al. (2010) used
a data-driven PSF model based on PCA when studying the COSMOS field.

The early severe aberrations of HST’s optical system were an important driver
of phase retrieval algorithms. Several efforts to characterise HST were
published in an Applied Optics special issue (Breckinridge and Wood, 1993).
Solving the phase retrieval problem provides a reliable approach to
characterise the optical system accurately. Fienup (1993) studied new phase
retrieval algorithms and Fienup et al. (1993) present several results
characterising HST’s PSF.

Euclid’s VIS parametric model constitutes the primary approach for Euclid’s
PSF modelling. The model will soon be published and used to work with
observations from Euclid.

###  7.2 Data-driven (or non-parametric) PSF models

The data-driven PSF models, also known as non-parametric, only rely on the
observed stars to build the model in pixel space. They are blind to most of
the physics of the inverse problem. These models assume regularity in the
spatial variation of the PSF field across the FOV and usually differ in how
they exploit this regularity. Data-driven models can easily adapt to the
current state of the optical system. However, they have difficulties modelling
complex PSF shapes occurring in diffraction-limited settings. One limitation
shared by all the data-driven models is their sensitivity to the available
number of stars to constrain their estimation. A low star number implies that
there might be not enough stars to sample the spatial variation of the PSF.
When the number of stars in a FOV falls below some threshold, the model built
is usually considered unusable for WL purposes. This family of models has been
widely used for modelling ground-based telescope PSFs. Nevertheless, they are
not yet capable of successfully modelling the chromatic variations in addition
to the spatial variations and the super-resolution.

We proceed by describing several PSF models in chronological order. The first
models, described in more detail, were used to process real data from
different surveys, except for Resolved Component Analysis (Ngolè et al.,
2016). The latter models are worth mentioning but have yet to be used to
produce a WL shape catalogue with all the validation and testing it implies.

####  7.2.1 Shape interpolation

The first approach for PSF modelling consisted in estimating some parameters
of the PSF at the positions of interest. It was done for early studies in WL
and is closely related to the widely-employed galaxy shape measurement method
KSB (Kaiser et al., 1995). This precursor method only required the PSF’s
ellipticity and size at the positions of the galaxies. Therefore, a full-pixel
image of the PSF was unnecessary. Then, the KSB method can correct the galaxy
shape measurement for the effects of the PSF. The method to interpolate the
shape parameters to other positions is usually a polynomial interpolation. For
example, this was the case for the early WL study of CFHTLS (Fu, L. et al.,
2008). Gentile et al. (2013) reviewed the different interpolation methods and
studied their performance for WL studies. Viola et al. (2011) performed a
study showcasing different biases present in the KSB (Kaiser et al., 1995)
method. These biases are a consequence of problematic KSB assumptions: (i) KSB
gives a shear estimate per individual image and then takes an average, while
WL shear should be estimated from averaged galaxy images (ii) KSB assumes that
galaxy ellipticities are small but what it is small in WL is the change in
ellipticity due to the shear, and (iii) KSB gives an approximate PSF
correction that only holds in the limit of circular sources and does not
invert the convolution with the PSF. Recent WL studies no longer use this
approach. The WL studies have evolved to more sophisticated galaxy shape
measurement techniques that require a full pixel image of the PSF at the
position of galaxies.

####  7.2.2 Principal component analysis (PCA)

Principal component analysis is a widely-known method for multivariate data
analysis and dimensionality reduction. Let us start with a set of star
observations in $\mathbb{R}^{p^{2}}$ that we concatenate in a matrix
$\bar{\mathbf{I}}=[\bar{I}_{1},\ldots,\bar{I}_{n}]$. We have flattened the 2D
images into an array to simplify expressions. One would like to represent the
observations with $r$ components $\\{X_{i}\\}_{i=1}^{r}$ in
$\mathbb{R}^{p^{2}}$, where $r\geq n$, assuming that $p^{2}>n$. The PCA method
gives $r$ orthonormal components in $\mathbb{R}^{p^{2}}$ which define
directions in the $\mathbb{R}^{p^{2}}$ space where the variance of the dataset
$\bar{\mathbf{I}}$ is maximized.

If $n$ components are used to represent the observations, then the learned
components in the PCA procedure represent a basis of the subspace spanned by
the observations or the columns of $\bar{\mathbf{I}}$. The method can be
interpreted as a linear transformation to a new representation with orthogonal
components. As it is usual to observe regularity in the spatial variations of
the PSF, most of the dataset variability can be described with a few
components. Then, one can only use the first $r$ principal components and
achieve a dimensionality reduction of the observations. The dimensionality
reduction technique allows denoising the model as the observational noise
cannot be represented with $r$ components and only the PSF trends are well
described.

The PCA method was used to model the PSF for the SDSS (Lupton et al., 2001),
although it was referenced as the Karhunen-Loève transform. Jarvis and Jain
(2004) then proposed its use in a WL context. Jee et al. (2007) used PCA to
model the spatial and temporal variations of the HST PSF. Jee and Tyson (2011)
also used PCA to model the PSF in LSST simulations. HST’s COSMOS catalogue
(Schrabback, T. et al., 2010) used PCA to model the PSF. PCA showcased the
utility and robustness of data-driven methods and the importance of using a
pixel representation of the PSF and is the precursor of several of the
following models.

####  7.2.3 PSFEx

PSFEx999github.com/astromatic/psfex (Bertin, 2011) has been widely used in
astronomy for weak-lensing surveys, for example, DES year $1$ (Zuntz et al.,
2018), HSC (Mandelbaum et al., 2017), and CFIS (Guinot et al., 2022). It was
designed to work together with the SExtractor (Bertin, E. and Arnouts, S.,
1996) software which builds catalogues from astronomical images and measures
several properties of the observed stars. PSFEx models the variability of the
PSF in the FOV as a function of these measured properties. It builds
independent models for each CCD in the focal plane and works with
polychromatic observations. It cannot model the chromatic variations of the
PSF field. The model is based on a matrix factorisation scheme, where one
matrix represents PSF features and the other matrix the feature weights. Each
observed PSF is represented as a linear combination of PSF features. The
feature weights are defined as a polynomial law of the selected measured
properties. This choice allows having an easy interpolation framework for
target positions. In practice, the properties that are generally used are both
components of the PSF’s FOV position. The PSF features are shared by all the
observed PSFs and are learned in an optimisation problem. The PSF
reconstruction at a FOV position $(u_{i},v_{i})$ can be written as

$\bar{I}^{\text{PSFEx}}_{\text{star}}(\bar{u},\bar{v}|u_{i},v_{i})=F^{\text{PSFEx}}\underbrace{\left\\{\sum_{\begin{subarray}{c}p,q\geq
0\\\ p+q\leq
d\end{subarray}}u_{i}^{p}\,v_{i}^{q}\,S_{p,q}(\bar{u},\bar{v})\,+\,S_{0}(\bar{u},\bar{v})\right\\}}_{\bar{H}^{\text{PSFEx}}(\bar{u},\bar{v}|u_{i},v_{i})}\,,$
(35)

where
$\bar{I}^{\text{PSFEx}}_{\text{star},(\cdot|u_{i},v_{i})}\in\mathbb{R}^{p\times
p}$ is the PSFEx reconstruction of the observed star
$\bar{I}_{(\cdot|u_{i},v_{i})}$, $S_{p,q}\in\mathbb{R}^{P\times P}$ represents
the learned PSF features or eigenPSFs, $S_{0}\in\mathbb{R}^{P\times P}$
represents a first guess of the PSF, the polynomial law is defined to be of
degree $d$, and $F^{\text{PSFEx}}$ represents the degradations required to
match the model with the observations. The model’s PSF reconstruction is
represented by $\bar{H}^{\text{PSFEx}}$. The first guess can be computed as a
function of the median of all the observations. The dimensions $p$ and $P$
will be the same if no downsampling operation is included in
$F^{\text{PSFEx}}$.

The PSF features are learned in an optimisation problem that aims to minimise
the reconstruction error between the PSFEx model and the observations, which
reads

$\min_{\begin{subarray}{c}S_{p,q}\\\ \forall p,q\geq 0\,,\,p+q\leq
d\end{subarray}}\left\\{\sum_{i=1}^{n_{\text{obs}}}\left\|\frac{\bar{I}_{(\cdot|u_{i},v_{i})}-\bar{I}^{\text{PSFEx}}_{\text{star},(\cdot|u_{i},v_{i})}}{\hat{\sigma}_{i}}\right\|_{F}^{2}+\sum_{\begin{subarray}{c}p,q\geq
0\\\ p+q\leq
d\end{subarray}}\left\|T_{p,q}\,S_{p,q}\right\|_{F}^{2}\right\\}\,,$ (36)

where $\hat{\sigma}_{i}^{2}$ represent the estimated per-pixel variances,
$\bar{I}$ represents the noisy observations, and $\|\cdot\|_{F}$ the Frobenius
norm of a matrix. The second term in Equation 36 corresponds to a Tikhonov
regularisation, where $T_{p,q}$ represents some regularisation weights to
favour smoother PSF models. The PSF recovery at target positions is
straightforward. One needs to introduce new positions in the Equation 35 after
learning the PSF features $S_{p,q}$. The recovery at a new FOV position
$(u_{j},v_{j})$ simply writes

$\bar{H}^{\text{PSFEx}}(\bar{u},\bar{v}|u_{j},v_{j})=\sum_{\begin{subarray}{c}p,q\geq
0\\\ p+q\leq
d\end{subarray}}u_{j}^{p}\,v_{j}^{q}\,S_{p,q}(\bar{u},\bar{v})\,+\,S_{0}(\bar{u},\bar{v})\;,$
(37)

where $\bar{H}^{\text{PSFEx}}$ is the model’s PSF reconstruction, and $S_{0}$
and $S_{p,q}$ were learned during the training procedure.

####  7.2.4 Resolved component analysis (RCA)

RCA101010github.com/CosmoStat/rca (Ngolè et al., 2016) is a state-of-the-art
data-driven method designed for the space-based Euclid mission (Schmitz et
al., 2020). The model builds an independent model for each CCD, can handle
super-resolution, and, like PSFEx, is based on a matrix factorisation scheme.
However, there are three fundamental changes with respect to PSFEx. The first
difference is that, in RCA, the feature weights are defined as a further
matrix factorisation, and are also learned from the data. The feature weights
are constrained to be part of a dictionary111111In the signal processing
community, a dictionary is a collection of templates, or basic elements, used
to decompose a signal. built with different spatial variations based on the
harmonics of a fully connected undirected weighted graph. The graph is built
using the star positions as the nodes and a function of the inverse distance
between the stars to define the edge weights. The rationale of the graph-
harmonics dictionary is to capture localised spatial variations of the PSF
that occur in space-based missions exploiting the irregular structure of the
star positions. The RCA reconstruction of an observed star is then

$\displaystyle\bar{I}^{\text{RCA}}_{\text{star}}(\bar{u},\bar{v}|u_{i},v_{i})=F^{\text{RCA}}\left\\{\bar{H}^{\text{RCA}}(\bar{u},\bar{v}|u_{i},v_{i})\right\\},\text{
where}\quad$ (38)
$\displaystyle\bar{H}^{\text{RCA}}(\bar{u},\bar{v}|u_{i},v_{i})=\sum_{r=1}^{N_{\text{comp}}}S_{r}(\bar{u},\bar{v})\,A[r,i]=\sum_{r=1}^{N_{\text{comp}}}S_{r}(\bar{u},\bar{v})\,(\alpha
V^{\top})[r,i]\,,$

where
$\bar{I}^{\text{RCA}}_{\text{star},(\cdot|u_{i},v_{i})}\in\mathbb{R}^{p\times
p}$, $F^{\text{RCA}}$ corresponds to the degradation model of RCA, including
downsampling, intra-pixel shifts among others, $S_{r}\in\mathbb{R}^{D\,p\times
D\,p}$ corresponds to the data-driven feature, i.e., eigenPSF, $r$ from a
total of $N_{\text{comp}}$ features, $D$ is the upsampling factor in case a
super-resolution step is required, and $A[r,i]$ is the $r$-th feature weight
of the $i$-th star. The feature weight matrix $A$ is decomposed into a sparse
matrix $\alpha$ and a dictionary matrix $V^{\top}$ of graph-based spatial
variations; see Ngolè et al. (2016) for more information.

To regularise the inverse problem, RCA enforces a low-rank solution by fixing
$N_{\text{comp}}$, a positivity constraint on the modelled PSF, a denoising
strategy based on a sparsity constraint in the starlet (Starck et al., 2015)
domain, which is a wavelet representation basis, and a constraint to learn the
useful spatial variations from the graph-harmonics-based dictionary. The
optimisation problem that the RCA method targets is

$\displaystyle\min_{S_{k},\alpha_{k}}$
$\displaystyle\Bigg{\\{}\frac{1}{2}\sum_{i=1}^{n_{\text{obs}}}\left\lVert\bar{I}_{(\cdot|u_{i},v_{i})}-F^{\text{RCA}}\left\\{\bar{H}^{\text{RCA}}_{(\cdot|u_{i},v_{i})}\right\\}\right\rVert_{F}^{2}+\sum_{r=1}^{N_{\text{comp}}}\|w_{r}\odot\Phi
S_{r}\|_{1}+$ (39)
$\displaystyle+\iota_{+}\left(\bar{H}^{\text{RCA}}_{(\cdot|u_{i},v_{i})}\right)+\iota_{\Omega}(\alpha)\Bigg{\\}}\,\quad\text{s.t.
}\quad\bar{H}_{i}^{\text{RCA}}=\sum_{r=1}^{N_{\text{comp}}}S_{r}\,(\alpha
V^{\top})[r,i],$

where $w_{r}$ are weights, $\Phi$ represents a transformation allowing the
eigenPSFs to have a sparse representation, e.g., a wavelet transformation,
$\odot$ denotes the Hadamard product, $\iota_{+}$ is the indicator function of
the positive orthant, and $\iota_{\Omega}$ is the indicator function over a
set $\Omega$, which is defined as a set of sparse vectors and is used to
enforce sparsity on $\alpha$.

The second difference with respect to PSFEx corresponds to the regularisations
used in the objective function from Equation 39, which ends up as being non-
convex due to the matrix factorisation and non-smooth due to the
$\|\cdot\|_{1}$ constraint. The optimisation is solved through a block
coordinate descent, as it is a multi-convex problem, and proximal optimisation
algorithms that tackle the non-smooth subproblems (Beck and Teboulle, 2009;
Condat, 2013).

The third difference is handling the PSF recovery at a new position
$(u_{j},v_{j})$. The recovery is carried out by a radial basis function (RBF)
interpolation of the learned columns of the $A$ matrix, issuing a vector,
$\hat{\mathbf{a}}_{j}\in\mathbb{R}^{N_{\text{comp}}}$, see Schmitz et al.
(2020) for more details. This way, the spatial constraints encoded in the $A$
matrix are preserved when estimating the PSF at galaxy positions. The
interpolated feature weights $\hat{\mathbf{a}}_{j}$ can be introduced in the
Equation 38 formulation to generate the PSF at the new $j$ position.

The RCA model has yet to be used to generate a WL shape catalogue from real
observations. Liaudat et al. (2021a) showed that RCA is not robust enough to
handle real ground-based observations from CFIS as some CCDs exhibited
significant errors of the PSF shape.

####  7.2.5 The Multi-CCD PSF model (MCCD)

MCCD121212github.com/CosmoStat/mccd (Liaudat et al., 2021a) is a state-of-the-
art data-driven method originally designed for the ground-based CFIS from
CFHT. MCCD can model the full focal plane at once by incorporating the CCD
mosaic geometry into the PSF model. The MCCD model rationale is explained by
the limitations of other PSF models that build independent PSF models for
every CCD, e.g., RCA and PSFEx, are: (i) fundamentally limited in the
possible’s model complexity due to the lack of constraining power of a reduced
number of star observations, and (ii) the difficulty of modelling PSF spatial
variations spanning the entire focal plane, i.e., several CCDs, from
independently modelled CCDs. MCCD overcomes these issues by building a PSF
model containing two types of variations, global and CCD-specific. Both
variations are modelled by a matrix factorisation approach, building over the
success of PSFEx and RCA. The global features are shared between all CCDs, and
the local CCD-specific features aim to provide corrections for the global
features. The feature weights are defined as a combination of the polynomial
variations from PSFEx and the graph-based variations from RCA. The model’s
optimisation is more challenging than the previous models and is based on a
novel optimisation procedure based on iterative schemes involving proximal
algorithms (Parikh and Boyd, 2014).

The MCCD model has proven robust enough to handle real observations from CFIS
(Liaudat et al., 2021a), giving state-of-the-art results. MCCD has been
incorporated into the recent ShapePipe shape measurement pipeline (Farrens, S.
et al., 2022) originally designed to process the CFIS survey and generate a WL
shape catalogue. The first version of the shape catalogue (Guinot et al.,
2022), spanning $1700~{}\text{deg}^{2}$, from ShapePipe used PSFEx. However,
the next version, spanning $\sim 3500~{}\text{deg}^{2}$, used the MCCD PSF
model and will be released soon.

####  7.2.6 lensfit

lensfit (Miller et al., 2007; Kitching et al., 2008; Miller et al., 2013)
refers to a Bayesian galaxy shape measurement method for WL surveys. It also
includes a data-driven PSF model that will also be referenced as lensfit and
is sparsely described throughout the different publications involving the
shape measurement lensfit (Miller et al., 2013; Kuijken et al., 2015; Giblin,
Benjamin et al., 2021). This method has been used with real data to produce
the WL shape catalogues of CFHTLenS (Erben et al., 2013; Miller et al., 2013),
KiDS-+VIKING-450 (Wright, Angus H. et al., 2019), KiDS-450 (Hildebrandt et
al., 2016; Fenech Conti et al., 2017), KiDS-1000 (Giblin, Benjamin et al.,
2021), and VOICE (Fu et al., 2018). However, the code is not publicly
available.

This PSF model differs from the previous ones. PSFEx and RCA learn some
features or eigenPSFs that all the PSFs share. The lensfit model is fitted on
a pixel-by-pixel basis. Each pixel is represented as a polynomial model of
degree $d$ of the FOV positions. The lensfit model can use all the
observations in one exposure, meaning that it uses several CCDs at once. The
model uses the low-order polynomials, up to degree $n_{\text{c}}<d$, to be
fitted independently for each CCD and the rest of the monomials are fitted
using the observations from all the CCDs. This multi-CCD modelling is a
significant change with respect to previous methods that built independent
models for each CCD. The total number of coefficients per pixel is

$N_{\text{coeff}}=\frac{1}{2}\left((d+1)(d+2)+(N_{\text{CCD}}-1)(n_{\text{c}}+1)(n_{\text{c}}+2)\right)\;,$
(40)

where $N_{\text{CCD}}$ is the total number of CCDs in the camera, $d$
represents the degree of the polynomial varying in the full FOV, and $n_{c}$
the polynomial that is CCD-dependent. We can write the description of the
pixel $(\bar{u},\bar{v})$ of the PSF model for a FOV position $(u_{j},v_{j})$
in CCD $k$ as follows

$\bar{H}^{\text{lensfit}}(\bar{u},\bar{v}|u_{j},v_{j})=\sum_{\begin{subarray}{c}p,q\geq
0\\\ p+q\leq
n_{\text{c}}\end{subarray}}u_{j}^{p}\,v_{j}^{q}\,a_{(p,q),(\bar{u},\bar{v})}^{k}+\sum_{\begin{subarray}{c}p+q>n_{\text{c}}\\\
p+q\leq d\end{subarray}}u_{j}^{p}\,v_{j}^{q}\,b_{(p,q),(\bar{u},\bar{v})}\;,$
(41)

where $a_{(p,q),(\bar{u},\bar{v})}^{k}$ is the coefficient specific of CCD
$k$, pixel $(\bar{u},\bar{v})$, and polynomial $(p,q)$ to be fitted to the
observations. The coefficient $b_{(p,q),(\bar{u},\bar{v})}$ is shared by all
the CCDs.

One thing to notice in this approach is that as the modelling of the PSF is
done pixel-by-pixel, then every observation should share the same pixel grid
of the PSF. There is no guarantee that an observation will have its centroid
aligned with the chosen pixel grid. Therefore, the PSF model has to be aligned
with the observations. Other methods, like PSFEx and RCA, interpolate the
model to the observation’s centroids. However, lensfit interpolates all the
observations to the model’s pixel grid with a sinc function interpolation
which implies interpolating noisy images. This procedure is described in
Kuijken et al. (2015).

For the KiDS DR2 (Kuijken et al., 2015), the hyperparameters used by lensfit
are $n_{c}=1$, $d=3$, and $N_{\text{CCD}}=32$ (from CFHT’s OmegaCAM
instrument), where the images used belong to a $32\times 32$ pixel grid. When
fitting the model’s parameters, each star is given a weight that is a function
of its SNR with the following empirical formula

$w_{i}=\frac{s_{i}^{2}}{s_{i}^{2}+50^{2}}\;,$ (42)

where $s_{i}$ is the measured SNR of the star $i$.

####  7.2.7 PSFs In the Full Field-of-View (PIFF)

PIFF131313github.com/rmjarvis/Piff (Jarvis et al., 2020) is a recently
developed PSF model that was used for the DES year $3$ WL shape catalogue
(Gatti et al., 2021) replacing PSFEx that was used for the DES year $1$
release. The PIFF model targets the LSST survey. Some improvements of PIFF
with respect to PSFEx are the ability to use the full focal plane to build the
model and modelling the PSF in sky coordinates rather than pixel coordinates.
PIFF offers a modular and user-friendly design that will enable further
improvements. The change of modelling coordinates was motivated by the strong
tree ring detector effect observed in the DES instrument, Dark Energy Camera,
which introduces astrometric distortions that are easier to correct in sky
coordinates. Pixel coordinates refer to the coordinates defined on the pixel
grid of the instrument. In contrast, sky coordinates refer to the angles in
the celestial sphere, known as right ascension (RA) and declination (DEC). The
geometric transformations that allow going back and forth between the pixel
and sky coordinates are known as WCS.

Being a modular PSF modelling code, PIFF allows choosing between different
options for the PSF model and the interpolation method. For example, the model
can be an analytical profile like a Gaussian, a Moffat or a Kolmogorov
profile, or a more general non-parametric profile called PixelGrid. The
interpolation method can be a simple polynomial interpolation, K-nearest
neighbours method, a Gaussian process (also known as Kriging), or a Basis-
function polynomial interpolation. Let us clarify the difference between the
first and last mentioned interpolations. The simple polynomial interpolation
first fits the PSF model’s parameters $\mathbf{p}$ for each observed star.
Then, it fits the coefficients of a polynomial of the 2D star positions that
will later be used to interpolate. In the Basis-function polynomial
interpolation, the position polynomial’s interpolation coefficients are fitted
simultaneously with the model’s parameters using all the available stars (from
a single CCD or the entire exposure). If this last option is used with the
PixelGrid model, it comes closer to the approaches of PSFEx and RCA without
the specific characteristics of each model. We have only mentioned position
polynomials, but, as in PSFEx, the interpolation polynomial can be built on
any parameter of the PSF, as, for example, a colour parameter.

The current PIFF PSF model includes an outlier check after the algorithm has
converged. The outlier check is based on the chi-squared, $\chi^{2}$, pixel
residual error between the model and the observations. The model implements an
iterative refining approach which means that after the model has converged,
one (or more) outlier star(s) is(are) removed from the observations. A new
iteration then starts with the model being recomputed. Although this approach
effectively removes outlier stars not representative of the PSF (because they
are binary stars or have some contamination), it can be very computationally
demanding. The computing time increases linearly with the number of iterations
used, which might be problematic depending on the total area to analyse or the
available computing resources. We refer the reader to Jarvis et al. (2020) for
more details.

The DES year $3$ shape catalogue (Gatti et al., 2021) used the PIFF model. The
model was a PixelGrid with Basis-function polynomial interpolation using a
$3$rd order polynomial. Even if PIFF has the potential to build a model across
several CCD chips, in practice, each model was built independently for each
CCD.

####  7.2.8 WaveDiff

The WaveDiff141414github.com/CosmoStat/wf-psf (Liaudat et al., 2023, 2021b)
PSF model was recently developed targeting space telescopes, in particular the
Euclid mission (Laureijs et al., 2011). WaveDiff proposes a paradigm shift for
data-driven PSF modelling. Instead of building a data-driven model of the PSF
in the pixel space as the previous models, WaveDiff builds its model in the
wavefront error (WFE) space. This change relies on a differentiable optical
forward model that allows propagating the wavefront from the pupil plane to
the focal plane and then computing the pixel PSF. The model is based on two
components, a parametric WFE and a data-driven WFE. The parametric WFE can be
based on optical simulations, characterisations of the optical system, or
complementary measurements such as phase diversity observations. The
parametric part should aim to model very complex dependencies that cannot be
inferred from the degraded star observations. Liaudat et al. (2023) proposes a
parametric model built using fixed features, namely the Zernike polynomials
(Noll, 1976). For several reasons, the parametric WFE model obtained often
cannot accurately represent the observed PSF, e.g., the telescope changes over
time, there are errors in the parametric model built, and relevant effects
were not considered or neglected. Consequently, the data-driven WFE should be
able to correct the aforementioned mismatches. This data-driven part is based
on a matrix factorisation with spatial variations inspired from PSFEx and RCA.
It is crucial to model smooth variations that have a reliable generalisation
of the PSF to different positions in the FOV. An overview of the model is
presented in Figure 10.

Figure 10: Overview of the WaveDiff approach to model the PSF with a
differentiable optical forward model. Figure adapted from Liaudat et al.
(2023).

Estimating the model’s parameters in the WFE space is a challenging, ill-posed
inverse problem known as phase retrieval Fienup et al. (1993); Fienup (1993);
Shechtman et al. (2015), i.e., estimating a complex signal from intensity-only
observations. The optimisation problem is non-convex and non-smooth, the star
observations are not very informative, and there is no guarantee that the WFE
model’s structure can represent the underlying ground truth WFE. (Liaudat et
al., 2023) shows that under the aforementioned conditions, targeting the
estimation of the ground truth WFE is not the best option. The PSF model’s
objective is to have a good pixel representation of the PSF. It is, therefore,
possible to estimate a WFE manifold far away from the underlying WFE but very
close in the pixel space. The data-driven features, the basis of the WFE
manifold, are estimated with a stochastic gradient descent method widely used
for estimating neural network parameters.

The WaveDiff model can handle spatial variations, super-resolve the PSF, and
model chromatic variations thanks to the WFE formulation and the optical
forward model, which also considers more general degradations as in Equation
30. To the best of our knowledge and at the time of writing, this is the only
data-driven PSF model able to cope with the spectral variations of the PSF.
WaveDiff shows a breakthrough in performance for data-driven PSF model in a
simplified Euclid-like setting (Liaudat et al., 2023). WaveDiff is flexible
enough to be adapted to different space telescopes. The framework proposed
shows an exciting research direction for future data-driven PSF modelling. The
WaveDiff model has yet to be tested with real space-telescope observations and
needs to incorporate detector-level effects, which are more naturally modelled
in pixel space. We refer the reader to Liaudat et al. (2023) and Liaudat
(2022) for more details.

####  7.2.9 Other PSF models

  * •

Shapelets: Refregier (2003) proposed a framework to analyse images based on a
series of localised basis functions of different shapes named shapelets.
Images can then be decomposed using these basis functions. Refregier and Bacon
(2003) continued the work proposing the shapelet framework to be used for
building shear estimates and modelling the PSF. The PSF modelling consists of
decomposing the star images in the shapelet basis and then performing an
interpolation of the coefficients to positions of interest. Essentially, it is
an extension of the approach seen in shape interpolation. Expressing the image
in shapelet coefficients allows denoising the star images and provides an
easier framework for the galaxy-PSF deconvolution. However, capturing all the
PSF structures in a finite expansion over analytical functions is not always
possible, leading to lost information. Massey and Refregier (2005) extended
the framework from Cartesian to polar shapelets.

  * •

Moffatlets and Gaussianlets: Li et al. (2016) proposed two other series of
basis functions to decompose the PSF named Moffatlets and Gaussianlets. Li et
al. compared the PSF reconstruction using the aforementioned basis with a PCA-
based method on LSST-like images. Using analytical basis functions leads to
more denoised models, as expected. Nie et al. (2021) continued the approach
and forced the principal components being learned in the PCA-like algorithm to
be built using the Moffatlets basis. This choice avoids the principal
components of learning noise as the Moffatlets basis avoids it. Both analyses
are missing a comparison with a reference PSF model like PSFEx to have a
reference performance. In addition, the performance comparison is made at the
same position as the observed stars, so the model’s generalisation to other
positions, a fundamental task of the PSF model, is not studied.

  * •

Fourier-based methods: Zhang (2007) proposed a Fourier-based method for
directly measuring the cosmic shear taking into account the PSF, which was
further developed in several publications (Zhang, 2011; Zhang et al., 2015; Lu
et al., 2017; Zhang et al., 2019). The method is based on the quadrupole
moments of the galaxy images (described in detail in subsection 9.2) but is
measured in Fourier space. The handling of the PSF is also done in Fourier
space. Lu et al. (2017) explores different interpolation approaches for the
PSF in the aforementioned Fourier framework. The 2D power spectrum of the
observed PSFs is interpolated to target positions. The interpolation is done
pixel-by-pixel, and the best-performing method is a well-parametrised
polynomial interpolation. An advantage of the Fourier interpolation is that
the 2D power spectrum is automatically centred in Fourier space, simplifying
the handling of images with intra-pixel shifts or, what is the same, different
centroids. Another Fourier-based shear measurement method is Bayesian Fourier
Domain (BFD) (Bernstein and Armstrong, 2014; Bernstein et al., 2016) built on
the Bayesian formalism. However, it does not include a specific PSF model.

  * •

Optimal transport (OT)-based methods: There exist two approaches involving OT
(Peyré and Cuturi, 2019) to tackle the PSF modelling problem. Ngolè and Starck
(2017) proposes to use Wasserstein barycenters as a non-linear geometric-aware
interpolation procedure of a low-dimensional embedding representation of the
PSFs. Although elegant, the performance of the approach does not seem to
justify its computational burden. In the comparison method, an RBF
interpolation of the principal components obtained through PCA achieves a
similar performance. The performance of the PCA method is better in terms of
ellipticity but slightly worse in terms of pixel error and PSF FWHM. Schmitz
(2019) worked on a data-driven PSF model based on RCA that would be able to
model the chromatic variations of the PSF through the use of Wasserstein
barycenters that were previously developed in Schmitz et al. (2018). The OT-
based PSF model coined $\lambda$RCA was compared to RCA. The comparison showed
a lower pixel and size error for $\lambda$RCA, although the ellipticity error
was similar or better for RCA. This method assumes that the PSF’s chromatic
variation is smooth over all the passband. This assumption holds if the only
chromatic effect of the PSF is due to the diffraction phenomena, which
exhibits a smooth variation with the $1/\lambda$ dependence in the WFE that
was already presented in Equation 26. However, if this is not the case and
another non-smooth chromatic variation is present, currently occurring in
Euclid (Venancio et al., 2016; Baron et al., 2022), there is no
straightforward way to adapt the $\lambda$RCA model to account for it.

  * •

Wavefront approach: Soulez et al. (2016) proposed to model the propagation of
light through the mirrors of the optical system. The PSF modelling problem is
recast into a phase retrieval problem. The article is a proof-of-concept as
there are only qualitative results, and many PSF-modelling difficulties remain
unaddressed.

  * •

Exploit out-of-focus images: Some instruments, like the Dark Energy Camera
(DECam) (Flaugher et al., 2015), are equipped with wavefront sensors that are
helpful for focus, alignment and adaptive optic system (AOS) (Roddier and
Roddier, 1993; Roodman, 2010, 2012). The LSST camera (LSST Science
Collaboration et al., 2009) will also be equipped with wavefront sensors
(Manuel et al., 2010; Xin et al., 2015, 2016; Claver et al., 2012). Roodman et
al. (2014) proposed to use the data from the wavefront sensors to constrain
the optical contribution of the PSF. The work was continued by Davis et al.
(2016) that proposed a wavefront-based PSF model for the DECam instrument
using out-of-focus observations. The PSF model is based on Zernike polynomials
fitted to out-of-focus stars, also called doughnuts, that contain considerably
more wavefront information than in-focus stars. Then, a new fit is done for
each exposure based on the measured quadrupole moments of the in-focus star
images. It is not easy to understand at which point the quadrupole moments
constrain the proposed PSF model and at which point it is the base physical
wavefront measured from the out-of-focus images that are the only part driving
the performance of the model. Snyder et al. (2016) proposed using the AOS
measurements to characterise the atmospheric turbulence in terms of a Zernike
decomposition.

  * •

Deep learning approaches: A model coined PSF-NET was proposed by Jia et al.
(2020) and is based on two convolutional neural networks (CNNs) trained
jointly. One network has to transform high-resolution images into low-
resolution images, while the other has to do the opposite. The CNNs are
trained in a supervised way expecting that the first network will learn a PSF
manifold. However, it is unclear how the approach handles the spatial
variation of the PSF, and it has not been tested for WL purposes. Jia et al.
(2020) proposed another approach for PSF modelling based on denoising auto-
encoders, but the spatial variation of the PSF remains untackled. Another
approach is followed by Herbel et al. (2018), where the PSF profile is
modelled by a parametric function consisting of a base profile of two Moffat
profiles and several parametrised distortions to increase the expressivity. A
CNN is trained in a supervised manner to predict the parameters of the
parametric profile from a noisy star observation. The neural network provides
a good estimation of the parameters, but the spatial variation of the PSF is,
again, not addressed. Having a PSF model that can model the observations is
important. However, in PSF modelling for WL analysis, a crucial part is to
capture the spatial variations of the PSF and that the model outputs the PSF
at different positions in the FOV.

##  8 Desirable properties of PSF models

In the previous section, we reviewed some of the most relevant PSF models
developed so far. After studying many models, we can conclude on desirable PSF
model properties. The PSF model should:

  * (a)

Have an accurate modelling of the PSF light profile. This modelling is
essential for any target task, as the light profile is the convolutional
kernel for a given position. The smoothness and structure in the PSF profile
are a consequence of the PSF being the Fourier transform, in Fraunhofer’s
approximation, of a particular finite-length aperture that limits the
frequency content of the PSF. This frequency limitation prevents us from
having a Dirac distribution as a PSF, as it would require an infinite
frequency content to build it. One difficulty is accurately modelling the
PSF’s wings, or outer region, which is often below the noise level. In ground-
based telescopes, the effect of the atmosphere can be interpreted as a low-
pass filter for the PSF, smoothing the PSF light profile.

  * (b)

Produce noiseless estimations of the PSF. The presence of noise in the PSF
estimations is an issue for purely data-driven models, which sometimes tend to
overfit noisy observations. Some regularisations have to be introduced in the
PSF model parameter optimisation to avoid producing noisy PSFs. A seemingly
straightforward solution to this problem is to rely on PSF models based on
fixed basis functions like Shapelets or Moffatlets. These models will always
output denoised PSFs as their basis functions cannot reproduce the noise.
However, they will introduce modelling errors if they cannot accurately model
the observed PSF light profiles.

  * (c)

Capture the PSF field’s variations. It is often the case that a good quality
PSF is required at the position of a certain object where no direct
information of the PSF is available. The PSF model first needs to capture most
of the relevant information from the observations at other positions and
wavelengths. Then, the model exploits this information and predicts the PSF at
the required position and wavelength. The PSF model relies upon its
generalisation power as it requires to exploit the PSF field information from
other positions.

  * (d)

Be able to exploit the structure of the PSF field variations. This desired
property is related to the previous point (c). An exciting approach to obtain
a good generalisation power is to learn the structure of the PSF field
variations. This structure is a consequence of the physical properties of the
telescope’s optical system. The subsection below provides a physical
understanding of the PSF field structure, which imposes a certain smoothness
to the variations. The spatial variations are also structured due to the
atmosphere if we study the PSF field of a ground-based telescope. A data-
driven PSF model should use a low-complexity representation of the PSF field,
which would be able to learn its structure and spatial variations.

  * (e)

Handle the CCD discontinuities of the PSF field. The CCDs misalignments are a
source of discontinuity of the PSF’s spatial variations. The PSF field is
piece-wise continuous, and the borders delimiting the discontinuity are well
known as the geometry of the focal plane is accurately known. A
straightforward way to handle the discontinuities is to model the PSF for each
CCD independently, e.g., PSFEx. Although this is simple to implement, it
limits the number of stars available to constrain the PSF model. Another more
difficult but potentially more powerful approach is to build a PSF model for
the entire focal plane, accounting for the focal plane discontinuities, e.g.,
MCCD.

  * (f)

Be robust to outliers and contaminations of the star sample. Contaminations
can come from the selection of stars, the fact that the objects classified as
stars are good representations of the PSF and are not small galaxies or binary
stars (Kuntzer and Courbin, 2017). Outliers can come from imperfect image
preprocessing where detector effects, like CTI, have not been adequately
removed. In addition, the model should be robust to different observation
conditions such as spatial distributions of the observed stars, SNRs, and the
number of observed stars.

  * (g)

Work appropiately with the target task. The PSF model can be exploited
differently according to the target task’s objective, e.g., estimating a
deconvolved object or estimating some summary statistic of the deconvolved
image. The model should be developed with the task in mind, as each task might
be more or less susceptible to different kinds of errors.

  * (h)

Be fast. Upcoming surveys will process a vast amount of observations.
Consequently, they put significant pressure on the computing time of PSF
models as they need to cope with the data intake. The requirements in terms of
computing time can drive many design choices in a PSF model, preventing the
use of costly physical simulations.

Once the PSF model has been developed with all the aforementioned properties
in mind, it is essential to validate the model’s performance. The validation
should ensure that the expected performance of the model is achieved and help
to identify sources of problems and provide directions for the improvement of
the model. In the next section, we give an overview of validation methods for
PSF models.

##  9 Validation of PSF models

The validation of PSF models is a challenging problem. To derive a validation
method, it would be necessary to quantify the impact of PSF modelling errors
on the final objective of our analysis. We consider, as an example, a weak-
lensing-based cosmological analysis, where the objective is to derive
constraints on the parameter of the cosmological model under analysis. This
exercise is challenging, given the analysis’ complexity and the large data
volume. Nevertheless, with some simplifying assumptions, it was carried out to
set the PSF modelling requirements for the Euclid mission as shown in
subsection 9.3. In this analysis, some assumptions on the PSF shape used do
not always hold for the high complexity of the PSF in a space-based mission
like Euclid. Even though it is essential to derive requirements for the PSF
model, these do not give much information on the nature of the errors and
possible problems the PSF model has. Therefore, it is necessary to derive
different diagnoses or null tests. Jarvis et al. (2016) proposed a set of null
tests for the DES WL shear catalogues science verification, including the PSF
model validation.

The most basic rule for any validation of PSF models is to separate the
observations in the FOV into two datasets for estimation and testing, i.e.,
validation, which could be $80\%$ and $20\%$, respectively. The first one
should be used to estimate the PSF model. The second one should help validate
its performance and not be used in the PSF model estimation. This rule tests
the PSF model’s generalisation power to unseen positions in the FOV. Next, we
will describe the most used PSF diagnosis that will help us to validate the
performance of the PSF models.

###  9.1 Pixel-based metrics

The most straightforward diagnostic we can think of is to compute the pixel
residual of our PSF model. Once trained, the model is used to recover the PSF
at test positions. We can then compute the RMSE of the pixel reconstruction
residuals. The PSF model can ideally predict the observed test stars without
error, and the reconstruction residual would only contain the observational
noise. If we work with simulations, we can produce noiseless stars for our
testing set, and the RMSE will directly indicate the pixel reconstruction
error. Even though pixel-based metrics can give insight into the PSF model
performance, they are not easily interpretable regarding scientific impact.
Errors in the PSF core or the PSF wings can impact the observed galaxy’s
estimated shape differently. With the existing methodology, it is difficult to
translate a pixel-based metric into a scientifically meaningful quantity in
terms of error propagation.

When working with real data, the PSF model’s validation with pixel-based
metrics becomes more complicated. The different noise levels in the data can
hide the pixel reconstruction error, making it difficult to compare different
PSF models or even assess the performance of a single one. Liaudat et al.
(2021a) proposed pixel-based reconstruction metrics for real observations. Let
us denote with
$I_{\text{star}}(\bar{u},\bar{v}|u_{i},v_{i}),I_{\text{model}}(\bar{u},\bar{v}|u_{i},v_{i})\in\mathbb{R}^{p\times
p}$ a test star and the predicted PSFs, respectively, at the FOV position
$(u_{i},v_{i})$, where $p^{2}$ is the total number of pixels in the image. To
simplify notation we write $I(\bar{u},\bar{v}|u_{i},v_{i})=I_{i}$. In star
images, most of the PSF flux is concentrated in the centre of the square
image, and the noise level can be considered constant in the image. Therefore,
we can mask the image to only consider the central pixels within a given