# Practical considerations for high-fidelity wavefront shaping experiments

Bahareh Mastiani  Biomedical Photonic Imaging Group, Faculty of Science and
Technology, University of Twente, P.O. Box 217, 7500 AE Enschede, The
Netherlands These authors contributed equally to this work.  e-mail:
<EMAIL_ADDRESS>Daniël W. S. Cox  Biomedical Photonic Imaging
Group, Faculty of Science and Technology, University of Twente, P.O. Box 217,
7500 AE Enschede, The Netherlands These authors contributed equally to this
work.  e-mail<EMAIL_ADDRESS>Ivo M. Vellekoop  Biomedical Photonic
Imaging Group, Faculty of Science and Technology, University of Twente, P.O.
Box 217, 7500 AE Enschede, The Netherlands

(February 2024)

###### Abstract

Wavefront shaping is a technique for directing light through turbid media. The
theoretical aspects of wavefront shaping are well understood, and under near-
ideal experimental conditions, accurate predictions for the expected signal
enhancement can be given. In practice, however, there are many experimental
factors that negatively affect the outcome of the experiment. Here, we present
a comprehensive overview of these experimental factors, including the effect
of sample scattering properties, noise, and response of the spatial light
modulator. We present simple means to identify experimental imperfections and
to minimize their negative effect on the outcome of the experiment. This paper
is accompanied by Python code for automatically quantifying experimental
problems using the OpenWFS framework for running and simulating wavefront
shaping experiments.

##  1 Introduction

The ability to focus light through scattering media is critical for a wide
range of applications, including deep tissue microscopy. However, it is
challenging to focus light deep inside opaque media due to inhomogeneities in
the refractive index of the medium. These inhomogeneities cause aberrations
and scattering that prevent light from forming a diffraction-limited focus,
thereby reducing the imaging resolution and contrast [16].

Figure 1: Schematic illustration of wavefront shaping. SLM: spatial light
modulator. OBJ: objective. a. Focusing inside a non-aberrating sample. b.
Focusing inside an aberrating sample results in a distorted focus. c. Focusing
inside an aberrating sample with a WFS microscope. By shaping the incoming
wavefront with a spatial light modulator, the focus can be restored.

Wavefront shaping (WFS) is a powerful method for focusing light at target
positions inside or behind scattering materials. The basic idea of WFS for
microscopy is illustrated in Fig. 1. Here, it is assumed that we use a
scanning microscope, and we are interested in focusing the light inside the
sample. In the ideal case (Fig. 1a), the sample is transparent and the
microscope objective focuses the light at the desired spot. We can describe
this focusing as constructive interference, where each part of the converging
wavefront interferes constructively at the focal point. When trying to focus
light into an aberrating or scattering sample, (Fig. 1b), light coming from
different parts of the incident wave picks up different random phase
deviations, causing random interference rather than constructive interference.
As a result, the intensity in the focus is reduced, and the focus is spread
out spatially. In the case of strong scattering, a disordered speckle pattern
is formed.

Wavefront shaping works by spatially modulating the incident wavefront to
exactly compensate these random phase deviations, restoring constructive
interference, causing the light to form a high-contrast, high-intensity focus
inside the scattering sample (Fig. 1c). Since the first demonstration of
focusing light through opaque scattering objects [38], the field has seen
rapid development towards applications in fundamental research, microscopy,
endoscopy, communication, and cryptography [15, 26, 8, 2, 11].

The key goal in WFS is to find the unique incident wavefront that maximizes
constructive interference in the target focus. A wide range of techniques are
already available (see e.g. Refs. [15, 12]), and exciting new approaches based
on digital twins [31, 44, 24, 46] or deep learning [8] are actively being
developed.

One can easily compute the theoretical expected improvement of the focus
quality under ideal circumstances [38]. In an actual experiment, however, the
measured improvement may be significantly lower. Here, we address the
following questions:

  * •

How much can I expect WFS to improve the focus?

  * •

Why does my WFS experiment perform worse than expected?

  * •

What can I do about it?

This paper is supported by our Python package OpenWFS [34], which includes
code for _automated analysis and troubleshooting_ of WFS experiments. The
package is available through the standard PyPI package repository [35], and
documentation can be found at the Read the Docs website [36]. Throughout this
paper, we will refer to the relevant functions for detecting specific problems
in WFS experiments.

For most of the analysis, we will assume that a ‘classical’ WFS approach is
used to find the optimized wavefront, that is, a phase-stepping technique that
uses feedback from the desired focal point [33, 25, 30, 21, 20]. Although this
covers only a fraction of the approaches in use, these classical algorithms
have the advantage that they can be described analytically, allowing one to
prove the convergence, analyze the effect of noise, and compute the expected
improvement under a large range of circumstances. Moreover, many of the
considerations here are generally valid regardless of the WFS algorithm used
to find the optimal wavefront.

After introducing the basic concepts that are needed to understand, run,
quantify, and predict the outcome of a WFS experiment, we cover the most
common issues influencing the outcome. We will explain how each of these
issues affect the enhancement, how we detect or quantify them, and how to
reduce their impact. In addition, we will show that the effect of most of
these issues can be quantified from regular phase-stepping measurements
directly. This allows for an approach where troubleshooting metrics are
conveniently built into the WFS algorithm at no extra cost. In addition, we
provide a troubleshooter for detecting issues that do require specific
measurement procedures to quantify them. We explain how to run the
troubleshooter and analyze the troubleshooting metrics in section 10.

##  2 Theory

We can describe propagation of light through an arbitrary object with a
matrix:

$E_{b}=\sum_{a}t_{ba}E_{a}$ (1)

where $E_{a}$ is the incident field for mode $a$, and the summation runs over
all optical modes. $E_{b}$ is the field at point $b$ where we want to focus
the light. When focusing light through a scattering object, $t$ is the
transmission matrix. Note, however, that this formalism is valid for any
linear system, so $t_{ba}$ may also refer to propagation to a point inside a
sample, or reflection, or in fact any linear system. Also note that the
formalism does not depend on how we represent these modes (e.g. as a Hadamard
or Fourier basis). We do restrict ourselves to orthonormal sets of modes, so
that we have the simple expression
$P_{\mathrm{in}}=\sum_{a}\left|E_{a}\right|^{2}$ for the the total incident
power.

The aim of WFS is to find the incident field $E_{a}$ that maximizes the
intensity111In this paper, we use the term ‘intensity’ for the power in a
single mode. When using a continuous representation of the transmission matrix
(see Ref. [22]) $I_{b}\equiv\left|E_{b}\right|^{2}$ in mode $b$, without
changing the total incident power, $I$ is a true intensity with unit
$W/m^{2}$, and all sums are replaced by surface integrals.). The answer to
this optimization problem is given by the Cauchy-Schwartz inequality, which
states that:

$I_{b}=\left|\sum_{a}t_{ba}E_{a}\right|^{2}\leq\sum_{a}\left|t_{ba}\right|^{2}\sum_{a}\left|E_{a}\right|^{2}$
(2)

The first term on the right hand side is a property of the sample. It can be
understood in terms of a phase conjugation experiment: imagine we would place
a point source in $b$, and observe the fraction of the total emitted power
that propagates to all modes $a$ together. In this case, this fraction is
exactly given by $\sum_{a}\left|t_{ba}\right|^{2}$, which we will denote
$T_{b}$ from here forward.

For a given sample and incident power, the right hand side of this inequality
is fixed. It directly follows that at most a fraction $T_{b}$ of the incident
light can be focused at $I_{b}$. This maximum is reached when both sides of
Eq. (2) are equal, which is the case when $E_{a}\propto t^{*}_{ba}$, leading
to an optimal value of

$I_{b}=\sum_{a}\left|t_{ba}\right|^{2}\sum_{a}\left|E_{a}\right|^{2}=T_{b}P_{in}$
(3)

All wavefront shaping algorithms for focusing light, in one way or the other,
are aimed at determining $t_{ba}$ and subsequently displaying $t_{ba}^{*}$ on
the SLM.

##  3 Experiment

###  3.1 Setup

A typical setup for WFS experiments is shown in Fig. 2. In this setup, HeNe
laser light (Thorlabs, $2\text{\,}\mathrm{m}\mathrm{W}$, $\lambda$ =
$632.8\text{\,}\mathrm{n}\mathrm{m}$) is expanded by a beam expander (Thorlabs
BE15M-A), and then modulated by a phase-only spatial light modulator (SLM,
Hamamatsu LCoS-SLM, X13138 SLM) using a 50$\%$ beam splitter (Thorlabs
CM1-BS013). A polariser ($P_{1}$, Thorlabs LPVISE100-A) is used to rotate the
polarization to the linear polarization required by the SLM. A half wave plate
(HWP, Thorlabs WPH10M-633) is used to control the total power entering the
system.

A 4f system images the SLM onto the back focal plane of a microscope objective
(Zeiss A-Plan 100$\times$/0.8 M27), which focuses the light onto the sample.
After transmitting through the sample, the light is collected by another
objective lens (Zeiss A-Plan 100$\times$/0.8 M27), and recorded by a CMOS
camera (Basler acA640-750um). A 4f system images the back surface of the
sample onto a CMOS camera. This camera provides the feedback signal for WFS.

Note that there are many other options to provide the feedback signal. A
mechanism to obtain a feedback signal that is correlated with localized light
intensity is referred to as a ‘guidestar’. This guidestar can be a localized
physical object (e.g. a fluorescent particle), as well as a localized physical
phenomenon (e.g. multi-photon fluorescence or the photo-acoustic effect). [12]

Figure 2:  Schematic of a simple WFS setup. HWP: half wave plate, P:
polarizer, M: mirror, $\times$15: beam expander, SLM: spatial light modulator,
BS: 50$\%$ non-polarizing beam splitter, obj: objective lens, $L_{1}:f=200$
mm, $L_{2}:f=75$ mm, $L_{3}:f=75$ mm, $L_{4}:f=150$ mm, $L_{5}:f=50$ mm, CMOS:
complementary metal oxide semiconductor camera.

###  3.2 Phase-stepping interferometry

Each element of the transmission matrix $t_{ba}$ corresponds to an initially
unknown complex optical field that must be measured. The phase of optical
fields can typically not be measured directly. However, we can indirectly
measure the relative phase of the field with a technique known as _phase-
stepping_ or _phase-shifting_ interferometry [17].

We will explain the use of phase-stepping interferometry with the _Stepwise
Sequential Algorithm_ (SSA) as an example, but the technique is generally
applicable. For SSA, the SLM pixels are divided into segments. We iterate over
each segment $a$ to measure its contribution $M_{a}=t_{ba}E_{0}$ to the field
at feedback location $b$. $E_{0}$ is the field amplitude at the SLM surface,
which is assumed to be equal for all segments. If the SLM is illuminated
unevenly, we can include the illumination profile in the matrix element
$t_{ba}$ (c.f. Section 6.2). The rest of the segments are used as static
reference and produce the field $R$ at $b$.

In each iteration, we phase-modulate the segment of interest by a set of
phases $\phi_{p}$ to produce the field $M_{a}e^{i\phi_{p}}$ at $b$. In other
words, the relative phase between the fields is varied between $0$ (including)
and $2\pi$ (excluding) in $P\geq 3$ steps. Due to the interference between
light originating from the modulated segment $M_{a}e^{i\phi_{p}}$ and light
originating from the other segments $R$, the feedback intensity will respond
sinusoidally:

$I^{(p)}_{ba}=|M_{a}e^{i\phi_{p}}+R|^{2}=|M_{a}|^{2}+|R|^{2}+M_{a}^{*}Re^{-i\phi_{p}}+M_{a}R^{*}e^{i\phi_{p}}$
(4)

where

$\phi_{p}=\frac{p}{P}2\pi,\qquad p=1,...,P$ (5)

$I^{(p)}_{ba}$ corresponds to the intensity in mode $b$, measured when
applying phase $\phi_{p}$ to the modulated input mode $a$. Note that the terms
with $M_{a}R^{*}$ and $M_{a}^{*}R$ carry a factor $e^{\pm i\phi_{p}}$. By
combining the measurements $I^{(p)}_{ba}$ in a _Discrete Fourier Transform_
(DFT), these terms can be isolated:

$F_{k,a}\equiv\frac{1}{P}\sum_{p=1}^{P}I^{(p)}_{ba}e^{-ik\phi_{p}},$ (6)

where $k$ is an integer in the range:

$k=\begin{cases}-P/2,\>\ldots,\>P/2-1,&\text{ if $P$ is even}\\\
-(P-1)/2,\>\ldots,\>(P-1)/2,&\text{ if $P$ is odd}\end{cases}$ (7)

It follows that $F_{1,a}=M_{a}R^{*}=t_{ba}E_{0}R^{*}\propto t_{ba}$. So we can
find the optimal input wavefront as:

$E_{a}^{\text{opt}}\propto\left(F_{1,a}\right)^{*}$ (8)

In theory $I^{(p)}_{ba}$ follows a sinusoid with an offset, hence only
$F_{-1,a},F_{0,a}$ and $F_{1,a}$ are non-zero. However, in real measurements,
the values for $\left|k\right|>1$ are never exactly zero, due to noise in the
measurements and imperfect SLM response calibration (see Section 8.1).

##  4 Enhancement

One way to quantify the improvement of the focus, is the fraction of the
incident power that is focused into mode $b$, i.e., the fraction $S\equiv
I_{b}/P_{in}$. In adaptive optics, it is usually assumed that all light is
scattered in the forward direction, so that $T_{b}=1$, and $S$ is called the
_Strehl ratio_ [16]. From Eq. (3), it is direct that in this case the optimal
wavefront achieves perfect diffraction-limited focusing: $S=1$.

For a scattering sample, it is often not possible to determine exactly what
fraction of the incident power reached the focus. Therefore, it makes more
sense to use a relative measure $\eta\equiv I_{b}/I_{0}$, which is called the
_enhancement_ [38]. The reference intensity $I_{0}$ is technically defined as
the intensity in point $b$, averaged over all possible ‘similar’ samples for a
given optimized incident wavefront [33], i. e.

$I_{0}=\left\langle\left|\sum_{a}\tilde{t}_{ba}E_{a}\right|^{2}\right\rangle,$
(9)

where $\tilde{t}_{ba}$ is the transmission matrix for a different sample with
the same scattering strength, thickness, etc., and
$\left\langle\cdot\right\rangle$ denotes averaging over all such samples.

For strongly scattering samples, $\tilde{t}_{ba}$ and $t_{ba}$ will be
completely uncorrelated. Moreover, we assume that all elements $t_{ba}$ are
statistically independent and identically distributed for all modes $a$ [9],
so that $\left\langle
t_{ba}^{*}t_{ba^{\prime}}\right\rangle=\delta_{aa^{\prime}}\tau$ where
$\delta_{aa^{\prime}}$ is the Kronecker delta, and the value
$\tau=\left\langle\left|t_{ba}\right|^{2}\right\rangle$ is an average ‘mode-
to-mode’ transmission coefficient which does not depend on mode index $a$. We
can now use these statistical properties and insert $E_{a}=E_{0}t_{ba}^{*}$
into Eq. (12) to find

$\displaystyle I_{0}$
$\displaystyle=E_{0}^{2}\left\langle\sum_{a,a^{\prime}}\tilde{t}^{*}_{ba}\tilde{t}_{ba^{\prime}}t^{*}_{ba}t_{ba^{\prime}}\right\rangle$
(10)
$\displaystyle=E_{0}^{2}\sum_{a,a^{\prime}}\delta_{aa^{\prime}}\tau\;t^{*}_{ba}t_{ba^{\prime}}$
(11) $\displaystyle=E_{0}^{2}\tau\sum_{a}\left|t_{ba}\right|^{2}=\tau
P_{\text{in}}$ (12)

In practice, $I_{0}$ can be determined by translating the sample, while
keeping the shaped wavefront constant, and averaging measurements taken at
different positions. Alternatively, the average speckle background can often
be used [33]. Expressions for the expected enhancement as a function of the
number of modes, and for phase-only light modulation are given in Section 6.1
and Section 6.2, respectively.

##  5 Fidelity

In an experimental setting, the generated wavefront always deviates from the
optimal one, so that the fraction of light that reaches the focus is always
lower than $T_{b}$. To account for these imperfections, we describe the
incident field as a superposition of the optimal field and a residual field
that is orthogonal to the optimal one [40]:

$E_{a}=\gamma
E_{a}^{\text{opt}}+\sqrt{1-\left|\gamma\right|^{2}}E_{a}^{\text{res}}.$ (13)

Here, the incident fields are normalized such that
$\sum_{a}\left|E_{a}^{\text{opt}}\right|^{2}=\sum_{a}\left|E_{a}^{\text{res}}\right|^{2}$.
The factor $\gamma$ denotes the fraction of the incident field that is shaped
correctly. It is defined as the correlation coefficient of the ideal wavefront
$E_{a}^{\text{opt}}$ and the actual wavefront $E_{a}$ [40]:

$\gamma\equiv\frac{\sum_{a}\left(E_{a}^{\text{opt}}\right)^{*}E_{a}}{\sqrt{\sum_{a}\left|E_{a}\right|^{2}}\sqrt{\sum_{a}\left|E_{a}^{\text{opt}}\right|^{2}}}.$
(14)

By definition, the residual field does not contribute to the target intensity
at all. Therefore, the fraction of light that reaches the focus equals
$\left|\gamma\right|^{2}T_{b}$. The factor $\left|\gamma\right|^{2}$ is called
the _fidelity_ , and equals the fraction of the incident power that is shaped
correctly. For ideal control over the incident wavefront
$\left|\gamma\right|^{2}=1$. In practice, however, there are many effects that
reduce the fidelity. In this paper, we will treat the following effects:

  * •

The effect of a finite number of modes $\left|\gamma_{N}\right|^{2}$

  * •

The effect of phase-only modulation $\left|\gamma_{A}\right|^{2}$

  * •

The effect of non-modulated light $\left|\gamma_{0}\right|^{2}$

  * •

The effect of imperfect phase modulation $\left|\gamma_{\phi}\right|^{2}$

  * •

The effect of phase jitter $\left|\gamma_{J}\right|^{2}$

  * •

The effect of a finite signal-to-noise ratio $\left|\gamma_{S}\right|^{2}$

  * •

The effect of drift and sample decorrelation $\left|\gamma_{D}\right|^{2}$

For each of these values, we will estimate the _expected_ value, averaged over
noise and over all ‘statistically similar’ samples; that is, samples that are
macroscopically the same, but differ in the microscopic placement of the
scattering particles. This averaging is denoted by
$\left\langle\cdot\right\rangle$.

We assume that each of these effects is uncorrelated with the other, meaning
that we can multiply the fidelities and simply estimate the expected intensity
in the target focus as

$\left\langle
I_{b}\right\rangle=\left\langle\left|\gamma_{D}\right|^{2}\right\rangle\left\langle\left|\gamma_{0}\right|^{2}\right\rangle\left\langle\left|\gamma_{\phi}\right|^{2}\right\rangle\left\langle\left|\gamma_{J}\right|^{2}\right\rangle\left\langle\left|\gamma_{S}\right|^{2}\right\rangle\left\langle\left|\gamma_{A}\right|^{2}\left|\gamma_{N}\right|^{2}T_{b}\right\rangle
P_{in}.$ (15)

Here, all factors were assumed to act independently, except for the terms
$\left|\gamma_{A}\right|^{2}\left|\gamma_{N}\right|^{2}T_{b}$, which were
averaged together, as described in Sections 6.1 and 6.2. Clearly, it is
essential that all fidelities are as close to $1$ as possible. For each of
these factors, we now describe how to estimate it from the experimental data,
and what measures can be taken to increase it if it is too far below $1$.

##  6 Limited control

###  6.1 Finite number of modes

So far, we considered $\sum_{a}$ a summation over _all_ optical modes. In
practice however, a WFS experiment only controls a fraction of the incident
angles and positions, and often only a single polarization. We now analyze how
this limitation affects WFS fidelity.

First, consider the case of _ideal phase and amplitude modulation_ , where the
SLM generates the optimal field $E_{a}=E_{0}t^{*}_{ba}$ for a finite number of
modes $a\leq N$, and that $E_{a}=0$ for $a>N$. Inserting this field into Eq.
(14) gives

$\left|\gamma_{N}\right|^{2}=\frac{\left[\sum_{a}^{N}E_{0}\left|t_{ba}\right|^{2}\right]^{2}}{\sum_{a}^{N}E_{0}^{2}\left|t_{ba}\right|^{2}\sum_{a}\left|t_{ba}\right|^{2}}=\frac{\sum^{N}_{a}\left|t_{ba}\right|^{2}}{\sum_{a}\left|t_{ba}\right|^{2}}.$
(16)

Unsurprisingly, this ratio equals $1$ when all optical modes are considered,
and it reduces when $N$ decreases. Intuitively, one can interpret
$\left|\gamma_{N}\right|^{2}T_{b}=\sum_{a}^{N}\left|t_{ba}\right|^{2}$ as the
fraction of the light that would propagate from a light source at point $b$
into the $N$ modes controlled by the SLM in a ‘time-reversal’ experiment.

For each additional mode, we need additional measurements to determine
$t_{ba}$. When fewer modes are considered, the WFS process will be faster, at
the cost of reducing the fidelity. In general, not all modes contribute
equally. In this case, it is beneficial to choose a set of modes that
contributes most to the intensity, and discard all other modes. For example,
for a sample that mainly scatters light in the forward direction, one can
choose a basis consisting of smooth modes that correspond to small angular
deviations [20]. For a homogeneously illuminated SLM, we can thus use the same
expression as for ideal phase and amplitude modulation (Eq. (16)).

A similar approach is taken in adaptive optics, where a set of low-order
Zernike modes may be used [16]. There is, however, a key difference: in
adaptive optics the modes are used to express the phase of the light rather
than the field. This results in a highly non-linear response for phase
variations $>1\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}$, causing many of the
algorithms used in adaptive optics to fail if aberrations are strong [27].

Figure 3: Measured beam location on the SLM. Each segment’s displayed
brightness corresponds to the measured contribution $\left|t_{ba}\right|^{2}$
for that segment.

In practice, the set of modes should be chosen correctly even for a strongly
scattering sample. For example, in the pupil-conjugate geometry shown in Fig.
2, the SLM is mapped to the back pupil of a microscope objective. If we choose
modes that correspond to (groups of) pixels, some modes may correspond to
pixels that are blocked by the back pupil of the objective (see Fig. 3).
Clearly, including these modes in the WFS process costs time and does not
contribute to the enhancement of the focus.

####  6.1.1 Enhancement in strongly scattering samples

For strongly scattering samples, we can use the assumptions described in
Section 2 to compute the estimated enhancement. Firstly, we use the fidelity
in Eq. (16) to find the optimized intensity of

$I_{b}=\left|\gamma_{N}\right|^{2}T_{b}P_{\text{in}}=\sum_{a}^{N}\left|t_{ba}\right|^{2}P_{\text{in}}.$
(17)

Together with the reference intensity $I_{0}=\tau P_{\text{in}}$ (Eq. (12)),
we get an enhancement of

$\eta\equiv\frac{I_{b}}{I_{0}}=\frac{\sum_{a}^{N}\left|t_{ba}\right|^{2}}{\tau},$
(18)

Which averages to

$\left\langle\eta\right\rangle=N.$ (19)

Therefore, for strongly scattering samples and ideal WFS, the expected
enhancement is equal to the number of segments. Due to the experimental
factors discussed below, in practice a lower value (around 0.5 N for strongly
scattering samples in our typical experiments) results.

####  6.1.2 Automatic detection and quantification

After running a WFS algorithm, we have access to the measured values of
$\left|t_{ba}\right|^{2}$ for each mode $a$ (see Section 3.2), which we call
the _contribution_ of that mode. By plotting the contribution for each mode,
as in Fig. 3, we can verify if all modes contribute approximately equally, or
whether some modes may be omitted. In practice, the modes that contribute will
often have a different contribution depending on their position on the SLM or
shape of the generated optical mode. In Section. 6.2, we explain how to
quantify the effect of this distribution automatically.

####  6.1.3 Mitigation

If the center and the size of the SLM pattern do not match the beam center and
size observed in the above measurement, these parameters must be manually
corrected. In OpenWFS the pattern can be translated and scaled with the
transform argument of the SLM object. In addition, the algorithm may be
written to skip modes that do not contribute strongly. Large gains can
selecting only the modes that contribute most, such as the modes in a Fourier
basis that correspond to small angle scattering [20].

When modes correspond to physical (groups of) pixels on a phase-only SLM, care
must be taken not to choose the pattern too small. If there is a significant
contribution of non-modulated light (from pixels outside the pattern), this
will also negatively affect the fidelity. See Section 7.1.

###  6.2 Phase only modulation

Often, a phase-only SLM is used in experiments. As a result, it is not
possible to match the amplitude distribution of the incident field to that of
$t_{ba}$, resulting in a reduction of the fidelity. We now analyze how this
effect reduces the fidelity beyond the reduction caused by selecting a finite
number of modes, $\left|\gamma_{N}\right|^{2}$ . For ideal phase-only light
modulation, we have an incident field
$E_{a}=E_{0}t_{ba}^{*}/\left|t_{ba}\right|$, in a basis where indices $a$
correspond to pixels on the SLM. Using the definition of the fidelity (Eq.
(14)) again, we find

$\left|\gamma_{A}\right|^{2}=\frac{\left[\sum_{a}^{N}E_{0}t_{ba}\;E_{0}t^{*}_{ba}/\left|t_{ba}\right|\right]^{2}}{\left(\sum_{a}E_{0}^{2}\right)\left(\sum_{a}^{N}E_{0}^{2}\left|t_{ba}\right|^{2}\right)}=\frac{\left[\sum^{N}_{a}\left|t_{ba}\right|\right]^{2}}{N\sum_{a}^{N}\left|t_{ba}\right|^{2}}=\frac{\overline{\left|t_{ba}\right|}^{2}}{\overline{\left|t_{ba}\right|^{2}}},$
(20)

where $\overline{\phantom{x}}$ denotes the averaging over all pixels of the
spatial light modulator.

Note that we assumed here that the SLM is illuminated homogeneously. However,
the illuminated beam on the modulator usually has Gaussian beam profile
resulting in a lower fidelity. We can simply include this illumination profile
in $t_{ba}$ so that $E_{0}$ is a constant, and we can still use Eq. (20).

####  6.2.1 Enhancement in strongly scattering samples

If all modes contribute equally on average, we can multiply the enhancement
$\eta$ in Eq. (18) by the fidelity in Eq. (20) to find the enhancement for
phase-only modulation

$\eta=\frac{\left[\sum_{a}^{N}\left|t_{ba}\right|\right]^{2}}{N\tau}.$ (21)

To compute the expected value for $\eta$ we use the fact that the elements
$t_{ba}$ for a strongly scattering sample follow a circular Gaussian
distribution [9]. This distribution has the property that
$\left\langle\left|t_{ba}\right|\right\rangle^{2}=\frac{\pi}{4}\left\langle\left|t_{ba}\right|^{2}\right\rangle=\frac{\pi}{4}\tau$.
Expanding the Eq. (25) , and ensemble averaging, we find

$\displaystyle\eta$
$\displaystyle=\frac{1}{N\tau}\left\langle\left[\sum_{a}^{N}\left|t_{ba}\right|\right]^{2}\right\rangle$
(22) $\displaystyle=\frac{1}{N\tau}\left(\left\langle\sum_{a,a^{\prime}\neq
a}^{N}\left|t_{ba}\right|^{2}\right\rangle+\left\langle\sum_{a}^{N}\left|t_{ba}\right|^{2}\right\rangle\right)$
(23) $\displaystyle=\frac{1}{N\tau}\left(\frac{\pi}{4}N(N-1)\tau+N\tau\right)$
(24)

giving the well known result [38]

$\left\langle\eta\right\rangle=\frac{\pi}{4}(N-1)+1.$ (25)

Note that the conditions leading to this result are almost never fulfilled
exactly. If the SLM is not illuminated evenly, or not all modes contribute
equally, we have to fall back to Eq. (20). For example, even if the sample is
completely random and $\left\langle\left|t_{ba}\right|^{2}\right\rangle$ does
not depend on $a$, we have to take into account the illumination profile on
the SLM, giving [32]

$\left\langle\left|\gamma_{A}\right|^{2}\right\rangle\equiv\frac{\overline{A_{a}}^{2}}{\overline{A_{a}^{2}}},$
(26)

with $A_{a}$ the field amplitude of the illumination at mode $a$.

Finally, note that for weak, forward scattering, we can approximate
$t_{ba}=e^{i\phi_{ba}}$. This means that there is no need to modulate the
intensity of the light and $\left|\gamma_{A}\right|^{2}=1$.

####  6.2.2 Automatic detection and quantification

This fidelity $\left|\gamma_{A}\right|^{2}$ in Eq. (20) is computed in OpenWFS
as the fidelity_amplitude attribute in the WFSResult object returned by a WFS
algorithm. If this value drops below $\pi/4\approx 0.79$ significantly, this
is an indication that not all modes contribute equally, which may be caused by
an uneven illumination of the SLM.

####  6.2.3 Mitigation

Figure 4:  SLM screen with a random wavefront displayed on the circular
geometry. Each segment’s displayed brightness corresponds to the applied
phase. The segment sizes are adjusted to compensate for the beam intensity
profile.

Ref. [21] pioneered the use of a circular geometry to mitigate the effect of a
non-uniform illumination of the SLM (see Fig. 4). The idea is to adjust the
segment size to the intensity of the illumination, using larger segments in
regions where the intensity is lower. This way, we are able to have each
segment, on average, contribute equally to the WFS procedure. For this
geometry, the segment size can be changed by adjusting the number of segments
per ring and the ring radii. In OpenWFS, this segment distribution is
available as the function slm.geometry.circular.

##  7 Imperfections of the spatial light modulator

In this section, we describe the typical limitations of liquid crystal on
silicon (LCoS) phase only SLMs, describe how these limitations affect the
fidelity, and provide strategies to mitigate these effects.

###  7.1 Non-modulated light

Figure 5:  The main transmitted/reflected beams in the simplified SLM layer
structure.

Most liquid-crystal SLM devices operate in reflection mode. As a result of the
layer structure of the SLM there are two main reflections. Fig. 5 shows the
simplified layer structure and the predominant transmitted/reflected beams.
The first reflection ($I_{1}$) is generated on the front surface of the
protective glass. The main beam ($I_{2}$), is reflected at the back surface,
thus traversing the liquid-crystal (LC) layer, and is the most intense
reflection. Finally, internal reflections inside the protective glass can
result in dim higher order reflections ($I_{3}$) [19].

The front surface reflection ($I_{1}$) of the SLM is not modulated and thus
does not contribute to the shaped wave. The reflection coefficient at the air-
glass interface depends on the SLM model and the wavelength of the incident
beam. For instance, the amplitude reflection coefficient at the air-glass
interface of the SLM, model X10468-08 from Hamamatsu, is measured as 0.198 for
$\lambda=$647\text{\,}\mathrm{n}\mathrm{m}$$ [19].

If we assume that the non-modulated light does not contribute to the intensity
in the focus, we can find the fidelity

$\left\langle\left|\gamma_{0}\right|^{2}\right\rangle=\frac{I_{\text{modulated}}}{I_{\text{total}}}.$
(27)

Note that for nearly transparent samples, the non-modulated light may still
contribute to the focus, so that this equation is not generally valid.
Finally, it should be noted that the popular parallel-aligned (PAL-LCoS) SLMs
only modulate one of the polarization states. Any component of the light in
the perpendicular polarization state contributes to non-modulated light.

####  7.1.1 Automatic detection and quantification

OpenWFS provides a simple means to estimate the amount of non-modulated light
originating from the front surface reflection of the SLM. It simply performs
phase stepping with a flat wavefront. This results in an intensity of

$I_{b}^{(p)}\propto|E_{0}\sqrt{1-\gamma_{0}^{2}}+E_{0}\gamma_{0}e^{i\phi}|+I_{\text{offset}},$
(28)

with $\gamma_{0}$ the fraction of the field that is not modulated, and
$I_{\text{offset}}$ the offset in the detector signal caused by non-coherent
sources, such as amplifier bias or ambient light. We first determine
$I_{\text{offset}}$ by closing the shutter of the laser, and subsequently
subtract this value from the signal. The amount of non-modulated light can now
be estimated from the phase-stepping data as

$\frac{\left|F_{1}\right|}{F_{0}}=\gamma_{0}\sqrt{1-\gamma_{0}^{2}}.$ (29)

Assuming that $\gamma_{0}^{2}>\frac{1}{2}$, we find

$\gamma_{0}^{2}=\frac{1+\sqrt{1-4\frac{\left|F_{1}\right|^{2}}{F_{0}^{2}}}}{2},$
(30)

which is output by the troubleshoot function of OpenWFS in the property
fidelity_non_modulated. Note that this function produces incorrect results in
the extreme case that more than half of the light is not modulated, such as
when the SLM is accidentally turned off.

####  7.1.2 Mitigation

High ratio’s of non-modulated light can be caused by a pattern of SLM segments
that underfills the illuminated SLM area, which can be solved by adjusting the
mapping of segments to the SLM surface (see Section 6.1.2).

A different source of non-modulated can be the previously explained front
reflection (see Fig. 5). The reflectivity of the surface depends on the model
of the SLM, and on the wavelength-specific anti-reflection coating applied to
it. If the front reflection is too strong, it is possible to tilt the SLM, and
aim the front reflection $I_{1}$ onto a beam dump. Of course this also
introduces tilt in the main reflection $I_{2}$, but this tilt can conveniently
be compensated for with a phase gradient on the SLM. Note that this technique
is typically unsuitable for pulsed light. The steep phase gradient typically
surpasses well beyond $2\pi$ and must therefore be phase-wrapped. Hence, the
pulse is delayed by different amounts on different parts of the SLM due to
phase wrapping.

Lastly, non-modulated light may be caused by the light not having the correct
polarization. PAL-SLMs only modulate one polarization component of the light.
By choosing an appropriate state of the incident polarization, we can minimize
the non-modulated light and thus optimize the fidelity.

###  7.2 Imperfect SLM calibration

In order to display an image on the SLM, the SLM hardware receives the image
and converts pixel values to voltages with the use of an internal lookup
table. Each SLM requires calibration to attain a linear phase response between
0 to 2$\pi$ for the wavelength of the incident laser light. Any deviation from
this linearity introduces a phase error for each input mode, and reduces the
fidelity of the WFS experiment. Therefore, it is crucial to calibrate the SLM
properly and load the correct lookup table to the SLM hardware. Alternatively,
OpenWFS provides the option to specify a software lookup table, which is
accelerated by the GPU using OpenGL.

Figure 6 shows an example of the measured phase response for the SLM
(Hamamatsu X13138-07) averaged over three trials. The measured phase response
closely matches the ideal linear phase response (dashed line) between 0 to
$2\pi$, which demonstrates that the proper lookup table is loaded on the SLM
for the chosen wavelength of incident laser light.

Figure 6:  Measured SLM phase response and the ideal phase response (dashed
line) versus the gray values displayed on the SLM.

If the response of the SLM deviates from a linear phase response, the fidelity
of displaying the wavefront reduces. To account for this effect, we assume
that the light modulator only modulates the phase, which is appropriate for
parallel-aligned (PAL) SLMs, but only approximately so for the older twisted-
nematic (TN) SLMs. Let $f(\phi_{p})$ denote the actual phase as a function of
the intended modulation phase $\phi_{p}$. The response of the light modulator
can now be written as a Fourier series

$e^{if(\phi_{p})}=\sum_{k}c_{k}e^{ik\phi_{p}},$ (31)

where the sum is over the range of $k$ as in Eq. (7). Since we assumed phase-
only modulation, the coefficient $c_{0}=0$. Note that this coefficient, if
present, would not effect the quantification measurements described below
anyway, and the effect of this non-modulated light is already accounted for in
$\left|\gamma_{0}\right|^{2}$, as explained in Section 7.1. The $c_{1}$
component corresponds to the intended phase response, and all other components
correspond to orthogonal modes that on average do not contribute to the focus.
We can thus quantify this effect with the fidelity

$\left\langle\left|\gamma_{\phi}\right|^{2}\right\rangle=\frac{\left|c_{1}\right|^{2}}{\sum_{k}\left|c_{k}\right|^{2}}.$
(32)

Note that for a perfect calibration $f(\phi)=\phi$, $c_{1}=1$ and all other
$c_{k}=0$. Also note that an imperfect phase calibration also reduces the
signal-to-noise ratio (SNR) in a measurement, causing an additional drop in
the fidelity (see Section 8.1).

####  7.2.1 Automatic detection and quantification

The effect of an imperfect calibration on the phase-stepping response Eq. (4)
can be written as:

$\displaystyle I^{pa}$
$\displaystyle=|M_{a}e^{if(\phi_{p})}+R|^{2}=|M_{a}|^{2}+|R|^{2}+M_{a}R^{*}e^{if(\phi_{p})}+M_{a}^{*}Re^{-if(\phi_{p})}$
(33)
$\displaystyle=|M_{a}|^{2}+|R|^{2}+M_{a}R^{*}\sum_{k}c_{k}e^{ik\phi_{p}}+M_{a}^{*}R\sum_{k}c_{k}^{*}e^{-ik\phi_{p}}$
(34)

The response is no longer a single sinusoid. This fact emerges in the form of
non-zero higher frequency components $F_{k,a}$ for $\left|k\right|>1$. Taking
the DFT as defined in Eq. (6) now yields:

$\displaystyle F_{k,a}=\begin{cases}|M_{a}|^{2}+|R|^{2}&\text{ for $k=0$}\\\
c_{k}M_{a}R^{*}+c_{-k}M_{a}^{*}R&\text{ otherwise}\end{cases}$ (35)

For phase stepping, we only use $F_{1,a}$. Note that if non-modulated light is
present, we can include it in the $R$ term and still choose $c_{0}=0$. Also
note that the signal is affected by the $c_{-1}$ component of the phase
response. The effect of a non-zero $c_{-1}$ component cannot be separated from
the measured response directly. Here, we assume that $c_{-1}\ll c_{1}$, so
that we can neglect it. In this case, $F_{1,a}\approx c_{1}M_{a}R^{*}$ and

$\displaystyle
F_{k,a}=\frac{c_{k}}{c_{1}}F_{1,a}+\frac{c_{-k}}{c_{1}^{*}}F^{*}_{1a}\qquad\text{
for $k\neq 0$}.$ (36)

To isolate the $c_{k}$ components, we perform a linear-least-squares fit. This
is done by constructing a matrix with two columns: $F_{1,a}$ and
$F_{1,a}^{*}$, and computing the Moore-Penrose pseudoinverse of this matrix:
$F^{+}$. The components of $c$ are then estimated from the measurements as

$\begin{bmatrix}c_{k}/c_{1}\\\
c_{-k}/c_{1}^{*}\end{bmatrix}=F^{+}\begin{bmatrix}F_{k,1}\\\ F_{k,2}\\\
\ldots\\\ F_{k,N}\end{bmatrix}$ (37)

giving a fidelity

$\left\langle\left|\gamma_{\phi}\right|^{2}\right\rangle=\frac{1}{\sum_{k\neq
0}{\left|c_{k}/c_{1}\right|^{2}}},$ (38)

which equals $1$ for a perfectly calibrated SLM. Computing
$\left\langle\left|\gamma_{\phi}\right|^{2}\right\rangle\leq 1$ provides a way
to check for systematic phase errors in the SLM response directly from the WFS
measurements. Note the pseudoinverse optimally uses the data of all
measurements, making it robust even in the presence of noise. Lastly, note
that if $c_{-1}$ is not negligible, some of the other components $c_{k}$ need
to be non-zero as well, in order to fulfill the phase-only constraint
$\left|e^{if(\phi)}\right|=1$. Therefore, we find a fidelity $<1$ in this case
too, so we can still detect a problem with the SLM calibration. The only
exception is when $c_{-1}=1$, meaning that the SLM modulates the phase in the
exact opposite direction we expect. This case cannot be detected with this
method. For feedback-based WFS this is not a problem at all, since the phase
sign flips occur both in the measurements and the construction of the focus
and thus cancel out. For model-based WFS, however, the sign of the phase
modulation is important and should be checked independently[31].

The value of $\left\langle\left|\gamma_{\phi}\right|^{2}\right\rangle$ is
computed automatically from the WFS measurements if the number of phase steps
is larger than three. It is available in the WFSResults structure as
fidelity_calibration.

####  7.2.2 Mitigation

Figure 7:  Schematic of the setup for calibrating the SLM lookup table. M:
mirror, BS: beam splitter.

A common method to calibrate the phase response of the SLM is using a
Michelson-Morley interferometer. Figure 7 shows a simplified schematic of the
calibration setup. The laser beam reflected from a mirror (M) interferes with
the modulated beam reflected from the SLM. The SLM is divided horizontally in
two halves. The phase of the top half (or the bottom half) is modulated from 0
to $2\pi$. This causes the fringe pattern in the top half to shift with
respect to the fringes in the bottom half. Since both top and bottom halves of
the beam follow the same path, this approach is insensitive to drift in the
path length of one of the arms of the interferometer. To define the phase
shift value in the fringes, we first perform Fourier transforms of the top and
bottom quarters of the fringes image. Neglecting the offset value (first
Fourier component), we find the Fourier component having the maximum amplitude
which corresponds to the spatial frequency of the fringes. The phase
difference of this maximum-amplitude component defines the phase of the
fringes. Finally, the measured phase response is inverted and loaded into the
SLM hardware as a lookup table. Lastly, note that several other calibration
procedures exist. [18]

###  7.3 SLM phase jitter

For LCoS SLMs, the effective phase applied by each pixel is not perfectly
stationary. Typically, the phase fluctuates periodically with a frequency
equal to the SLM refresh rate or a multiple thereof. This effect is known as
_phase jitter_ or _flicker_ [7], and it causes the displayed phase pattern to
deviate periodically from the phase pattern that was sent to the SLM. For a
uniformly distributed phase error between $-\Delta\phi$ and $\Delta\phi$, the
fidelity reduction is given by [3]

$\left\langle\left|\gamma_{J}\right|^{2}\right\rangle=\left(\frac{\sin\Delta\phi}{\Delta\phi}\right)^{2}.$
(39)

It can be seen that small errors in the phase do not contribute significantly
to the overall fidelity. However, this phase jitter also causes noise during
the measurements, which is especially prominent if the algorithm needs to
isolate a small variation on top of a large background of modulated light. The
phase jitter also adds noise to the WFS measurements, causing a further
reduction of $|\gamma|^{2}$ (see Section 8.1).

####  7.3.1 Automatic detection and quantification

The SLM phase jitter can be distinguished from regular noise by physically
turning off the SLM and observing the effect on the noise. The phase jitter
may be quantified by performing the interference calibration explained in
Section 7.2, but instead of modulating the phase of one of the halves, both
SLM halves remain constant and are observed at a high frame rate. Hence, for
this method the required observation frame rate is several times higher than
the refresh rate of the SLM.

####  7.3.2 Mitigation

The effect of phase jitter on the measurements can be reduced by setting the
exposure time of a camera, or the integration time of a detector, to a
multiple of the period of the jitter. This way, the jitter largely averages
out over the measurement. However, this does not reduce the effect of jitter
in constructing the wavefront, so a factor of
$\left\langle\left|\gamma_{J}\right|^{2}\right\rangle$ still arises there.

To reduce the amount of phase jitter, we may use the lowest possible voltages
on the SLM that still cause a full $2\pi$ phase shift, and choose zero grey
value for flat reference wavefronts whenever possible because the jitter is
typically higher at higher voltages. It may also be higher for specific gray
values due to the pulse code modulation scheme applied in many SLMs. The
amount of jitter is also dependent on the model and vendor of the SLM, and can
be further reduced by cooling the SLM [7].

###  7.4 SLM timing

An important aspect in optimizing SLM performance is correct timing of the
measurements. If the measurement is started too early, the liquid crystals are
not yet settled into their new state. If the measurement is started too late,
it wastes photon budget and needlessly increases the duration of the
experiment.

The image on the modulator is updated with the specified refresh rate. For
some SLMs, a hardware clock output is available, or may be isolated from the
video signal. OpenWFS, instead, automatically syncrhonizes the measurements
using a software triggering mechanism, where OpenGL is used to wait for the
_vertical retrace_ event, a historical name for the start of a new frame, and
update the image during that event. It is essential to only update the SLM
image during the vertical retrace: if it is updated before or after _tearing_
occurs, causing the SLM to display part of the old frame and part of the new
frame simultaneously.

After the vertical retrace, it takes some time for the image information to be
transferred from the PC to the SLM hardware, and it may also take some
additional time for the SLM hardware to start updating the voltages on the SLM
pixels. During this _latency_ , or _idle time_ , there is no optical response
of the SLM at all. We found that an idle time of 1-3 frame periods is typical,
depending on the SLM hardware.

The second interval, the settle time, is the time that is needed for the
liquid crystals of the SLM to stabilize. We should not perform measurements
during the settle time, and we need to take into account both time intervals
for proper synchronization.

If the timing of the SLM is known exactly, we can partially overcome the
latency by updating the image on the SLM _before_ the current measurement has
finished. Once the SLM starts responding (that is, after the latency), the
measurement will have completed. With these ‘pipelined’ measurements, enabled
by default in OpenWFS, an experiment can be run significantly than without,
especially if the settle time is short [39].

Figure 8:  Measured SLM response (a) when the SLM grey value is changed from 0
to 255, and (b) when the grey value is changed from 255 to 0.

####  7.4.1 Automatic detection and quantification

To measure idle and settle time interval, the phase of the SLM screen is
switched between two phase values while recording the reflected beam from the
SLM using a photodiode. Figure 8 shows an example of the detected signal
measured for Hamamatsu X13138-07 SLM. The image on this SLM is updated with a
refresh rate of 60 frames per second. When we change the grey value from 0
(black screen) to 255 (white screen) the SLM responds as shown in Fig. 8(a)
(red curve). The blue pulse represents the input signal whose peak indicates
the moment that the new grey value is sent to the SLM. From Fig. 8(a) it can
also be estimated that the signal oscillates with a period of
$8.3\text{\,}\mathrm{m}\mathrm{s}$ (half a frame) representing the phase
jitter period.

After sending the new grey value (white screen) on the SLM, the image on the
SLM does not change for at least $33\text{\,}\mathrm{m}\mathrm{s}$, which
corresponds to 2 frame periods of latency. After this latency, the SLM takes
between 20 and $40\text{\,}\mathrm{m}\mathrm{s}$ to stabilize (the _settle
time_). It can be seen that the settle time is different for switching from
white to black or from black to white.

####  7.4.2 Mitigation

The idle time and settle time should be adjusted to match the measured times
(See Fig. 8). Comparing the responses in Fig. 8(a) and Fig. 8(b), we can
conclude that the SLM (Hamamatsu X13138-07) responds faster when the grey
value is switched from 0 to 255 than when the grey value is changed from 255
to 0. Such information about the SLM timing helps us to design the fast WFS
more effectively. For instance, if we arrange the phase stepping for WFS in
ascending order, equivalently changing the grey values from 0 to 255 with a
chosen number of steps. It causes the grey value to fall from 255 to 0 in the
first iteration of the stepping which takes longer than when the grey value
rises from 0 to 255. Therefore, arranging the phase steps in descending order
results in avoiding the grey value to fall from 255 to 0, consequently making
the WFS faster.

##  8 Imperfections in the measurements

###  8.1 Noise in the feedback signal

As explained in Section 3.2, the WFS measurements consist of modulating part
of the SLM segments and finding the required phase essentially by fitting
sinusoids. In practice, the recorded feedback signal includes the effect of
different types of noise such as read-out noise, shot noise, and laser excess
noise. In the presence of noise, the measured feedback signal has deviations
from the ideal sinusoid, causing the fidelity to decrease by a factor of [45]

$\left\langle\left|\gamma_{N}\right|^{2}\right\rangle=\frac{SNR}{1+SNR}=\frac{P_{\text{sig}}}{P_{\text{sig}}+P_{\text{noise}}},$
(40)

where SNR is the signal-to-noise ratio, and $P_{\text{sig}}$ and
$P_{\text{noise}}$ are the signal and noise components of the measurement,
respectively.

####  8.1.1 Automatic detection and quantification

In the DFT of the phase-stepping experiments (see Section 3.2) component
$F_{0}$ refers to the signal offset, and $F_{-1}$ and $F_{1}$ correspond to
the sinusoidal response. All remaining components in the Fourier transform
correspond to noise, imperfections in the calibration, or non-linear response
from the sample. Here, we assume that the response is perfectly linear, so
that the only effect comes from noise. We furthermore assume white noise,
meaning that the noise is distributed equally over all Fourier components. In
particular, this means that the total power of the noise in the 1st Fourier
component, which also includes the signal, can be estimated as

$P_{\text{noise}}=\frac{1}{P-3}\sum_{\left|k\right|>1}\sum_{a}^{N}\left|F_{k,a}\right|^{2},$
(41)

where the sum is over all $P$ frequency components of the Fourier transform,
excluding the $-1$ and $1$ terms (which include the signal) and the $0$
component (which is the signal offset). Assuming that the noise is independent
of the signal, we can estimate the total signal power as

$P_{\text{sig}}=\sum_{a}^{N}\left|F_{1,a}\right|^{2}-P_{\text{noise}},$ (42)

which can be inserted into Eq. (40) to compute the fidelity reduction caused
by noise. A further refinement can be done by not treating $F_{1,a}$ as the
signal, but to take into account the calibration error and also include all
other components $c_{k}$ in Eq. (36). We do so by constructing a ‘filtered’
signal $\tilde{F}$

$\begin{bmatrix}\tilde{F}_{k,1}\\\ \tilde{F}_{k,2}\\\ \ldots\\\
\tilde{F}_{k,N}\end{bmatrix}=FF^{+}\begin{bmatrix}F_{k,1}\\\ F_{k,2}\\\
\ldots\\\ F_{k,N}\end{bmatrix}$ (43)

using the matrix $F$ and its pseudo-inverse $F^{+}$ as defined in Eq. (37). We
now treat $\tilde{F}$ as the signal instead of $F_{k,1}$. With a perfect
calibration of the SLM, these two approaches are equivalent. With a poorly
calibrated SLM, however, this approach allows separating the effect of noise
and poor calibration to avoid ‘double counting’ of the calibration fidelity.
OpenWFS computes $\left\langle\left|\gamma_{N}\right|^{2}\right\rangle$ from
$\tilde{F}$ and stores the value as fidelity_noise.

####  8.1.2 Mitigation

One technique to reduce the effect of noise is to increase the size of the
area where the feedback signal is collected from. Since more light is
collected, the SNR of the measurement increases. However, if we take the
feedback area too large, the enhancement will be distributed over multiple
targets and decreases as $|\gamma|^{2}=1/M$, where $M$ is the number of
targets [41].

A different technique to improve the SNR is to perform a pre-optimization with
a limited number of optical modes [38, 45, 30]. This technique increases the
strength of the reference field, thereby increasing the signal strength. The
drawback of this technique is that it only works when wavefront shaping for a
single target. Many algorithms allow measuring the optimal wavefront for many
targets $b$ simulaneously, for example when a camera is used as a detector
[25].

It was shown in Ref. [25] that modulating multiple segments at a time can be
used to greatly increase the SNR of the phase stepping measurements. In these
approaches the number of controlled pixels is equal to the number of reference
pixels. It leads to an optimally balanced interference of modulated and
reference field. This 50/50 choice for the reference and controlled part
guarantees the maximal interferometric visibility and thus maximizes
$P_{\text{sig}}$ described above. An example of this technique is to use part
of the SLM as reference, and perform measurements of the remaining part of the
SLM by grouping the pixels in a Hadamard basis [25]. Since the fixed reference
basically corresponds to non-modulated light, this technique does reduce the
fidelity a bit. This reduction can be overcome using a dual-reference
technique [21], of which a version working in a Fourier basis [20] is
available as the FourierDualReference algorithm.

###  8.2 Instability

A low fidelity may also be the result of an instability in the setup, or a
time varying sample. As in any optical setup, there are possible sources of
drift, vibrations, and turbulence. The transmission matrix of the sample
itself may also vary on very short timescales, such as when the sample is
perfused tissue [13].

####  8.2.1 Automatic detection and quantification

In a transmission geometry, a straightforward way to check the stability of
the setup is to capture the transmitted light through the sample before and
after the measurements. An exact value for the fidelity cannot easily be
determined from these measurements only, since it strongly depends on the
nature of the sample. However, it was found that for strongly scattering
samples the correlation coefficient between the two images gives a good
indication of the reduction of the fidelity due to the instability of the
setup [13].

In OpenWFS, we estimate the effect of decorrelation using the Pearson
correlation coefficient,

$\left\langle\left|\gamma_{D}\right|^{2}\right\rangle=\frac{\left\langle\left(I_{1}-\left\langle
I_{1}\right\rangle_{b}\right)\left(I_{2}-\left\langle
I_{2}\right\rangle_{b}\right)\right\rangle_{b}}{\sqrt{\left\langle\left(I_{1}-\left\langle
I_{1}\right\rangle_{b}\right)^{2}\right\rangle_{b}}\sqrt{\left\langle\left(I_{2}-\left\langle
I_{2}\right\rangle_{b}\right)^{2}\right\rangle_{b}}},$ (44)

where $\left\langle\cdot\right\rangle_{b}$ denotes averaging over all pixels
of the output image. For strongly scattering samples, where the intensity has
an exponential distribution, this approach is equivalent to that in Ref. [13].
The advantage of Eq. (44) is that
$\left\langle\left|\gamma_{D}\right|^{2}\right\rangle\leq 1$ regardless of the
intensity distribution of the transmitted light, and that the effect of any
background intensity $I_{\text{offset}}$ is removed. This estimate is
available as fidelity_decorrelation in the result of the troubleshoot
function.

####  8.2.2 Mitigation

There are some simple ways to increase the stability of the setup: 1) Covering
the setup using boxes around the setup, with heat sources placed outside the
boxes, decreases the air turbulence and temperature-induced drift. 2) In the
case of using a water immersion objective consider the possibility of bubbles
formation and avoid it. 3) Since heating of the SLM can cause drift, send
randomly changing patterns to the SLM for some time to ‘pre-heat’ it. 4) In
order to define the source of instability, separating the probable source from
the rest of the setup can help, for instance, turning off all equipment that
may cause vibrations or heating, and comparing the long-term stability with
each of these devices turned on individually.

Lastly, it should be noted that one should either perform the full WFS
algorithm faster than the persistence time of the dynamic sample [5], or use a
dynamic algorithm designed to continuously follow the scattering behaviour of
the sample [41, 6].

###  8.3 Photobleaching

In fluorescence microscopy, another reason that the enhancement can seem lower
than it should be, is due to bleaching of the fluorophores. Due to many WFS
measurements done on the same volume in the sample, significant photobleaching
can occur in the measured spot during the WFS measurements. For this reason,
the signal can appear to not have improved all that much, while in reality,
the sample simply photobleached. Moreover, a very successful WFS experiment
can in some cases improve the focus so much that the sample rapidly
photobleaches as soon as the computed phase correction pattern is displayed on
the SLM.

####  8.3.1 Automatic detection and quantification

Photobleaching can easily be detected by measuring the same frame before and
after the WFS measurements. If the fluorescence signal has significantly
decreased due to photobleaching, the frame taken after the measurements will
show a lower signal. Note that both frames should be taken with same phase
pattern (typically a flat phase pattern) on the SLM, in order to make a fair
comparison. To assess the performance of the WFS algorithm, enhancement should
be quantified based only on frames frames taken after the WFS measurements.

####  8.3.2 Mitigation

Photobleaching can be reduced by simply reducing the excitation power. Of
course, this comes at the cost of a lower signal to noise ratio, meaning that
the mitigation strategies presented in Section 8.1 may need to be applied.

##  9 Other considerations

###  9.1 Multi-target feedback

Throughout this paper, we considered the task of optimizing the intensity in a
single mode $b$. In practice, however, the feedback often comes from a finite-
size region that effectively spans multiple optical modes. In this case, WFS
finds a solution that simultaneously maximizes the intensity in multiple modes
simultaneously. It can be shown [41] that in this case the enhancement is
distributed over these modes, effectively causing the total enhancement to
reduce with a factor of the number of modes $M$.

In practice, this means there is a trade-off between the size of the feedback
area and the enhancement. Obtaining feedback from a larger area (e.g. by
averaging over multiple camera pixels, or using a larger fluorescent guide
star) may result in a better SNR, but choosing this area too large causes the
enhancement (and also the resolution of the focus) to decrease.

###  9.2 Non-linear feedback

Feedback in a WFS experiment can come in different forms. An interesting form
of feedback is that of multi-photon excitation fluorescence (multi-PEF, or
more specifically 2-PEF, 3-PEF, etc.). In multi-PEF, a fluorophore is excited
by multiple photons before emitting a fluorescence photon. A 2-photon
excitation spectrum of a fluorophore is typically similar to the single-photon
excitation spectrum but with the wavelength doubled [43]. In multi-PEF, the
fluorescence signal scales non-linearly with the excitation intensity. If the
fluorescence is confined to a small volume (a ‘guide star’), a narrower spot
size will produce a significantly higher non-linear signal. This can be used
as feedback signal for WFS [29].

Despite the non-linear characteristic of the feedback, Phase stepping can
still be used to retrieve the relative phase with Eq. (6), but ideally $I_{p}$
follows the shape of a $\cos^{2n}(\phi_{p}/2)$ function, with $n$ denoting the
order of non-linearity, as explained in Ref. [28]. This non-sinusoidality of
$I_{p}$ causes predictable non-zero $F_{k}$ for ${\left|k\right|\geq 2}$, even
in the absence of noise and for a perfect SLM response.

Lastly, if the fluorescence signal is not coming from a single point, the
exact relationship between the order of the non-linearity and the achievable
enhancement is complicated and depends on properties of the excitation source
and the structure of the fluorescent material in the sample. [28, 4, 42] With
2-PEF in a fluorescent dye pool, for example, forming a single focus does not
optimize the total fluorescence. Instead, WFS finds a complicated pattern of
hot-spots throughout the pool [14]. If the labelling density is sparse enough,
however, iterative WFS converges to a tight focus [14, 23]. Lastly, for non-
linear feedback, some amount of non-sinusoidality of $I_{p}$ is to be
expected. See also Section 9.2.

###  9.3 Amplitude-only light modulators

Digital micro-mirror devices (DMDs) are light modulators that modulate only
the amplitude of the incident light [5, 1]. DMDs have a much faster switching
rate than liquid crystal SLMs, which is an advantage. However, for strongly
scattering samples the amplitude-only modulation reduces the fidelity to a
maximum of $|\gamma|^{2}=1/(2\pi)\approx 0.16$, [1] which is significantly
lower than the maximum theoretical fidelity of $\pi/4\approx 0.79$ for phase-
only SLMs. For strongly scattering samples, the fidelity can be increased by
roughly a factor of $2$ by simply blocking the 0-th order reflection of the
DMD [3, 37]. Alternatively, the DMD can be used in a configuration that
enables phase modulation, or even phase-and-amplitude modulation, at the cost
of reducing it’s effective resolution [5, 10]

##  10 OpenWFS for automatic troubleshooting

⬇

import numpy as np

import astropy.units as u

from openwfs.processors import SingleRoi

from openwfs.simulation import SLM, Microscope, Shutter

from openwfs.algorithms import StepwiseSequential, troubleshoot

from openwfs.utilities import set_pixel_size

# === Define virtual devices for a WFS simulation ===

# Note: these virtual-device classes may be replaced

# with real-device classes to run a real WFS experiment

# Define aberration as a pattern of random phases at the pupil plane

aberrations = np.random.uniform(size=(40, 40)) * 2 * np.pi

# Define specimen as an image with several bright pixels

specimen = np.zeros((120, 120))

specimen[60, 60] = 1e4

specimen = set_pixel_size(specimen, pixel_size=200 * u.nm)

# Simulate an SLM with incorrect phase response

# Also simulate a shutter that can turn off the light

# The SLM is conjugated to the back pupil plane

slm = SLM(shape=(100, 100),

phase_response=(np.arange(256) / 128 * np.pi) * 1.2)

shutter = Shutter(slm.field)

# Simulate a WFS microscope looking at the specimen

sim = Microscope(source=specimen, incident_field=shutter,
aberrations=aberrations, wavelength=800 * u.nm)

# Simulate a camera device with gaussian noise and shot noise

cam = sim.get_camera(analog_max=1e4, shot_noise=True, gaussian_noise_std=4.0)

# Define feedback as circular region of interest in the center of the frame

roi_detector = SingleRoi(cam, radius=0.1)

# === Run wavefront shaping experiment ===

# Use the stepwise sequential (SSA) WFS algorithm

alg = StepwiseSequential(feedback=roi_detector, slm=slm, n_x=10, n_y=10,
phase_steps=16)

# Define a region of interest to determine background intensity

roi_background = SingleRoi(cam, radius=10)

# Run WFS troubleshooter and output a report to the console

trouble = troubleshoot(algorithm=alg, background_feedback=roi_background,

frame_source=cam, shutter=shutter)

trouble.report()

Figure 9: Example code for running a WFS simulation in OpenWFS.
(examples/troubleshooter_demo.py)

⬇

=== Feedback metrics ===

number of modes (N): 100.000

fidelity_amplitude: 0.726

fidelity_noise: 0.719

fidelity_non_modulated: 0.931

fidelity_phase_calibration: 0.884

expected enhancement: 30.154

measured enhancement: 15.920

=== Frame metrics ===

fidelity_decorrelation: 0.702

average signal offset: 1.63

photobleaching ratio: 1.001

Figure 10: Example output from the troubleshoot function.

Figure 11: Simulated wavefront shaping with OpenWFS. a. A simulated speckle
pattern. b. A focus acquired after simulated wavefront shaping under imperfect
circumstances.

In this section we show how to run a wavefront shaping experiment with our
OpenWFS framework [34]. In the Python script shown in Fig. 9, we define
virtual hardware devices to run a WFS simulation. In this example, we simulate
an SLM with an incorrect phase response by defining the SLM object with a
custom phase_response function, and a virtual camera with significant noise.
Before running the WFS algorithm, the virtual camera sensor shows a speckle
pattern (Fig. 11a). After running the WFS algorithm, the virtual camera sensor
shows a focus (Fig. 11b). Due to several detrimental effects, the focus is
significantly lower than the theoretical maximum. We show that our
troubleshooter can detect and approximately quantify the degraded quality of
the result and identify the detrimental effects.

The troubleshooter extracts much of its information from the WFS measurements
directly. These results are saved as properties of the WFSResult object that
is returned by the algorithm. We can currently simulate the detrimental
effects of unequal illumination profile, noise, incorrect phase calibration
and non-modulated light. The code is generic such that it can be used both for
simulations and for experiments, by replacing the virtual hardware classes
with real hardware classes.

We show an example output of the troubleshoot function in Fig. 10. An
enhancement of almost 16 has been reached in this experiment. Although this is
already a usable enhancement, this is not quite as high as the theoretical
value of 100, equal to the number of measured modes. The metrics give insight
into what caused the reduced enhancement. The fidelity_amplitude is
significantly lower than 1, indicating that the SLM segments were not equally
illuminated. The fidelity_noise is significantly lower than 1, indicating a
significant effect of noise. The fidelity_phase_calibration is significantly
lower than 1, indicating a problem with the SLM phase calibration. Lastly, all
these fidelity estimations are combined to compute an order of magnitude
estimate of the expected enhancement. These numbers indicate which are the
likely causes of the lower-than-ideal enhancement, and thus help in mitigating
their detrimental effects.

##  11 Conclusion

The field of WFS is full of exciting new advances, from interesting new use
cases to state-of-the-art algorithms. However, the most scientific
publications only give limited attention to the the technical details of
implementing these methods in practice.

In this paper, we described many of the problems that can degrade the
performance of a WFS experiment. Anyone starting in this field is likely to
encounter one or more of these problems, giving WFS a steep ‘learning curve’
that may frustrate researchers and needlessly slow down progress in this
exciting field.

We gave a comprehensive overview of the most common issues that may plague WFS
experiments, covering issues caused by imperfections in the feedback
measurements, imperfections of the SLM, or simply properties of the modulation
scheme.

We found that many of these issues can be detected and quantified
automatically, often using the data that is already available from the WFS
experiments themselves. These automatic troubleshooters are integrated in our
open source OpenWFS Python package [34], which also provides all the essential
boilerplate code to start working with spatial light modulators, WFS
algorithms, and for simulating WFS under realistic experimental conditions.
This way, we hope to lower the threshold for starting in this field, and to
reduce much of the frustration that may be accompanied with troubleshooting
experimental issues.

## Acknowledgments

This work was supported by Nederlandse Organisatie voor Wetenschappelijk
Onderzoek (Vidi grant 14879).

## Data and code availability

openwfs is available on https://github.com/IvoVellekoop/openwfs and on the
standard PyPI Python package repository. For this manuscript version 0.1.0rc1
was used.

## References

  * [1] Duygu Akbulut, Thomas J Huisman, Elbert G van Putten, Willem L Vos, and Allard P Mosk.  Focusing light through random photonic media by binary amplitude modulation.  Optics express, 19(5):4017–4029, 2011.
  * [2] Zizheng Cao, Xuebing Zhang, Gerwin Osnabrugge, Juhao Li, Ivo M Vellekoop, and Antonius MJ Koonen.  Reconfigurable beam system for non-line-of-sight free-space optical communication.  Light: Science & Applications, 8(1):69, 2019.
  * [3] Sri Nivas Chandrasekaran, Hans Ligtenberg, Wiendelt Steenbergen, and Ivo M Vellekoop.  Using digital micromirror devices for focusing light through turbid media.  In Emerging Digital Micromirror Device Based Systems and Applications VI, volume 8979, pages 18–27. SPIE, 2014.
  * [4] Li-Chung Cheng, Nicholas G. Horton, Ke Wang, Shean-Jen Chen, and Chris Xu.  Measurements of multiphoton action cross sections for multiphoton microscopy.  Biomedical Optics Express, 5(10):3427–3433, October 2014.
  * [5] Donald B Conkey, Antonio M Caravaca-Aguirre, and Rafael Piestun.  High-speed scattering medium characterization with application to focusing light through turbid media.  Optics express, 20(2):1733–1740, 2012.
  * [6] Alexandra d’Arco, Fei Xia, Antoine Boniface, Jonathan Dong, and Sylvain Gigan.  Physics-based neural network for non-invasive control of coherent light in scattering media.  Optics Express, 30(17):30845–30856, August 2022.
  * [7] Jorge García-Márquez, Victor López, Arturo González-Vega, and Enrique Noé.  Flicker minimization in an LCoS spatial light modulator.  Optics Express, 20(8):8431–8441, April 2012.
  * [8] Sylvain Gigan, Ori Katz, Hilton B De Aguiar, Esben Ravn Andresen, Alexandre Aubry, Jacopo Bertolotti, Emmanuel Bossy, Dorian Bouchet, Joshua Brake, Sophie Brasselet, et al.  Roadmap on wavefront shaping and deep imaging in complex media.  Journal of Physics: Photonics, 4(4):042501, 2022.
  * [9] Joseph W Goodman.  Statistical optics.  John Wiley & Sons, 2015.
  * [10] Sebastianus A. Goorden, Jacopo Bertolotti, and Allard P. Mosk.  Superpixel-based spatial amplitude and phase modulation using a digital micromirror device.  Opt. Express, 22(15):17999–18009, Jul 2014.
  * [11] Sebastianus A Goorden, Marcel Horstmann, Allard P Mosk, Boris Škorić, and Pepijn WH Pinkse.  Quantum-secure authentication of a physical unclonable key.  Optica, 1(6):421–424, 2014.
  * [12] Roarke Horstmeyer, Haowen Ruan, and Changhuei Yang.  Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue.  Nature Photonics, 9(9):563–571, September 2015.
  * [13] Mooseok Jang, Haowen Ruan, Ivo M Vellekoop, Benjamin Judkewitz, Euiheon Chung, and Changhuei Yang.  Relation between speckle decorrelation and optical phase conjugation (opc)-based turbidity suppression through dynamic scattering media: a study on in vivo mouse skin.  Biomedical optics express, 6(1):72–85, 2015.
  * [14] Ori Katz, Eran Small, Yefeng Guan, and Yaron Silberberg.  Noninvasive nonlinear focusing and imaging through strongly scattering turbid layers.  Optica, 1(3):170–174, Sep 2014.
  * [15] Joel Kubby, Sylvain Gigan, and Meng Cui.  Wavefront shaping for biomedical imaging.  Cambridge University, 2019.
  * [16] Joel A Kubby.  Adaptive optics for biological imaging.  CRC Press, 2013.
  * [17] Guanming Lai and Toyohiko Yatagai.  Generalized phase-shifting interferometry.  JOSA A, 8(5):822–827, May 1991.
  * [18] Rujia Li and Liangcai Cao.  Progress in Phase Calibration for Liquid Crystal Spatial Light Modulators.  Applied Sciences, 9(10):2012, January 2019.
  * [19] José Luis Martínez, Ignacio Moreno, María del Mar Sánchez-López, Asticio Vargas, and Pascuala García-Martínez.  Analysis of multiple internal reflections in a parallel aligned liquid crystal on silicon SLM.  Optics express, 22(21):25866–25879, 2014.
  * [20] Bahareh Mastiani, Gerwin Osnabrugge, and Ivo M Vellekoop.  Wavefront shaping for forward scattering.  Optics Express, 30(21):37436–37445, 2022.
  * [21] Bahareh Mastiani and Ivo M Vellekoop.  Noise-tolerant wavefront shaping in a hadamard basis.  Optics express, 29(11):17534–17541, 2021.
  * [22] AP Mosk.  Exploiting disorder for perfect focusing.  Nature photonics, 4(5):320–322, 2010.
  * [23] Gerwin Osnabrugge, Lyubov V Amitonova, and Ivo M Vellekoop.  Blind focusing through strongly scattering media using wavefront shaping with nonlinear feedback.  Optics express, 27(8):11673–11688, 2019.
  * [24] Martin Plöschner, Tomáš Tyc, and Tomáš Čižmár.  Seeing through chaos in multimode fibres.  Nature Photonics, 9:529–535, 2015.
  * [25] SM Popoff, G Lerosey, R Carminati, M Fink, AC Boccara, and S Gigan.  Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media.  Physical review letters, 104(10):100601, 2010.
  * [26] Stefan Rotter and Sylvain Gigan.  Light fields in complex media: Mesoscopic scattering meets wave control.  Reviews of Modern Physics, 89(1):015005, 2017.
  * [27] Marijn E. Siemons, Naomi A. K. Hanemaaijer, Maarten H. P. Kole, and Lukas C. Kapitein.  Robust adaptive optics for localization microscopy deep in complex tissue.  Nature Communications, 12(1):3407, June 2021.
  * [28] David Sinefeld, Hari P. Paudel, Dimitre G. Ouzounov, Thomas G. Bifano, and Chris Xu.  Adaptive optics in multiphoton microscopy: Comparison of two, three and four photon fluorescence.  Optics Express, 23(24):31472–31483, November 2015.
  * [29] Jianyong Tang, Ronald N. Germain, and Meng Cui.  Superpenetration optical microscopy by iterative multiphoton adaptive compensation technique.  Proceedings of the National Academy of Sciences, 109(22):8434–8439, May 2012.
  * [30] Xiaodong Tao, Tuwin Lam, Bingzhao Zhu, Qinggele Li, Marc R Reinig, and Joel Kubby.  Three-dimensional focusing through scattering media using conjugate adaptive optics with remote focusing (CAORF).  Optics express, 25(9):10368–10383, 2017.
  * [31] Abhilash Thendiyammal, Gerwin Osnabrugge, Tom Knop, and Ivo M. Vellekoop.  Model-based wavefront shaping microscopy.  Optics Letters, 45(18):5101–5104, September 2020.
  * [32] Eibert Gerjan van Putten.  Disorder-Enhanced imaging with spatially controlled light.  PhD thesis, University of Twente, 2011.
  * [33] Ivo M Vellekoop.  Feedback-based wavefront shaping.  Optics express, 23(9):12189–12206, 2015.
  * [34] Ivo M Vellekoop.  OpenWFS on Github - A framework for conducting and simulating wavefront shaping experiments in Python, July 2023.  https://github.com/IvoVellekoop/openwfs. For this manuscript version 0.1.0rc1 was used.
  * [35] Ivo M Vellekoop.  OpenWFS on PyPI, July 2023.  https://pypi.org/project/openwfs. For this manuscript version 0.1.0rc1 was used.
  * [36] Ivo M Vellekoop.  OpenWFS on readthedocs, July 2023.  https://openwfs.readthedocs.io/en/latest.
  * [37] Ivo M Vellekoop, Meng Cui, and Changhuei Yang.  Digital optical phase conjugation of fluorescence in turbid tissue.  Applied physics letters, 101(8), 2012.
  * [38] Ivo M Vellekoop and Allard P Mosk.  Focusing coherent light through opaque strongly scattering media.  Optics letters, 32(16):2309–2311, 2007.
  * [39] Ivo Micha Vellekoop.  Controlling the propagation of light in disordered scattering media.  PhD thesis, University of Twente, 2008.
  * [40] Ivo Micha Vellekoop and Allard P Mosk.  Universal optimal transmission of light through disordered materials.  Physical review letters, 101(12):120601, 2008.
  * [41] Ivo Micha Vellekoop and AP Mosk.  Phase control algorithms for focusing light through turbid media.  Optics communications, 281(11):3071–3080, 2008.
  * [42] C Xu, W Zipfel, J B Shear, R M Williams, and W W Webb.  Multiphoton fluorescence excitation: New spectral windows for biological nonlinear microscopy.  Proceedings of the National Academy of Sciences, 93(20):10763–10768, October 1996.
  * [43] Chris Xu, Rebecca Williams, Warren Zipfel, and Watt Webb.  Multiphoton cross-sections of molecular fluorophores.  Bioimaging, 4:198–207, September 1996.
  * [44] Kazushi Yamaguchi, Kohei Otomo, Yuichi Kozawa, Motosuke Tsutsumi, Tomoko Inose, Kenji Hirai, Shunichi Sato, Tomomi Nemoto, and Hiroshi Uji-i.  Adaptive Optical Two-Photon Microscopy for Surface-Profiled Living Biological Specimens.  ACS Omega, 6(1):438–447, January 2021.
  * [45] Hasan Yılmaz, Willem L Vos, and Allard P Mosk.  Optimal control of light propagation through multiple-scattering media in the presence of noise.  Biomedical optics express, 4(9):1759–1768, 2013.
  * [46] Kevin C Zhou, Ruobing Qian, Simone Degan, Sina Farsiu, and Joseph A Izatt.  Optical coherence refraction tomography.  Nature Photonics, 13:794–802, 2019.