# Synthetic aperture holographic third harmonic generation microscopy

Yusef Farah  Morgridge Institute for Research, Madison, WI 53715 USA  Gabe
Murray  Physics Department, Colorado State University, Fort Collins, CO 80523
USA  Jeff Field  Maxine Xiu  Electrical and Computer Engineering Department,
Colorado State University, Fort Collins, CO 80523 USA  Lang Wang  Morgridge
Institute for Research, Madison, WI 53715 USA  Olivier Pinaud  Mathematics
Department, Colorado State University, Fort Collins, CO 80523 USA  Randy
Bartels  Morgridge Institute for Research, Madison, WI 53715 USA  Biomedical
Engineering Department, University of Wisconsin - Madison, Madison, WI 53715
USA<EMAIL_ADDRESS>
###### Abstract

Third harmonic generation (THG) provides a valuable, label-free approach to
imaging biological systems. To date, THG microscopy has been performed using
point scanning methods that rely on intensity measurements lacking phase
information of the complex field. We report the first demonstration of THG
holographic microscopy and the reconstruction of the complex THG signal field
with spatial synthetic aperture imaging. Phase distortions arising from
measurement-to-measurement fluctuations and imaging components cause optical
aberrations in the reconstructed THG field. We have developed an aberration-
correction algorithm that estimates and corrects for these phase distortions
to reconstruct the spatial synthetic aperture THG field without optical
aberrations.

††preprint: APS/123-QED

##  I Introduction

Nonlinear microscopy has found widespread use for imaging inside of complex
scattering environments, most notably in biological tissues. Deeper imaging in
these tissues is made possible by using longer wavelengths, and corresponding
increases in the scattering length of light, compared with those used in
linear fluorescent imaging. The nonlinear scattering process of third harmonic
generation (THG) [1, 2] has emerged as a powerful label-free imaging modality
complementary to second harmonic generation (SHG) [3] and multiphoton
absorption fluorescent microscopy. [4] THG microscopy can be used to study
important cellular structures, such as lipid bodies, cell membranes,
osteocytes, large protein aggregates, biogenic crystals, myelin, and
microglial activation, and to differentiate subtypes of cancer. [5, 6, 7, 8,
9, 10, 11, 12, 13] The structural properties that govern the intrinsic
contrast mechanism of THG are ideal for complex, heterogeneous biological
samples such as a porous bone matrix and retinal layers. [14, 15, 16]
Furthermore, THG imaging has been shown to accurately quantify malignant
tumors in breast tissue by capturing irregularities in lipid bodies that vary
with microglia activation associated with different cancer subtypes and even
glioma infiltration. [6, 17, 18] Particularly as suitable laser systems become
more readily available, THG is increasingly playing an important role in
biological imaging.

Like other nonlinear microscopy methods, THG requires a sufficiently high
laser light intensity to activate nonlinear signal generation due to the weak
light–matter interaction strength. [19, 20, 21, 22] To achieve sufficient
signal intensity, short laser pulses are a prerequisite for nonlinear
microscopy. Specifically, THG imaging is nearly always implemented in a point-
scanning configuration, where a laser beam is tightly focused to a
diffraction-limited point inside of a specimen. Images are built from
assigning a fraction of the emitted power from each known focal point
location. However, point scanning comes with several drawbacks. Traditional
point-scanning multiphoton microscopes have a limited field of view (FOV). In
turn, beam scanning requires an extensive amount of time to form an image with
a large FOV. And as THG microscopy relies on measurements of the THG scattered
power, THG microscopes are not able to capture the phase of the THG field.
This results in a loss of information of the complex-valued (i.e., amplitude
and phase) THG fields. To date, no THG microscopes have provided access to the
amplitude and phase information of THG light scattered by a sample. The
development and expansion of THG imaging methods is necessary to push the
capabilities of THG microscopy for biological imaging.

In this article, we introduce holographic imaging to directly record complex-
valued THG fields, providing the first demonstration of amplitude and phase-
sensitive THG holographic microscopy. This microscope is configured as a
widefield, stage-scanning THG holography that enables large FOV imaging thanks
to the large illumination area and stage scanning that is not limited by
optical components. Use of the widefield detection scheme with a larger
illumination beam than point focusing reduces the time needed for recording
images, while use of a lower local illumination intensity avoids photo
damaging the sample. [23, 24, 25, 26] Further, we implement off-axis digital
holography that provides heterodyne amplification to the weak THG signal and
preserves the phase information of the complex-valued THG field for each
illumination position.

This article also redresses a major obstacle for conventional THG imaging:
that there is no recourse for correcting THG field distortions arising from
aberrations that lead to poor image quality. These aberrations originate from
scattered light and wavefront phase distortions introduced by imperfections in
the optics and inhomogeneities in the specimen. We have developed a
computational adaptive optics [27, 28, 29, 30, 31] approach to estimate and
correct aberrations of the THG image field while preserving the signal’s phase
information, allowing us to coherently combine reconstructed holograms to
produce spatial synthetic aperture aberration-free THG images. We test our
system on a diverse representation of samples, demonstrating its versatility.
Aside from using the THG phase for aberration estimation and correction, we
expect that the acquisition of phase information will prove valuable to THG
three-dimensional tomography [25] and to study the physical origin of complex-
valued nonlinear susceptibility in both biological systems and materials.

##  II Third harmonic generation holography

Our synthetic aperture THG field imaging system is built as an off-axis
holographic microscope that is similar in design to SHG holography and
tomography systems. [23, 24, 25, 31] The experimental arrangement shown in
Fig. 1 is an interferometric technique that allows us to quickly reconstruct
complex images while simultaneously applying a heterodyne amplification [24]
to boost weak THG signals.

Because THG holography requires measuring interferograms with a low noise
camera, it allows for a widefield imaging scenario. We record the intensity
pattern formed through interference between a reference field, $E_{\rm
ref}(\mathbf{x}_{m})$ and a signal field $E_{\rm sig}(\mathbf{x}_{m})$. These
fields are provided by splitting the output of an ultrafast laser pulse from a
Yb:fiber laser-amplifier system into a fundamental beam that illuminates the
sample and a beam used to generate the reference field. The fundamental field
is gently focused (to a beam diameter of $\sim 7\mu$m, or roughly $15\times$ a
typical point focus) to illuminate the sample. The signal field is generated
through THG scattering produced by the quasi-widefield illumination of the
sample by the fundamental beam. The THG scattered field is imaged with a 4-f
microscope consisting of a UV objective lens and a 200-mm tube lens. The
coherent THG reference beam is generated by tightly focusing to the back
surface of a 0.5-mm thick TiO2 window, where strong THG light is generated.
[1, 2, 32] The THG reference light is collimated and then combined with the
THG signal field thin film polarizer positioned at 45 $\deg$. The combined
beam passes through a polarizer and the interference pattern is detected on an
EMCCD camera. Both the signal and reference fields are at an optical
frequency, $\omega_{3}=3\,\omega_{1}$ that is three times the fundamental
frequency $\omega_{1}$. These two fields interfere at the camera surface to
form a holographic intensity pattern $I_{h}(\mathbf{x}_{m})=|E_{\rm
ref}(\mathbf{x}_{m})|^{2}+|E_{\rm sig}(\mathbf{x}_{m})|^{2}+E_{\rm
ref}^{*}(\mathbf{x}_{m})\,E_{\rm sig}(\mathbf{x}_{m})+E_{\rm
ref}(\mathbf{x}_{m})\,E_{\rm sig}^{*}(\mathbf{x}_{m})$. We numerically isolate
the interferometric pseudo-intensity term $\tilde{I}_{h}=E_{\rm
ref}^{*}(\mathbf{x}_{m})\,E_{\rm sig}(\mathbf{x}_{m})$ to estimate the complex
THG image field using standard protocols. [23, 24] For each sample, a
reference (i.e., ground truth) image is taken by illuminating the sample with
a UV LED lamp to record an incoherent brightfield widefield image of the
specimen to compare to the synthetic aperture THG images. A complete
description of the experimental setup is detailed in the Supplemental
Materials.

Each measured hologram at a sample position $\mathbf{x}^{(j)}_{s}$ produces a
two-dimensional (2D) complex-valued THG field map, $E^{(j)}_{\rm
sig}(\mathbf{x}_{m})$. An example set obtained by scanning $\mathbf{x}_{s}$ is
illustrated in Fig. 1. The spatial coordinates of the camera, defined by the
pixel pitch, are mapped to measurement coordinates, $\mathbf{x}_{m}$, by
dividing by the magnification, $M$, of the THG 4-f imaging system. This
magnification was measured to be $M=47.8$ by imaging a known image test
target.

The synthetic spatial aperture is produced by scanning the specimen and
recording the THG hologram for each sample position $\mathbf{x}_{s}$. Datasets
of THG holograms are acquired by scanning the sample to both extend the
spatial region that is imaged and to ensure sufficient overlap between
adjacent hologram images, which in turn ensures convergence of our
computational adaptive optics algorithm. Sufficient overlap of the hologram is
obtained by translating the stage by roughly half of the illumination spot
size. A stack of THG 2D field image data is formed by recording the THG field
from specimen positions scanned over a large spatial aperture of the sample.
This stack is then stored for data processing.

Each 2D THG image, $E^{(j)}_{\rm sig}(\mathbf{x}_{m})$, in the data stack is
flattened into a linear vector that is inserted into a column of a matrix,
$D(\mathbf{x}_{m},\mathbf{x}_{s})$ as illustrated in Fig. 1. The column for
each specimen scan position, $\mathbf{x}^{(j)}_{s}$ is set to match the
lexicographical order of the $\mathbf{x}_{m}$ coordinates when the 2D field
data are flattened. Because we are moving the specimen relative to the fixed
THG holography image system, the matrix formed from the THG hologram data is a
direct recording of the THG nonlinear distortion operator (NDO). [31]

Another consequence of the fact that we move the specimen is that the
illumination fundamental field, and thus the effective third-order
illumination field for the THG scattering process, is the same for each
measurement. Moreover, the imaging model for the collection of the THG
scattered fields that is imaged onto the holography camera is also identical
for each measurement. Thus, we have an intrinsically isoplanatic imaging
model, provided that the imaging system for the THG scattered field is
isoplanatic, which is exploited later for aberration correction.

However, a complication arises in the recording of the THG holograms. Due to
random phase fluctuations arising from air density fluctuations and mechanical
vibrations in the relative path between the signal and reference THG fields,
there is a random phase difference between the signal and reference fields for
each THG hologram measurement. Given the weak THG field strength, we require
an integration time for each measured THG hologram that is long compared to
the timescale for the relative phase fluctuations. If we directly integrate
the hologram, the fringes will shift during the camera acquisition time,
drastically reducing fringe visibility and thus degrading the THG field
signal.

This phase fluctuation problem is mitigated by taking many THG holograms with
short integration times for each position of the sample. Each THG field is
extracted and cropped. The set of extracted fields are each flattened then
stacked sequentially into columns of a matrix. Because each THG field is
identical, except for a random mean phase difference between the signal and
reference field and for random noise acquired on each measurement, an
averaged, high signal-to-noise ratio THG field can be estimated from the
dominant singular vector of the measurement matrix [31]; the averaged field
extracted from the dominant singular vector constitutes the measured SHG
field, $E^{(j)}_{\rm sig}(\mathbf{x}_{m})$, for the jth sample position
$\mathbf{x}^{(j)}_{s}$ with a mean phase $\phi_{r}$.

##  III THG holographic aberration correction

Phase jumps between hologram measurements and pupil plane phase distortions
contribute to phase aberrations that prevent the extraction of a clean THG
field image. Distortions to the estimated images are illustrated by the “No
correction” smiley faces in Fig. 1. Without any corrections, the image is a
poor-quality reconstruction of the THG field. Therefore, we have developed an
algorithm to estimate and correct the aberrations present in the synthetic
aperture THG holographic imaging process.

Figure 1: Diagram illustrating the collection method and aberration correction
algorithm. The fundamental beam is shown in red and the generated THG is in
purple. The collected THG signal and reference are overlapped on the detector
to generate the hologram. Reconstructed holograms are flattened, stacked into
the NDO, $D(\mathbf{x}_{m},\mathbf{x}_{s})$, then input to our aberration
correction algorithm. In the first step of the algorithm, the THG field images
are shifted, mapping $\mathbf{x}_{m}\rightarrow\mathbf{x}_{o}+\mathbf{x}_{m}$,
which transforms the NDO into the transmission matrix,
$T(\mathbf{x}_{o},\mathbf{x}_{s})$. THG field images synthesized with
$T(\mathbf{x}_{o},\mathbf{x}_{s})$ using $\hat{E}_{3}(\mathbf{x}_{m})=\int
T(\mathbf{x}_{o},\mathbf{x}_{s})d^{2}\mathbf{x}_{s}$ from the raw data are
highly distorted (no correction smiley face). The dominant relative phase
aberration, $\phi_{r}(\mathbf{x}_{s})$ is corrected with modified DASH
optimization algorithm to produce a corrected transmission matrix that
produces an improved THG field image (DASH corrected smiley face). Then, we
further correct pupil plane aberrations with the spatial frequency-spatial
NDO, $D(\mathbf{u}_{o},\mathbf{x}_{s})$, by estimating and correcting the
pupil plane distortions $\phi_{o}(\mathbf{u}_{o})$. The last step rescans the
corrected spatial-spatial NDO, $\tilde{D}(\mathbf{x}_{m},\mathbf{x}_{s})$ to
the coordinates of the corrected transmission matrix,
$\tilde{T}(\mathbf{x}_{o},\mathbf{x}_{s})$. Spatial filtering and truncation
of the singular value decomposition (SVD) of
$\tilde{T}(\mathbf{x}_{o},\mathbf{x}_{s})$ are used to further denoise the
final image. The cartoon smiley faces represent the evolution of the
synthesized THG image reconstruction as the algorithm proceeds, starting with
no correction, advancing to partial correction, and culminating with full
correction. ${{\mathcal{F}\\{\cdot\\}}}$ represents a fast Fourier transform
operation that is being performed in the transformation to the
$D(\mathbf{u}_{o},\mathbf{x}_{s})$. The summation steps that form the THG
images are also followed by a reshaping step (not shown) that reverses the
flattening of the image.

The first correction in our algorithm corrects for the dominant phase
distortion, $\phi_{r}(\mathbf{x}_{s})$ that arises from the hologram
measurement process. This phase is uniformly described over 2$\pi$, leading to
speckle on an estimated image. The first step to correcting
$\phi_{r}(\mathbf{x}_{s})$ maps the spatial-spatial NDO to the THG
transmission matrix $T(\mathbf{x}_{o},\mathbf{x}_{s})$. We built an
optimization algorithm that maximizes the cost function, $J=\int I^{c}_{\rm
THG}(\mathbf{x}_{o})\,d^{2}\mathbf{x}_{o}$, that described the power of the
image, $I^{c}_{\rm
THG}(\mathbf{x}_{o})=\left|\int\tilde{T}(\mathbf{x}_{o},\mathbf{x}_{s})\,d^{2}\mathbf{x}_{s}\right|^{2}$,
estimated from a phase-corrected transmission matrix,
$\tilde{T}(\mathbf{x}_{o},\mathbf{x}_{s})=T(\mathbf{x}_{o},\mathbf{x}_{s})\,\exp[-i\phi_{r}(\mathbf{x}_{s})]$.
Because aberration phase from the hologram phase difference induces a phase
shift between otherwise identical complex vectors for each $\mathbf{x}_{o}$
position of the THG field, a correction of $\phi_{r}(\mathbf{x}_{s})$
maximizes the estimated object intensity at each spatial point by exploiting
the fact that the spatial map is highly correlated along $\mathbf{x}_{o}$.
This is possible because the data shares the spatial structure of the sample
nonlinear susceptibility, $\chi^{(3)}({{\bf x}})$. We implemented an efficient
optimization algorithm that was modified from the dynamic adaptive scattering
compensation holography (DASH) that was originally developed for wavefront
shaping in scattering media. Details of the optimization algorithm are
provided in the supplemental information. With $\phi_{r}(\mathbf{x}_{s})$
correction applied to the THG field image, a clearer image can be
reconstructed, depicted in the “DASH corrected” smiley faces in Fig. 1.

In the next step of the algorithm, we take advantage of the isoplanatic nature
of the THG holographic imaging system to correct the image aberrations
introduced by the pupil plane phase distortions, $\phi_{o}(\mathbf{u}_{o})$,
by using the concept of the distortion operator [29] that is adapted for
nonlinear imaging. [31] To use the NDO for aberration correction, we map the
“DASH corrected” spatial-spatial transmission matrix,
$T(\mathbf{x}_{o},\mathbf{x}_{s})$, into the spatial frequency-spatial
distortion operator $D(\mathbf{u}_{o},\mathbf{x}_{s})$. To do so, we first
shift the spatial fields back from $\mathbf{x}_{o}$ to $\mathbf{x}_{s}$ to
produce $D(\mathbf{x}_{m},\mathbf{x}_{s})$. Then each 2D output field for each
$\mathbf{x}_{s}$ is Fourier transformed to the spatial frequency space, so
that the NDO provides a relationship between $\mathbf{x}_{s}$ and the output
spatial frequency coordinates $\mathbf{u}_{o}$. Because the THG holographic
imaging system is well described by an isoplanatic model, the correction phase
can be estimated from the dominate singular vector of
$D(\mathbf{u}_{o},\mathbf{x}_{s})$. Following a similar procedure for the
linear DO, [29] we take the phase and amplitude of the dominant singular
vector across the output spatial frequency coordinate and apply the inverse,
with a regularization parameter, to each $\mathbf{u}_{o}$ across the
$\mathbf{x}_{s}$ spatial positions. Then, an inverse Fourier transform of
$\mathbf{u}_{o}$ brings back the spatial-spatial NDO,
$\tilde{D}(\mathbf{x}_{m},\mathbf{x}_{s})$, which has been aberration
corrected for hologram measurement phase distortions and pupil plane phase
distortions.

In the final step of our algorithm, we apply two additional noise-reduction
measures to further improve the final THG field image, as depicted as the
“Fully corrected” smiley faces in Fig. 1. First, we transform the data back
into the transmission matrix, $T(\mathbf{x}_{o},\mathbf{x}_{s})$. Taking
advantage of the fact that the fundamental illumination field spans a finite
spatial support, set by the width of the fundamental illumination beam, we
implement a spatial filtering for each $\mathbf{x}_{s}$ location analogously
to smart-OCT filtering [28] to eliminate multiple scattering. Finally, we
apply a truncated SVD and retain the singular values that contribute greater
than 5% to the total image intensity variance to further reduce noise in the
image, resulting in the final corrected and complex-valued THG image, the
results of which are shown in the next section.

##  IV Results

We applied our synthetic aperture THG imaging system and aberration correction
algorithm to several samples. All reconstructed hologram images have been
propagated via angular spectrum propagation to ensure images are in focus, a
conventional method for hologram reconstruction. [23] Focal propagation
distances are estimated using a gradient sharpness criterion.

The first set of images shown in Fig. 2 are of a $6\mu$m thick mouse tail
cross section prepared on a glass slide (Triarch Incorporated). THG
reconstructed images are of the epidural region of the mouse tail cross
section, with the LED ground truth image shown in Fig. 2a. The mouse tail
sample demonstrates our THG holographic image reconstruction and algorithm
corrections for a biological sample with a continuous structure, in this case
the long and continuous strands of collagen fibers. Before any correction to
the transmission matrix, the THG intensity and phase images are poor
reconstructions of the mouse tail sample, as indicated in Fig. 2b,c. After
first correcting the hologram phase distortion $\phi_{r}(\mathbf{x}_{s})$, the
image quality improves significantly, as shown in Fig. 2d. The fully
corrected, aberration-free THG holographic images are shown in 2e,f.
Individual tissue strands in intensity image Fig. 2e are well resolved,
showing the fibers as continuous strands that wrap around the crimped region
near the middle of the FOV. The recorded phase of individual fibers of Fig. 2f
can be identified in the phase map of the aberration-free image, reporting
characteristics of sample morphology.

Figure 2: Mouse tail cross section. a) Brightfield image. b) THG intensity
image with no correction. c) Phase map of THG image with no correction. d)
Partially corrected THG intensity image with only hologram distortion phase
$\phi_{r}(\mathbf{x}_{s})$ corrected. e) Final THG intensity image after all
corrections have been applied by our aberration correction algorithm. f) Final
THG phase map with all corrections. Scale bars are 10 $\mu$m.

The next biological sample we demonstrate our THG holographic microscope on is
a developing bone sample, approximately $6\mu$m thick (Triarch Incorporated),
shown in Fig. 3. The THG reconstructed image shows several osteoblast cells,
commonly found in developing bone. Osteoblasts are discrete cells within the
bone and offer a good contrast to the continuous structure seen in the mouse
tail images. When we apply our algorithm, the final THG intensity image
reconstruction shown in Fig. 3c is of much higher quality compared to the
intensity image shown in panel (b). The osteoblast cells have a discernible
structure in the fully corrected THG intensity image, comparable to the
brightfield image. Further, the phase map of the reconstructed field is well
resolved and able to capture the signal’s phase from individual osteoblast
sites. The phase corrections calculated with our algorithm are shown in Fig.
3e, with the top panel’s phase map showing the phase corrections applied to
the transmission matrix to correct hologram phase distortions,
$\phi_{r}(\mathbf{x}_{s})$ and the bottom panel showing the pupil plane phase
corrections acquired by the NDO step in our algorithm.

Figure 3: Developing bone. a) Brightfield image. b) THG intensity image with
no correction. c) Fully corrected THG intensity image. d) Phase map of fully
corrected THG image. e) (Top) The first corrections applied in our algorithm
to estimate and correct $\phi_{r}(\mathbf{x}_{s})$. (Bottom) NDO pupil plane
correction phase map, $\phi_{o}(\mathbf{u}_{o})$. Scale bars are 10 $\mu$m.

The final THG reconstructed image presented here is that of monolayer
molybdenum disulfide (MoS2) (Fig. 4). Monolayer MoS2 is a two-dimensional
semiconductor with a thickness of $<$1 nm and has shown significant
enhancements to the third-order susceptibility through exciton resonance. In
our system, the THG wavelength (355 nm) is resonant with the C-exciton and,
therefore, provides a strong THG signal considering that the material is only
three atoms thick. The MoS2 monolayers were originally grown on a sapphire
substrate via chemical vapor deposition, and then mechanically transferred to
a glass microscope slide. During the mechanical transfer, wrinkles and slight
defects form in the monolayer and are present in the images shown here.
Similar to the biological images shown in Fig. 2 and Fig. 3, we see
significant improvement to the reconstructed THG images after applying our
aberration correction algorithm. Most interestingly, the phase map of the
fully corrected image shows spatial variation of phase, even though the sample
is atomically thin.

Figure 4: MoS2 monolayer. a) Brightfield image. b) THG intensity image with no
correction. c) Phase map of THG image with no correction. d) Fully corrected
THG intensity image. e) Phase map of fully corrected THG image. Scale bars are
10 $\mu$m.

##  V Discussion

Here, we have demonstrated the first amplitude and phase THG images taken with
the first THG holographic imaging experiment. While widefield SHG imaging and
holography have been implemented by many groups, THG microscopy has remained
confined to point-scanning microscopy except for a recent report on intensity
widefield THG with a low-repetition rate amplifier. [33] The phase of the THG
field obtained from these images provides access to a new physical
measurement. We are still developing a detailed theory of THG scattering and
holographic imaging to fully quantify the contributions to the THG phase, yet
some conclusions can be readily drawn.

To understand THG holographic image formation, we must understand the
properties of THG scattering and how the scattered field is modified by the
incident fundamental field, the spatial organization of the linear,
$\chi^{(1)}({{\bf x}})$, and nonlinear, $\chi^{(3)}({{\bf x}})$,
susceptibility, and the orientation of the harmonophores within the
interaction volume. Unlike SHG, the THG scattering process is symmetry allowed
in all media, meaning that THG light is always generated throughout the full
focal volume.

The differences between THG imaging with point scanning compared to widefield
illumination are best illustrated by comparing the THG scattering between
these two cases in a sample exhibiting a spatially uniform third-order
susceptibility ($\chi^{(3)}({{\bf x}})$). Because ($\chi^{(3)}$) is symmetry
allowed everywhere, in a uniform medium THG is generated everywhere at a rate
determined by the local strength of the electric field of the fundamental
laser beam. In the case of a tightly focused beam, the scattered THG light
accumulates a phase shift three times the fundamental field phase, and thus
the phase of the local THG scattering depends on the fundamental field. The
focused fundamental field contains the Gouy phase shift through the focus that
is transferred to the THG field. Because the Gouy phase shift is an odd
function about the focal plane of the beam, there is complete destructive
interference of THG light scattered after the focus with light generated
before the focus. Thus, in a uniform medium, no THG light is generated from a
tightly focused beam. [34] Conventional laser-scanning THG imaging produces no
THG signal from a homogeneous sample and only produces a signal when symmetry
of the spatial distribution of the $\chi^{(3)}({{\bf x}})$ is broken by either
an interface or a small inclusion along the direction of propagation of the
beam. [1] Images are formed in this configuration by recording some fraction
of the scattered THG power and such measurements lose all phase information.
In short, point-scanning THG imaging does not provide access to the phase of
the THG field.

By contrast, consider the case in which a plane wave fundamental beam
illuminates the uniform sample. This scenario will produce a plane wave at the
third harmonic frequency of the fundamental produced by THG scattering. As the
phase shift accumulated with propagation by the fundamental plane wave
advances linearly, there is a measurable THG field after the homogeneous
sample because destructive interface across the volume is avoided. This
behavior is in stark contrast to the tightly focused case. While in our
experiment we are gently focusing the fundamental beam, our sample thickness
is thinner than the confocal parameter of the beam, $k_{1}\,w_{0}^{2}\sim 0.5$
mm, where $k_{1}=n_{1}\,\omega_{1}/c$, $n_{1}$ is the refractive index (RI) at
the fundamental wavelength, and $w_{0}$ is the radius of the beam focus.
Within the confocal parameter, the beam propagation can be reliably modeled as
a plane wave.

Relative propagation effects between the fundamental and third harmonic field,
referred to as phase matching, influence the interference of THG fields
scattered at various spatial locations within the sample, thereby impacting
the overall amplitude and phase of the THG field. For these two cases
discussed above, we tacitly assume that there is no phase mismatch, i.e., that
the RI at the fundamental and third harmonic wavelengths, $n_{1}$ and $n_{3}$,
are the same. In normal materials, $\Delta n=n_{3}-n_{1}>0$ due to the
dispersion of the RI. The collinear phase mismatch parameter is defined as
$\Delta k=k_{3}-3k_{1}=\Delta n=n_{3}-n_{1}>0$. Here we have also written the
THG wavenumber as $k_{3}=n_{3}\,\omega_{3}/c$. The phase mismatch also affects
the magnitude of the THG scattered power in both geometries, but the overall
sensitivity of the THG imaging modalities to heterogeneity in
$\chi^{(3)}({{\bf x}})$ is not substantially modified.

Another source of phase variation in the THG holographic images occurs when
the THG field propagates through the specimen’s spatial RI heterogeneity,
reflected in $n_{3}({{\bf r}})$, that is imprinted onto the THG field. In
addition, spatial inhomogeneities in $n_{1}({{\bf r}})$ can further distort
the fundamental field that drives THG scattering, thereby altering the
amplitude and phase of the total THG field. These two effects arising from RI
inhomogeneities provide a complicated relationship between the THG field phase
and the sample spatial heterogeneity in the linear and nonlinear optical
susceptibility. Analysis of these effects will be the focus of future work.
However, when linear scattering can be neglected, the spatial distribution of
the third-order susceptibility can be estimated by using harmonic optical
tomography. [25]

Because the $\chi^{(3)}$ is a tensor, the polarization of the fundamental
field and the orientation of the harmonophore that generates the THG
scattering also influence the complex THG image field. In many samples, the
harmonophores that generate the nonlinear signals for THG are isotropically
oriented, which reduces the number of THG susceptibility tensor elements and
leads to a lack of angular momentum conservation for THG scattering when the
specimen is illuminated with a circularly polarized beam. This conservation of
angular momentum can be exploited for increased spatial resolution in THG
imaging [35], as can structured illumination [36]. When the harmonophores are
not isotopically oriented, such as with a birefringent, THG images can be
recorded with circularly polarized incident light crystal. [5, 37, 38]

THG phase variations can also arise from the physics of the third-order THG
optical tensor. A complex-valued $\chi^{(3)}$ will lead to phase shifts of
scattered THG light that will be reflected the phase of the imaged THG field.
Such a contribution to the THG phase is evident in the spatial variation of
the phase in THG field image of MoS2 in Fig. 4; as this sample is only three
atomic layers thick, there is no possibility of phase accumulation from
propagation. We speculate that the origin of this phase response arises from
the spatial localization of the C-excitation in these systems, [39, 40] but
confirmation of this conjecture requires further investigation. In general,
linear and multiphoton resonances will also lead to phase variations that
arise from chemical-specific absorption resonances in molecular systems. For
instance, previous spectroscopy measurements of THG scattering have revealed
two-photon absorption resonances in oxy- and deoxy-hemoglobin. [41]

Because $\chi^{(3)}$ is a weak quantity, we require an intense laser field to
drive the THG scattering process sufficiently strongly to observe THG images.
This weak interaction requires that we lightly focus the fundamental
illumination light, even in the case of widefield imaging. The result is that
our FOV is limited by practical considerations. By recording the THG field as
a hologram, we benefit from heterodyne amplification of the THG field,
allowing us to relax focusing conditions. However, to record images from a
sufficiently large FOV that we are able to observe spatial morphologies of
samples, we take a series of images as we scan the sample position. In this
way, we synthesize a FOV that is not restricted by the transverse spatial
extent of the fundamental illumination beam. In conventional THG imaging,
stitching together a larger FOV is quite simple because we simply need to add
together intensities on the boundaries. [42] In the case of holographic
imaging, because we have access to the field information, we are extremely
sensitive to phase errors when stitching together many images to synthesize a
large image FOV.

We have created an algorithm to both estimate and correct these random phase
errors that inevitably occur in any experimental system to realize the
synthesis of a large FOV for THG field imaging. The first step is to remove
the phase jumps that occur for each measured THG field. Because these fields
are recorded as a hologram, they are subject to random phase variations
between measurements. This random phase noise induces phase discontinuities
between different sample positions that ruin the correlations identified in
the NDO. [31] Without this first step, the NDO will not pick up on pupil plane
phase aberrations because phase discontinuities from different sample
positions hide the correlations in the NDO. The use of these phase corrections
is critical for the implementation of THG holography because we need to scan
the specimen to increase the FOV to observe meaningful spatial morphology;
however, without phase aberration corrections, the images are too distorted to
be useful. Once these aberrations are corrected, the synthetic aperture THG
image field is then obtained from the corrected spatial-spatial THG
transmission matrix.

##  VI Conclusion

The results here show the incredible capabilities of our aberration-free,
synthetic aperture THG holography for biological imaging. Not only have we
demonstrated THG holography for the first time, but we have also developed an
algorithm to reconstruct the THG field across a large FOV, preserve phase
information, and correct distortions and aberrations. By scanning the
specimen, we directly measure the NDO. The NDO for aberration correction is
directly accessible with our experimental setup of stationary illumination.
Measurements of the phase of the THG field are not accessible in conventional
point-scanning THG microscopy. This phase information provides new measurement
capabilities because the phase information reveals underlying resonances that
introduce an imaginary component to the nonlinear optical susceptibility,
i.e., $\mathrm{Im}\\{\chi^{(3)}({{\bf x}})\\}$. Several exciting avenues are
opened by this work. By acquiring the phase of the THG signal, for example, we
open a pathway for harmonic tomographic imaging with third harmonic light,
while the heterodyne amplification of the THG signal opens the path for
widefield THG imaging deep inside of scattering media.

##  VII Acknowledgements

We are grateful for funding support from the Chan Zuckerberg Initiative.
Olivier Pinaud is supported by NSF grant DMS-2006416. Alicia Williams from the
Morgridge Institute provided scientific editing support. Matt Stefely from the
Morgridge Institute made Fig. 1. The MoS2 sample was prepared on a glass slide
and provided to us by Prof. Justin Sambur’s group at Colorado State
University.

## References

  * Barad  _et al._ [2003] Y. Barad, H. Eisenberg, M. Horowitz, and Y. Silberberg, SPIE milestone series 175, 47 (2003).
  * Squier  _et al._ [1998] J. A. Squier, M. Müller, G. Brakenhoff, and K. R. Wilson, Optics express 3, 315 (1998).
  * James and Campagnola [2021] D. S. James and P. J. Campagnola, BME frontiers (2021).
  * Carriles  _et al._ [2009] R. Carriles, D. N. Schafer, K. E. Sheetz, J. J. Field, R. Cisek, V. Barzda, A. W. Sylvester, and J. A. Squier, Review of scientific instruments 80 (2009).
  * Oron  _et al._ [2004] D. Oron, D. Yelin, E. Tal, S. Raz, R. Fachima, and Y. Silberberg, Journal of structural biology 147, 3 (2004).
  * Débarre  _et al._ [2006] D. Débarre, W. Supatto, A.-M. Pena, A. Fabre, T. Tordjmann, L. Combettes, M.-C. Schanne-Klein, and E. Beaurepaire, Nature methods 3, 47 (2006).
  * Witte  _et al._ [2011] S. Witte, A. Negrean, J. C. Lodder, C. P. De Kock, G. Testa Silva, H. D. Mansvelder, and M. Louise Groot, Proceedings of the National Academy of Sciences 108, 5970 (2011).
  * Farrar  _et al._ [2011] M. J. Farrar, F. W. Wise, J. R. Fetcho, and C. B. Schaffer, Biophysical journal 100, 1362 (2011).
  * Lim  _et al._ [2014] H. Lim, D. Sharoukhov, I. Kassim, Y. Zhang, J. L. Salzer, and C. V. Melendez-Vasquez, Proceedings of the National Academy of Sciences 111, 18025 (2014).
  * Kuzmin  _et al._ [2016] N. Kuzmin, P. Wesseling, P. de Witt Hamer, D. Noske, G. Galgano, H. Mansvelder, J. Baayen, and M. Groot, Biomedical optics express 7, 1889 (2016).
  * Weigelin  _et al._ [2016] B. Weigelin, G.-J. Bakker, and P. Friedl, Journal of Cell Science 129, 245 (2016).
  * Tokarz  _et al._ [2017] D. Tokarz, R. Cisek, M. N. Wein, R. Turcotte, C. Haase, S.-C. A. Yeh, S. Bharadwaj, A. P. Raphael, H. Paudel, C. Alt, _et al._ , PloS one 12, e0186846 (2017).
  * Gavgiotaki  _et al._ [2020] E. Gavgiotaki, G. Filippidis, V. Tsafas, S. Bovasianos, G. Kenanakis, V. Georgoulias, M. Tzardi, S. Agelaki, and I. Athanassakis, Scientific Reports 10, 11055 (2020).
  * Genthial  _et al._ [2017] R. Genthial, E. Beaurepaire, M.-C. Schanne-Klein, F. Peyrin, D. Farlay, C. Olivier, Y. Bala, G. Boivin, J.-C. Vial, D. Débarre, _et al._ , Scientific Reports 7, 3419 (2017).
  * Genthial  _et al._ [2019] R. Genthial, M. Gerbaix, D. Farlay, L. Vico, E. Beaurepaire, D. Débarre, and A. Gourrier, PLoS One 14, e0209079 (2019).
  * Masihzadeh  _et al._ [2015] O. Masihzadeh, T. C. Lei, S. R. Domingue, M. Y. Kahook, R. A. Bartels, and D. A. Ammar, Molecular Vision 21, 538 (2015).
  * Gavgiotaki  _et al._ [2017] E. Gavgiotaki, G. Filippidis, H. Markomanolaki, G. Kenanakis, S. Agelaki, V. Georgoulias, and I. Athanassakis, Journal of biophotonics 10, 1152 (2017).
  * Zhang  _et al._ [2019] Z. Zhang, J. C. de Munck, N. Verburg, A. J. Rozemuller, W. Vreuls, P. Cakmak, L. M. van Huizen, S. Idema, E. Aronica, P. C. de Witt Hamer, _et al._ , Advanced Science 6, 1900163 (2019).
  * XU  _et al._ [2003] C. XU, W. ZIPFEL, J. B. SHEAR, R. M. WILLIAMS, and W. W. WEBB, SPIE milestone series 175, 501 (2003).
  * Erikson  _et al._ [2007] A. Erikson, J. Örtegren, T. Hompland, C. de Lange Davies, and M. Lindgren, Journal of biomedical optics 12, 044002 (2007).
  * Fuentes-Corona  _et al._ [2019] C. G. Fuentes-Corona, J. Licea-Rodriguez, R. Younger, R. Rangel-Rojo, E. O. Potma, and I. Rocha-Mendoza, Biomedical optics express 10, 6449 (2019).
  * Débarre and Beaurepaire [2007] D. Débarre and E. Beaurepaire, Biophysical journal 92, 603 (2007).
  * Masihzadeh  _et al._ [2010] O. Masihzadeh, P. Schlup, and R. A. Bartels, Optics express 18, 9840 (2010).
  * Smith  _et al._ [2013] D. R. Smith, D. G. Winters, and R. A. Bartels, Proceedings of the National Academy of Sciences 110, 18391 (2013).
  * Hu  _et al._ [2020] C. Hu, J. J. Field, V. Kelkar, B. Chiang, K. Wernsing, K. C. Toussaint, R. A. Bartels, and G. Popescu, Nature Photonics 14, 564 (2020).
  * Fantuzzi  _et al._ [2023] E. M. Fantuzzi, S. Heuke, S. Labouesse, D. Gudavicius, R. Bartels, A. Sentenac, and H. Rigneault, Nature Photonics 17, 1097 (2023).
  * Adie  _et al._ [2012] S. G. Adie, B. W. Graf, A. Ahmad, P. S. Carney, and S. A. Boppart, Proceedings of the National Academy of Sciences 109, 7175 (2012).
  * Badon  _et al._ [2016] A. Badon, D. Li, G. Lerosey, A. C. Boccara, M. Fink, and A. Aubry, Science Advances 2, e1600370 (2016).
  * Badon  _et al._ [2020] A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, and A. Aubry, Science Advances 6, eaay7170 (2020).
  * Kang  _et al._ [2017] S. Kang, P. Kang, S. Jeong, Y. Kwon, T. D. Yang, J. H. Hong, M. Kim, K.-D. Song, J. H. Park, J. H. Lee, M. J. Kim, K. H. Kim, and W. Choi, Nature Communications 8, 2157 (2017).
  * Murray  _et al._ [2023] G. Murray, J. Field, M. Xiu, Y. Farah, L. Wang, O. Pinaud, and R. Bartels, Opt. Express 31, 32434 (2023).
  * Wilson and Bartels [2012] J. W. Wilson and R. A. Bartels, JOSA B 29, 1875 (2012).
  * Bernard  _et al._ [2023] O. Bernard, A. Kraxner, A. Boukhayma, J. A. Squier, C. Enz, and Y. Bellouard, Optica 10, 774 (2023).
  * Ward and New [1969] J. Ward and G. New, Physical Review 185, 57 (1969).
  * Masihzadeh  _et al._ [2009] O. Masihzadeh, P. Schlup, and R. A. Bartels, Optics letters 34, 1240 (2009).
  * Field  _et al._ [2016] J. J. Field, K. A. Wernsing, S. R. Domingue, A. M. Allende Motz, K. F. DeLuca, D. H. Levi, J. G. DeLuca, M. D. Young, J. A. Squier, and R. A. Bartels, Proceedings of the National Academy of Sciences 113, 6605 (2016).
  * Débarre  _et al._ [2005] D. Débarre, W. Supatto, and E. Beaurepaire, Optics letters 30, 2134 (2005).
  * Olivier  _et al._ [2010] N. Olivier, M. A. Luengo-Oroz, L. Duloquin, E. Faure, T. Savy, I. Veilleux, X. Solinas, D. Débarre, P. Bourgine, A. Santos, _et al._ , Science 329, 967 (2010).
  * Woodward  _et al._ [2016] R. Woodward, R. Murray, C. Phelan, R. De Oliveira, T. Runcorn, E. Kelleher, S. Li, E. De Oliveira, G. Fechine, G. Eda, _et al._ , 2D Materials 4, 011006 (2016).
  * Echarri  _et al._ [2018] A. R. Echarri, J. D. Cox, R. Yu, and F. J. Garcia de Abajo, ACS Photonics 5, 1521 (2018).
  * Clay  _et al._ [2006] G. O. Clay, A. C. Millard, C. B. Schaffer, J. Aus-der Au, P. S. Tsai, J. A. Squier, and D. Kleinfeld, JOSA B 23, 932 (2006).
  * Young  _et al._ [2015] M. D. Young, J. J. Field, K. E. Sheetz, R. A. Bartels, and J. Squier, Advances in optics and photonics 7, 276 (2015).
  * Popoff  _et al._ [2010] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, Physical Review Letters 104, 10.1103/physrevlett.104.100601 (2010).
  * May  _et al._ [2021] M. A. May, N. Barré, K. K. Kummer, M. Kress, M. Ritsch-Marte, and A. Jesacher, Nature Communications 12, 10.1038/s41467-021-24666-9 (2021).

Supplemental Materials: Synthetic aperture holographic third harmonic
generation microscopy

###  VII.1 Experimental setup

A home-built Yb:fiber amplifier system produces $<$80 fs pulses centered at
1065 nm at a repetition rate of 40 MHz with 1 W of power. The output light
from the amplifier is sent to a Martinez style compressor to pre-compensate
for dispersion in the THG holographic microscope. From the compressor, the
beam passes through an isolator and proceeds to the microscope section of the
layout, shown in Fig. S1. A waveplate and polarizing beam splitter cube
directs the beam into two paths, one to be used for the signal and one for the
reference. The signal beam path is focused on the sample with quasi-widefield
illumination by underfilling the back aperture of an aspheric objective lens
(Newport, 0.5NA, 20x) with a fill ratio of 0.378. THG signal is collected in a
transmission geometry by a microscope objective (Optosigma, Near UV, 0.45 NA)
and reimaged to the camera (Teledyne, e2v EMCCD).

Figure S1: Experimental layout of the THG Holographic microscope for
collection in a transmission geometry. The nonlinear crystal used for the
reference is TiO2. The red line is the fundamental light and the purple line
is THG.

Meanwhile, after the polarizing beam splitter cube, the reference beam path is
focused on the distal side of a TiO2 single crystal to generate the THG
reference field centered at 355 nm. The reference THG is collimated, sent
through a mechanical delay stage for timing control, and is directed to the
camera to interfere with the signal. A 4f configuration is used to expand the
THG reference beam size to fill the camera chip and control the off-axis tilt
angle relative to the signal arm. Both signal and reference arms pass through
filters to select for the THG wavelength of 355 nm and polarizers to ensure
similar linear polarization for optimal interference. THG holographic images
are recorded using LightField® software package provided with the Teledyne
camera.

###  VII.2 THG holography and field estimation

The THG hologram is formed through interference between a signal field,
$E_{\rm sig}(\mathbf{x}_{m})$, and a reference field, $E_{\rm
ref}(\mathbf{x}_{m})$, and the intensity of that interference, where
uninteresting constants of proportionality are suppressed, is given by

$I_{\rm h}(\mathbf{x}_{m})=|E_{\rm ref}(\mathbf{x}_{m})+E_{\rm
sig}(\mathbf{x}_{m})|^{2}=I_{\rm ref}(\mathbf{x}_{m})+I_{\rm
sig}(\mathbf{x}_{m})+\tilde{I}_{\rm
h}(\mathbf{x}_{m})\,e^{i\,\phi_{r}}+\tilde{I}^{*}_{\rm
h}(\mathbf{x}_{m})\,e^{-i\,\phi_{r}}.$ (S1)

Here, we define the reference intensity as $I_{\rm
ref}(\mathbf{x}_{m})=|E_{\rm ref}(\mathbf{x}_{m})|^{2}$, the signal intensity
as $I_{\rm sig}(\mathbf{x}_{m})=|E_{\rm sig}(\mathbf{x}_{m})|^{2}$, and the
pseudo intensity, from which we estimate the complex THG signal field is
$\tilde{I}_{\rm h}(\mathbf{x}_{m})=E_{\rm ref}(\mathbf{x}_{m})\,E_{\rm
sig}^{*}(\mathbf{x}_{m})$. The signal and reference fields may be separated
into the amplitude, $A$, and phase, $\phi$, as $E_{\rm
sig}(\mathbf{x}_{m})=A_{\rm sig}(\mathbf{x}_{m})\,\exp[i\,\phi_{\rm
sig}(\mathbf{x}_{m})]$ and $E_{\rm ref}(\mathbf{x}_{m})=A_{\rm
ref}(\mathbf{x}_{m})\,\exp[i\,\phi_{\rm ref}(\mathbf{x}_{m})]$.

The intensity pattern formed during the recording of the hologram is an
interferogram with phase information encoded as intensity modulations in the
phase difference, $\Delta\phi=\phi_{\rm sig}(\mathbf{x}_{m})-\phi_{\rm
ref}(\mathbf{x}_{m})+\phi_{r}$. Holography relies on the preparation of a
well-known (ideally) reference field. In our case, we use off-axis holography
with a uniform reference field, where the amplitude, $A_{\rm
sig}(\mathbf{x}_{m})=A_{0}$ is a constant and the phase is that of a tilted
plane wave, $\phi_{\rm ref}(\mathbf{x}_{m})={{\bf k}}_{r}\cdot{{\bf x}}_{m}$.
Here ${{\bf k}}_{r}$ is the wavevector of the reference wave relative to the
optical axis, along $\hat{z}$, the $z$-direction, of the holographic imaging
system.

We have included an additional random phase, $\phi_{r}$, that arises from
relative phase differences accumulated in the holographic recording process.
Over a set of holograms, one may see that $\phi_{r}$ is a random variable
uniformly distributed over $[0,2\,\pi]$. In conventional digital holography,
this phase offset is not problematic, as the overall image amplitude and phase
is not affected by such a uniform phase. However, in our case of synthetic
spatial aperture holography, this random phase is extremely detrimental. One
of the steps of our aberration correction algorithm estimates and corrects for
this phase.

With a well-characterized reference field, i.e., with a known $E_{\rm
ref}(\mathbf{x}_{m})$, we can estimate the complex-valued signal electric
field. In this case, this estimated field is that of the THG field that we
have imaged to the camera. The process for estimating the THG signal field
follows standard protocols [23, 24, 25] and is illustrated in Fig. S2. The
estimation of $E_{\rm sig}(\mathbf{x}_{m})$ begins by taking the Fast Fourier
Transform (FFT) of the THG hologram data that is described mathematically as
Eq. S1, i.e., ${{\mathcal{F}\\{I_{\rm h}(\mathbf{x}_{m})\\}}}$. Because we
have implemented off-axis holography, the terms in the FFT of Eq. S1 separate
into the first two terms, ${{\mathcal{F}\\{I_{\rm ref}(\mathbf{x}_{m})\\}}}$
and ${{\mathcal{F}\\{I_{\rm sig}(\mathbf{x}_{m})\\}}}$, appearing centered at
the zero, i.e., “dc”, spatial frequency, and ${{\mathcal{F}\\{\tilde{I}_{\rm
h}(\mathbf{x}_{m})\\}}}$ centered at ${{\bf k}}_{r}$, while
${{\mathcal{F}\\{\tilde{I}^{*}_{\rm h}(\mathbf{x}_{m})\\}}}$ is centered at
$-{{\bf k}}_{r}$.

The standard algorithm is to multiply the sideband centered at ${{\bf k}}_{r}$
by a filter to isolate that component and compute the inverse FFT (IFFT) to
recover an estimate of $E_{\rm sig}(\mathbf{x}_{m})$. This field is corrupted
by two factors: aberrations and amplitude variations in $E_{\rm
ref}(\mathbf{x}_{m})$ and by aberrations in the optical system that images the
scattered THG field from the specimen to the camera. With our algorithm, we
estimate and correct for these aberrations. As noted in the article, random
fluctuations in $\phi_{r}$ can degrade our signal to noise ratio (SNR). SNR is
boosted by rapidly acquiring a stack of THG holograms, following the THG field
estimation procedure, stacking each THG field into the column of a matrix,
then finally estimating the average THG field from the set of holograms from
the dominate singular vector of the matrix. The phase of the dominant singular
vector now contains an offset phase taken from the average of the set of THG
fields for that particular scan position in the sample,
$\mathbf{x}^{(j)}_{s}$. This process is repeated for many sample positions,
$\mathbf{x}^{(j)}_{s}$’s to build a synthetic spatial aperture, enabling
extended FOV synthetic aperture microscopy.

Figure S2: Reconstruction of the THG field. The hologram is a recorded
interference pattern on the camera. The reconstruction is performed in MATLAB.
The circular filter has a diameter of the approximate frequency support of the
collection objective ($2$NA$/\lambda$). The scale bar is $7\mu$m.

Sample scanning is done by translating the stage in steps of 2.78 $\mu$m,
equivalent to a spatial shift of 10 pixels across the camera and slightly less
than half of the illumination spot size. The recorded holograms are shifted to
correspond to the sample scan position, then flattened into vectors and
stacked into a matrix, referred here as the transmission matrix. The columns
of the transmission matrix are the different scanned positions taken across
the sample plane, and the rows are the pixel positions. Adjacent widefield THG
holograms within the transmission matrix have redundant spatial information
because adjacent scan positions overlap. If holograms have offset phases from
one another, the fields will destructively interfere with one another when the
complex fields are combined.

###  VII.3 Imaging model and definition of the THG field transmission matrix

The experimental setup establishes the imaging system that describes our
experimental system. To properly correct for aberrations acquired in the
system, we first build a model of the THG synthetic aperture imaging system
that we then exploit in our aberration correction algorithm. Phase distortions
from aberrations jeopardize the reconstruction of the field across the FOV.
These aberrations originate from random offset phases between each hologram
measurement, along with imaging system aberrations imposed by the illumination
and image collection optics, and aberrations introduced by the specimen. Our
goal is to estimate and correct for the sources of these phase distortions in
order to produce an estimate of the complex-valued THG scattered field. To
estimate these aberrations, we adapt our recently introduced method of
aberration-free synthetic spatial frequency aperture SHG field imaging [31]
for spatially scanned extended aperture field imaging for THG.

As noted above, we obtain a scaled estimate of the scattered THG signal field,
$E_{\rm sig}(\mathbf{x}_{m})$, that is incident on the camera. This signal
field is generated through THG scattering produced by the quasi wide-field
illumination of the sample by the fundamental beam, $E_{1}({{\bf x}})$, that
is a pulsed laser beam oscillating at a center frequency $\omega_{1}$. The THG
field, $E_{3}(\mathbf{x}_{m})$, generated in the sample region and referenced
to the object plane conjugate to the camera surface, is imaged to the camera
surface with a 4-f microscope imaging system, which is well described by an
isoplanatic image model.

The model for the system begins with the fundamental illumination field
incident, $E_{1}({{\bf x}})$, on the spatial coordinates of the object ${{\bf
x}}$. This fundamental field is produced by a focusing optic with an input
fundamental field, $E_{i}(\mathbf{u}_{i})$, presented to the back focal plane
of the lens with coordinates specified by the input spatial frequencies
$\mathbf{u}_{i}$. The illuminating fundamental field is present in the back
focal plane of the illumination lens and is written as

$E_{1}({{\bf
x}})=\int\,E_{i}(\mathbf{u}_{i})\,P_{i}(\mathbf{u}_{i})\,e^{-i\,2\,\pi\,{{\bf
x}}\cdot\mathbf{u}_{i}}\,d^{2}\mathbf{u}_{i}.$ (S2)

Here we have included the possibility of aberrations in the illumination field
that are modeled as a pupil phase, $\phi_{i}(\mathbf{u}_{i})$, of the
illumination lens with the pupil
$P_{i}(\mathbf{u}_{i})=\exp(i\,\phi_{i}(\mathbf{u}_{i}))$, and we have assumed
that the incident field underfills the input lens pupil. In addition, we have
suppressed un-important scaling factors in front of the integral. We retain
this suppression of unimpactful constants in all of the expressions in this
derivation.

In the THG scattering process, the cube of the fundamental field drives a
nonlinear source term proportional to $\chi^{(3)}({{\bf x}})\,E_{1}^{3}({{\bf
x}})$. For a sample with a thickness less than the phase mismatch length,
$\lambda_{3}/(n_{3}-n_{1})$, the image of the radiated field produced in a 4-f
microscope can be written as

$E_{3}(\mathbf{x}_{m})=\int H_{o}(\mathbf{x}_{m}+{{\bf x}})\,\chi^{(3)}({{\bf
x}})\,E_{1}^{3}({{\bf x}})\,d^{2}{{\bf x}}\rightarrow E_{\rm
sig}(\mathbf{x}_{m}),$ (S3)

where $\mathbf{x}_{m}$ represents the spatial coordinates of the plane on the
camera detector surface. The coherent spread function, $H_{o}({{\bf x}})$, is
the inverse Fourier transform of the complex-valued pupil function
$P_{o}(\mathbf{u}_{o})=\mathcal{X}_{o}(\mathbf{u}_{o})\,\exp(i\phi_{o}(\mathbf{u}_{o}))$
in the output isoplanatic imaging pupil of the THG microscope. This pupil
consists of a support defined by a circular characteristic function,
$\mathcal{X}_{o}(\mathbf{u}_{o})$, centered on the origin of the output
spatial frequency coordinate $\mathbf{u}_{o}$ with a radius,
$\mathrm{NA}/\lambda_{3}$, defined by the THG wavelength, $\lambda_{3}$, and
the imaging system numerical aperture, NA. For a thin object, the THG field is
directly proportional to the product $E_{3}(\mathbf{x}_{m})=\chi^{(3)}({{\bf
x}})\,E_{1}^{3}({{\bf x}})$, and thus the desired object nonlinear
susceptibility, $\chi^{(3)}({{\bf x}})$, may be obtained from a set of $E_{\rm
sig}(\mathbf{x}_{m})$ measurements that are extracted from the THG hologram.
this coherent spread function for the THG 4-f imaging system is given by

$H_{o}({{\bf x}})=\int P_{o}(\mathbf{u}_{o})\,e^{-i\,2\,\pi\,{{\bf
x}}\cdot\mathbf{u}_{o}}\,d^{2}\mathbf{u}_{o}.$ (S4)

This measured THG field represents our signal field for each measurement.

In this work, we build a synthetic aperture measurement of the specimen by
taking a sequence of THG field measurements, as estimated with a holographic
microscopy, for a number of measurements where the object is translated
laterally by $\mathbf{x}_{s}$, leading to the set of data
$E_{3}(\mathbf{x}_{m};\mathbf{x}_{s})=E_{3}({{\bf x}}-\mathbf{x}_{s})$. The
shifted positions are a discrete set with the jth position denoted by
$E^{(j)}_{\rm sig}(\mathbf{x}_{m})\equiv
E_{3}(\mathbf{x}_{m};\mathbf{x}_{s}^{(j)})=E_{3}({{\bf
x}}-\mathbf{x}_{s}^{(j)})$. Adding the shifted data to Eq. S3, we obtain the
expression

$E_{3}(\mathbf{x}_{m};\mathbf{x}_{s})=\int H_{o}(\mathbf{x}_{m}+{{\bf
x}})\,\chi^{(3)}({{\bf x}}-\mathbf{x}_{s})\,E_{1}^{3}({{\bf x}})\,d^{2}{{\bf
x}}.$ (S5)

Making the variable substitution ${{\bf x}}^{\prime}={{\bf
x}}-\mathbf{x}_{s}$, we may write

$E_{3}(\mathbf{x}_{m};\mathbf{x}_{s})=\int
H_{o}(\mathbf{x}_{m}+\mathbf{x}_{s}+{{\bf x}}^{\prime})\,\chi^{(3)}({{\bf
x}}^{\prime})\,E_{1}^{3}(\mathbf{x}_{s}+{{\bf x}}^{\prime})\,d^{2}{{\bf
x}}^{\prime}.$ (S6)

The fields are mapped to a synthetic spatial aperture by defining the
coordinate $\mathbf{x}_{o}=\mathbf{x}_{m}+\mathbf{x}^{(j)}_{s}$, so that we
may write

$E^{(j)}_{\rm sig}(\mathbf{x}_{o})=\int H_{o}(\mathbf{x}_{o}+{{\bf
x}}^{\prime})\,\chi^{(3)}({{\bf
x}}^{\prime})\,E_{1}^{3}(\mathbf{x}^{(j)}_{s}+{{\bf x}}^{\prime})\,d^{2}{{\bf
x}}^{\prime}.$ (S7)

Unfortunately, the extracted THG signal fields, $E_{\rm sig}(\mathbf{x}_{m})$,
are corrupted by aberrations that manifest as phase distortions. These phase
distortions arise from the measurement phase differences, $\phi_{r}$, and the
phase distortions from the THG imaging system, $\phi_{o}(\mathbf{u}_{o})$, and
from the cube of the fundamental field, $E_{1}^{3}({{\bf x}})$ that drives the
THG source term in the sample. This fundamental field is also described with
an isoplanatic model, $E_{1}({{\bf x}})=\int
P_{i}(\mathbf{u}_{i})\,e^{-i\,2\,\pi\,\mathbf{u}_{i}\cdot\mathbf{x}_{i}}\,d^{2}\mathbf{x}_{i}$,
with
$P_{i}(\mathbf{u}_{i})=g(\mathbf{u}_{i})\,\exp(i\phi_{i}(\mathbf{u}_{i}))$,
where $g(\mathbf{u}_{i})$ is the input fundamental beam that is mildly focused
to produce the quasi-widefield fundamental illumination beam and
$\phi_{i}(\mathbf{u}_{i})$ are pupil-plane aberrations imposed onto the
illumination beam.

Estimation of these aberrations requires a set of THG fields, $E^{(j)}_{\rm
sig}(\mathbf{x}_{m})$, that exhibit spatial overlap, and thus redundancy of
the underlying $\chi^{(3)}({{\bf x}})$ that we wish to extract. These
measurements are indexed by $j$ in which the object is translated in the plane
by $\mathbf{x}^{(j)}_{s}$.

####  VII.3.1 Transmission matrix

The transmission matrix, $T(\mathbf{x}_{o},\mathbf{x}_{i})$, is an operator
used in linear optics [43] that provides a map from an input field,
$E_{i}(\mathbf{x}_{i})$, with coordinates $\mathbf{x}_{i}$ denoting the input
plane where a field is incident on an optics system. The field exiting the
system, $E_{o}(\mathbf{x}_{o})$, at an output plane with coordinates
$\mathbf{x}_{o}$ reads

$E_{o}(\mathbf{x}_{o})=\int
T(\mathbf{x}_{o},\mathbf{x}_{i})\,E_{i}(\mathbf{x}_{i})\,d^{2}\mathbf{x}_{i}.$
(S8)

In a linear system, the transmission operator for a system with a planar
scattering object, described by $\chi({{\bf x}})$, is defined by

$T(\mathbf{x}_{o},\mathbf{x}_{i})=\int H_{o}(\mathbf{x}_{o},{{\bf
x}})\,\chi({{\bf x}})\,H_{i}({{\bf x}},\mathbf{x}_{i})d^{2}{{\bf x}}.$ (S9)

The Green’s function for illumination at the input plane, $\mathbf{x}_{i}$,
and propagation to the scattering plane at ${{\bf x}}$, is defined as
$H_{i}({{\bf x}},\mathbf{x}_{i})$. Similarly the exiting Green’s function
mapping the field from ${{\bf x}}$ to $\mathbf{x}_{o}$ is denoted
$H_{o}(\mathbf{x}_{o},{{\bf x}})$. In the case of an isoplanatic (i.e., a
shift invariant system), the transmission operator becomes

$T(\mathbf{x}_{o},\mathbf{x}_{i})=\int H_{o}(\mathbf{x}_{o}+{{\bf
x}})\,\chi({{\bf x}})\,H_{i}({{\bf x}}+\mathbf{x}_{i})d^{2}{{\bf x}}.$ (S10)

Analogously to the case of second harmonic generation holography [31], we
define a transmission matrix for THG with an isoplanatic imaging system as

$T(\mathbf{x}_{o},\mathbf{x}_{i})=\int H_{o}(\mathbf{x}_{o}+{{\bf
x}})\,\chi^{(3)}({{\bf x}})\,H_{i}({{\bf x}}+\mathbf{x}_{i})d^{2}{{\bf x}}.$
(S11)

Comparing this expression to Eq. S7, we see that they are identical if we make
the following associations: $H_{i}({{\bf x}}+\mathbf{x}_{i})\rightarrow
E_{1}^{3}(\mathbf{x}^{(j)}_{s}+{{\bf x}}^{\prime})$, $\chi^{(3)}({{\bf
x}})\rightarrow\chi^{(3)}({{\bf x}}^{\prime})$, and
$H_{o}(\mathbf{x}_{o}+{{\bf x}})\rightarrow H_{o}(\mathbf{x}_{o}+{{\bf
x}}^{\prime})$. Now, we may write the transmission operator as

$T(\mathbf{x}_{o},\mathbf{x}_{s})=\int H_{o}(\mathbf{x}_{o}+{{\bf
x}})\,\chi^{(3)}({{\bf x}})\,E_{1}^{3}(\mathbf{x}_{s}+{{\bf x}})d^{2}{{\bf
x}},$ (S12)

and where we note that we can exchange $\mathbf{x}_{i}$ and $\mathbf{x}_{s}$.

At this point, it is important to note that this definition of a transmission
operator is not a general transmission operator as in the case of a linear
transformation, but rather a sub-space of a higher-order transmission tensor.
However, we can define this effective THG transmission operator that we will
use for estimation of the complex THG field, which is proportional to
$\chi^{(3)}({{\bf x}})$. Once defined over a set of the discrete spatial
coordinates, the transmission operator is then represented as a transmission
matrix as discussed in the article.

####  VII.3.2 Estimation of the synthetic THG field

Once we have obtained a corrected form of the transmission operator matrix, we
are in a position to estimate the THG field. To understand how this proceeds,
we note that once the data are transformed into the transmission matrix, Eq.
S7, we have numerically mapped the stationary cube of the fundamental
illumination field, $E_{1}^{3}({{\bf x}})$, to an effective scanning
illumination. Then, the THG source generated by this incident cubed field
generates a THG field that we image with the coherent spread function
$H_{o}({{\bf x}})$. Let us consider an idealized case where the fundamental
field is a normally-incident, unity-amplitude plane wave, i.e., $E_{1}({{\bf
x}})=1$. Then, the THG transmission operator becomes

$T(\mathbf{x}_{o},\mathbf{x}_{s})=\int H_{o}(\mathbf{x}_{o}+{{\bf
x}})\,\chi^{(3)}({{\bf x}})d^{2}{{\bf x}}.$ (S13)

We see that in this approximation, we obtain a low-pass filtered image of
$\chi^{(3)}({{\bf x}})$. In the limiting case where $H_{o}({{\bf
x}})\rightarrow\delta({{\bf x}})$, then we simply have

$T(\mathbf{x}_{o},\mathbf{x}_{s})\approx\chi^{(3)}(-\mathbf{x}_{o}).$ (S14)

However, because the illumination field has a finite spatial extent, the
estimate of the THG field is obtained by integrating over the scan coordinate,
$\mathbf{x}_{s}$, leading to

$\hat{E}_{3}(\mathbf{x}_{m})=\int
T(\mathbf{x}_{o},\mathbf{x}_{s})d^{2}\mathbf{x}_{s}.$ (S15)

Applying the definition of this estimator to Eq. S7, we obtain

$\hat{E}_{3}(\mathbf{x}_{m})=\int H_{o}(\mathbf{x}_{o}+{{\bf
x}})\,\chi^{(3)}({{\bf x}})\,F_{1}({{\bf x}})\,d^{2}{{\bf x}}.$ (S16)

where we have defined the synthetic illumination function

$F_{1}({{\bf x}})\equiv\int E_{1}^{3}(\mathbf{x}_{s}+{{\bf
x}})\,d^{2}\mathbf{x}_{s}.$ (S17)

This function gates the region of illumination of the sample. Under conditions
of ideal convolution, the $F_{1}({{\bf x}})\rightarrow\mathbbm{1}({{\bf x}})$,
indicating that the illumination function converges to a characteristic
function that defines the region of illumination in the specimen. The goal of
the aberration correction algorithm in the following section is to remove the
aberrations introduced by $E_{1}^{3}({{\bf x}})$ and
$H_{o}(\mathbf{x}_{o}+{{\bf x}})$. With the estimated THG field, we then
obtain the estimated intensity from $I_{\rm
THG}(\mathbf{x}_{o})=|\hat{E}_{3}(\mathbf{x}_{m})|^{2}$.

###  VII.4 Nonlinear distortion operator and aberration correction of THG
spatial synthetic aperture holography

Matrix-based approaches to optical imaging have been explored in recent years
to combat multiple scattering, wavefront distortions, and aberrations in
biological imaging.[43] We use the distortion operator approach inspired by
Badon et. al to correct for pupil plane phase distortions.[29] This approach
relies on the measurement of the distortion operator, $D$. Here, we show that
our experimental setup directly measures $D(\mathbf{x}_{m},\mathbf{x}_{s})$ in
the spatial-spatial domain, and we show how $D$ is related to the transmission
matrix. Here, we describe the THG nonlinear distortion operator (NDO) measured
in these experiments, which is a generalization of the NDO recently introduced
for SHG holographic imaging. [31] In the experiments, a set of THG field data
are acquired, one for each specimen location,
$\mathbf{x}_{s}=\mathbf{x}^{(j)}_{s}$, as given in Eq. S5. This expression is
the definition of the spatial-spatial THG distortion operator,
$D(\mathbf{x}_{m},\mathbf{x}_{s})$. As we saw in the previous section, the
transmission operator arises from a coordinate shift introduced to the data,
resulting in Eq. S7.

In the case of an isoplanatic system, aberrations are corrected by taking a
singular value decomposition of the distortion operator in the spatial-
frequency/spatial space, i.e., $D(\mathbf{u}_{o},\mathbf{x}_{s})$. [29] To
obtain this mixed-space operator, we transform the output spatial coordinate
to the conjugate spatial frequency coordinate, $\mathbf{u}_{o}$, using the
transform

$D(\mathbf{u}_{o},\mathbf{x}_{s})=\int
D(\mathbf{x}_{m},\mathbf{x}_{s})\,e^{i\,2\,\pi\,\mathbf{x}_{m}\cdot\mathbf{u}_{o}}\,d^{2}\mathbf{x}_{m}.$
(S18)

Using the definition of the spatial-spatial NDO in Eq. S5, we obtain

$D(\mathbf{u}_{o},\mathbf{x}_{s})=\iint H_{o}(\mathbf{x}_{m}+{{\bf
x}})\,e^{i\,2\,\pi\,\mathbf{x}_{m}\cdot\mathbf{u}_{o}}\,\chi^{(3)}({{\bf
x}}-\mathbf{x}_{s})\,E_{1}^{3}({{\bf x}})\,d^{2}{{\bf
x}}\,d^{2}\mathbf{x}_{m}.$ (S19)

Making the variable substitution ${{\bf x}}^{\prime}={{\bf x}}-\mathbf{x}_{s}$
and using $\mathbf{x}_{o}=\mathbf{x}_{s}+\mathbf{x}_{m}$ leads to

$D(\mathbf{u}_{o},\mathbf{x}_{s})=e^{-i\,2\,\pi\,\mathbf{x}_{s}\cdot\mathbf{u}_{o}}\,\int
e^{i\,2\,\pi\,\mathbf{x}_{o}\cdot\mathbf{u}_{o}}\,\left(\int
H_{o}(\mathbf{x}_{o}+{{\bf x}}^{\prime})\,\chi^{(3)}({{\bf
x}}^{\prime})\,E_{1}^{3}(\mathbf{x}_{s}+{{\bf x}}^{\prime})\,d^{2}{{\bf
x}}^{\prime}\right)\,d^{2}\mathbf{x}_{o}.$ (S20)

Identifying the quantity in the parentheses as the transmission operator,
$T(\mathbf{x}_{o},\mathbf{x}_{s})$, we obtain the standard definition of the
distortion operator [29, 31]

$D(\mathbf{u}_{o},\mathbf{x}_{s})=e^{-i\,2\,\pi\,\mathbf{x}_{s}\cdot\mathbf{u}_{o}}\,T(\mathbf{u}_{o},\mathbf{x}_{s}),$
(S21)

where we have defined the spatial-frequency/spatial transmission operator as

$T(\mathbf{u}_{o},\mathbf{x}_{s})=\int
e^{i\,2\,\pi\,\mathbf{x}_{o}\cdot\mathbf{u}_{o}}\,T(\mathbf{x}_{o},\mathbf{x}_{s})\,d^{2}\mathbf{x}_{o}.$
(S22)

Applying our assumption of an isoplanatic system, we may write

$T(\mathbf{u}_{o},\mathbf{x}_{s})=P_{o}(\mathbf{u}_{o})\,\int
e^{-i\,2\,\pi\,{{\bf x}}^{\prime}\cdot\mathbf{u}_{o}}\chi^{(3)}({{\bf
x}}^{\prime})\,E_{1}^{3}(\mathbf{x}_{s}+{{\bf x}}^{\prime})\,d^{2}{{\bf
x}}^{\prime}.$ (S23)

We may then write the distortion operator in the isoplanatic case as

$D(\mathbf{u}_{o},\mathbf{x}_{s})=P_{o}(\mathbf{u}_{o})\,\int\chi^{(3)}({{\bf
x}})\,E_{1}^{3}(\mathbf{x}_{s}+{{\bf
x}})\,e^{-i\,2\,\pi\,(\mathbf{x}_{s}+{{\bf
x}})\cdot\mathbf{u}_{o}}\,d^{2}{{\bf x}}.$ (S24)

In order to understand how and when the aberrations can be estimated, we
consider the correlation

$DD^{*}(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})=P_{o}(\mathbf{u}_{o})\,\Gamma(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})\,P_{o}^{*}(\mathbf{u}_{o}^{\prime}),$
(S25)

where

$\Gamma(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})=\iint\chi^{(3)}({{\bf
x}})\,\chi^{(3)*}({{\bf
x}}^{\prime})\,K(\mathbf{u}_{o}^{\prime}-\mathbf{u}_{o},{{\bf
x}}^{\prime}-{{\bf x}})\,e^{-i\,2\,\pi\,{{\bf
x}}\cdot\mathbf{u}_{o}}\,e^{i\,2\,\pi\,{{\bf
x}}^{\prime}\cdot\mathbf{u}_{o}^{\prime}}\,d^{2}{{\bf x}}\,d^{2}{{\bf
x}}^{\prime},$ (S26)

and with the illumination correlation function defined as

$K(\Delta{{\bf u}},\Delta{{\bf x}})=\int
E_{1}^{3}(\mathbf{x}_{s})\,E_{1}^{3*}(\mathbf{x}_{s}+\Delta{{\bf
x}})\,e^{-i\,2\,\pi\,\mathbf{x}_{s}\cdot\Delta{{\bf
u}}}\,d^{2}\mathbf{x}_{s},$ (S27)

where we have identified $\Delta{{\bf x}}={{\bf x}}^{\prime}-{{\bf x}}$ and
$\Delta{{\bf u}}=\mathbf{u}_{o}-\mathbf{u}_{o}^{\prime}$. The left singular
vectors of the SVD of the distortion operator $D$ will then yield the
eigenvectors of the correlation $DD^{*}$.

The form of $\Gamma$ depends on two important characteristic lengths: the
spatial support of the illumination, $\ell_{1}$, and a typical length for
$\chi^{(3)}({{\bf x}})$, $\ell_{0}$, which models either its spatial support
or its scale of oscillation, whichever is shortest. There are various scenarii
depending on which length is shortest, and we suppose as a start that
$\ell_{0}$ corresponds to the spatial support of $\chi^{(3)}$. Here,
$E_{1}^{3}(\mathbf{x}_{s})$ is a function of $\mathbf{x}_{s}/\ell_{1}$ and
$\chi^{(3)}({{\bf x}})$ is a function of ${{\bf x}}/\ell_{0}$. When
$\ell_{0}\ll\ell_{1}$, then the $\Delta x$ can be ignored in $K$ and replaced
by $0$, so that

$\Gamma(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})=\hat{\chi}^{(3)}(\mathbf{u}_{o})\,K(\mathbf{u}_{o}^{\prime}-\mathbf{u}_{o},0)\,\hat{\chi}^{(3)}(\mathbf{u}_{o}^{\prime})^{*}$
(S28)

and then

$DD^{*}(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})=P_{o}(\mathbf{u}_{o})\,\hat{\chi}^{(3)}(\mathbf{u}_{o})\,K(\mathbf{u}_{o}^{\prime}-\mathbf{u}_{o},0)\,\hat{\chi}^{(3)}(\mathbf{u}_{o}^{\prime})^{*}\,P_{o}^{*}(\mathbf{u}_{o}^{\prime}).$
(S29)

We need to introduce the scale of oscillation of the phase function
$\phi_{o}$, denoted $u_{o}$, to extract information from the equation above.
We suppose then that $\phi_{o}(\mathbf{u}_{o})$ is a function of
$\mathbf{u}_{o}/u_{o}$. Since $E_{1}^{3}$ is localized at the scale
$\ell_{1}$, $K$ localizes $\mathbf{u}_{o}$ around $\mathbf{u}_{o}^{\prime}$ at
the scale $\ell_{1}^{-1}$. Hence, for a nearly widefield imaging scenario when
$\ell_{1}$ is large compared to $u_{o}^{-1}$,
$K(\mathbf{u}_{o}^{\prime}-\mathbf{u}_{o},0)$ is narrow, and we can replace
$\hat{\chi}^{(3)}(\mathbf{u}_{o}^{\prime})^{*}\,P_{o}^{*}(\mathbf{u}_{o}^{\prime})$
by $\hat{\chi}^{(3)}(\mathbf{u}_{o})^{*}\,P_{o}^{*}(\mathbf{u}_{o})$,
resulting in

$DD^{*}(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})=|\hat{\chi}^{(3)}(\mathbf{u}_{o})|^{2}\,K(\mathbf{u}_{o}^{\prime}-\mathbf{u}_{o},0),$
(S30)

making $DD^{*}(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})$ nearly diagonal. In
that case it is not possible to extract the aberrations from the distortion
operator. We need then the illumination to be sufficiently focused, namely
satisfying $u_{o}^{-1}\geq\ell_{1}$ (but still satisfying
$\ell_{0}\ll\ell_{1}$), to extract the aberrations. In that case,
$DD^{*}(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})$ is more extended around the
diagonal and we have

$\displaystyle\int
DD^{*}(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})P_{o}(\mathbf{u}_{o}^{\prime})d^{2}\mathbf{u}_{o}^{\prime}$
$\displaystyle=$ $\displaystyle
P_{o}(\mathbf{u}_{o})\,\hat{\chi}^{(3)}(\mathbf{u}_{o})\,\int
K(\mathbf{u}_{o}^{\prime}-\mathbf{u}_{o},0)\,\hat{\chi}^{(3)}(\mathbf{u}_{o}^{\prime})^{*}\,d^{2}\mathbf{u}_{o}^{\prime}$
$\displaystyle\simeq$ $\displaystyle
P_{o}(\mathbf{u}_{o})\,|\hat{\chi}^{(3)}(\mathbf{u}_{o})|^{2}\int
K(\mathbf{u}_{o}^{\prime},0)\,d^{2}\mathbf{u}_{o}^{\prime},$

where we used the string of inequalities $u_{o}^{-1}\geq\ell_{1}\gg\ell_{0}$
and the change of variables
$\mathbf{u}_{o}^{\prime}-\mathbf{u}_{o}\to\mathbf{u}_{o}^{\prime}$. Since
$\hat{\chi}^{(3)}(\mathbf{u}_{o})$ varies at the scale $\ell_{1}^{-1}$, which
is much greater than that of $P_{o}(\mathbf{u}_{o})$, the function
$\hat{\chi}^{(3)}(\mathbf{u}_{o})$ is nearly constant compared to
$P_{o}(\mathbf{u}_{o})$ and as a consequence $P_{o}(\mathbf{u}_{o})$ can be
separated from $|\hat{\chi}^{(3)}(\mathbf{u}_{o})|^{2}$ and is approximately
an eigenvector. Now, for $P_{o}(\mathbf{u}_{o})$ to be estimated for all
$\mathbf{u}_{o}$, the largest wavenumber that can be measured by the
experimental system, denoted $u_{m}$, has to be much less than $\ell_{0}^{-1}$
since $|\hat{\chi}^{(3)}(\mathbf{u}_{o})|^{2}$ becomes small when the length
of $\mathbf{u}_{o}$ is greater than $\ell_{0}^{-1}$. In that scenario,
$DD^{*}$ is approximately of rank one and the leading eigenvector of $DD^{*}$
gives an estimation of $P_{o}(\mathbf{u}_{o})$.

In the opposite situation where $\ell_{1}\ll\ell_{0}$, the analysis is similar
as above with $\chi^{(3)}$ and $E_{1}^{3}$ swapped. Indeed,
$\mathbf{x^{\prime}}$ is now localized around $\mathbf{x}$ in
$\chi^{(3)*}({{\bf x}}^{\prime})$, which becomes $\chi^{(3)*}(\mathbf{x})$,
and changing variables as $\mathbf{x^{\prime}}\to\mathbf{x}-\mathbf{x}_{s}$,
we find

$\Gamma(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})=\left(\int|\chi^{(3)}(\mathbf{x})|^{2}e^{-i2\pi\mathbf{x}\cdot\Delta\mathbf{u}}d^{2}{{\bf
x}}\right)\,\int K(\mathbf{u}_{o}^{\prime}-\mathbf{u}_{o},{{\bf
x}}^{\prime}-\mathbf{x}_{s})\,e^{i2\pi\mathbf{u}_{0}^{\prime}\cdot(\mathbf{x}^{\prime}-\mathbf{x}_{s})}d^{2}{{\bf
x}}^{\prime}.$ (S31)

The first term on the right above localizes $\Delta\mathbf{u}$ around 0 at the
scale $\ell_{0}^{-1}$, while the integral of $K$ can be recast as

$\hat{E}_{1}^{3}(\mathbf{u}_{o})(\hat{E}_{1}^{3}(\mathbf{u}_{o}^{\prime}))^{*}$

where $\hat{E}_{1}^{3}$ is the Fourier transform of $E_{1}^{3}$. Based on the
analysis for the case $\ell_{1}\gg\ell_{0}$, the aberrations can then be
reconstructed provided $u_{o}^{-1}\geq\ell_{0}$, namely the support of
$\chi^{(3)}$ is sufficiently small compared to $u_{o}^{-1}$. The condition
$u_{m}\ll\ell_{1}^{-1}$ is also necessary to estimate $P_{o}(\mathbf{u}_{o})$
for all $\mathbf{u}_{o}$ since $\hat{E}_{1}^{3}(\mathbf{u}_{o})$ becomes small
when the length of $\mathbf{u}_{o}$ is greater than $\ell_{1}^{-1}$.

When $\ell_{0}$ corresponds to the scale of oscillation of $\chi^{(3)}$, we
rewrite $\chi^{(3)}$ as
$\chi^{(3)}(\mathbf{x})=\varphi(\mathbf{x}/\ell_{\varphi})e^{i2\pi\mathbf{u}_{\chi}\cdot\mathbf{x}}$
where $\varphi$ is an envelope function and the length of $\mathbf{u}_{\chi}$
is $\ell_{0}^{-1}$. Here $\ell_{0}\ll\ell_{\varphi}$, and we assume
$\ell_{\varphi}\gg\ell_{1}$. In that case, a short calculation shows that
$\Gamma$ is as in (S31) with $\chi^{(3)}$ replaced by $\varphi$ and the
integral of $K$ becomes

$\hat{E}_{1}^{3}(\mathbf{u}_{o}-\mathbf{u}_{\chi})(\hat{E}_{1}^{3}(\mathbf{u}_{o}^{\prime}-\mathbf{u}_{\chi}))^{*}.$

The correlation matrix is not too narrow when $u_{o}^{-1}\geq\ell_{\varphi}$.
The term $P_{o}(\mathbf{u}_{o})$ has then to be separated from
$\hat{E}_{1}^{3}(\mathbf{u}_{o}-\mathbf{u}_{\chi})$, which is possible when
$u_{o}^{-1}\gg\ell_{1}$, and $P_{o}(\mathbf{u}_{o})$ can be estimated for all
$\mathbf{u}_{o}$ provided $\ell_{1}^{-1}\gg u_{m}$ and the length of
$\mathbf{u}_{\chi}$, namely $\ell_{0}^{-1}$, is less than $u_{m}$.

In the scenario where $\ell_{0}$ and $\ell_{1}$ are comparable, the expression
of $\Gamma$ does not simplify, but the kernel $K$ still does localize
$\mathbf{u}_{o}$ around $\mathbf{u}_{o}^{\prime}$ at the scale
$\ell_{1}^{-1}$. Hence, as before, the aberrations can be reconstructed when
the illumination is focused enough so that $u_{o}^{-1}\geq\ell_{1}$.

The imaging system’s pupil phase distortions are estimated from the first left
singular value from the SVD of the mixed basis NDO,
D($\mathbf{u}_{o},\mathbf{x}_{s}$), which is the eigenvector of the
correlation of the NDO, $DD^{*}(\mathbf{u}_{o},\mathbf{u}_{o}^{\prime})$. This
is particularly the case in our setup when the illumination is stationary
relative to the sample, making our the detection scheme isoplanatic. An
example of the amplitude and phase of the SVD first left singular value of the
NDO is shown in Fig. S3.

Figure S3: The amplitude (left) and phase (right) of the SVD first left
singular value of the NDO. The images shown are in the spatial frequency
domain with the circular diameter being approximately the spatial frequency
support of the collection objective.

Our experimental setup allows us to directly measure the nonlinear distortion
operator, perform a Fourier transform to get the mixed basis NDO, and estimate
the pupil phase distortions and subsequently correct for those distortions.
However, as pointed out in the main text, correcting for the dominant phase
aberrations that occur from phase jumps in hologram-to-hologram measurements
is necessary before using the NDO aberration correction.

###  VII.5 Optimization algorithm for spatial phase estimation and correction

The phase of each hologram recorded from different spatial positions has a
random phase offset, $\phi_{r}(\mathbf{x}_{s})$, relative to each other due to
random phase fluctuations between measurements. This random phase difference
between the holograms of each $\mathbf{x}_{s}$ position in the transmission
matrix prevents a good reconstruction of the THG field across the FOV. It is
necessary to correct for these random phase fluctuations to coherently combine
the THG fields. A brute force optimization approach can correct the phase by
applying a piston phase shift to each reconstructed THG field in transmission
matrix and maximizing the integrated intensity of the transmission matrix.
Maximizing the intensity serves as a good parameter to optimize because when
the phase offset of each field is corrected, adjacent fields will
constructively combine and cause the intensity of the transmission matrix to
increase. However, this brute force approach becomes computationally time
consuming for large data sets. Therefore, we utilize a recently created
adaptive phase algorithm for optimizing the phase matching of reconstructed
THG fields in the transmission matrix.

Dynamic adaptive scattering compensation holography (DASH) was originally
designed for applying phase masks to a spatial light modulator (SLM) to
optimize beam focusing through scattering media [44]. The algorithm finds the
optimal phase mask to project onto the SLM by continuously updating and
improving the mask with every iteration. Here, we modify the DASH algorithm to
apply a basis of digital phase masks to operate on the transmission matrix,
rather than an SLM. The DASH phase optimized transmission matrix is defined as
the original transmission matrix by the phase correction,

$T^{c}(\mathbf{x}_{o},\mathbf{x}_{s})=T(\mathbf{x}_{o},\mathbf{x}_{s})\,\exp[-i\,\phi_{r}(\mathbf{x}_{s})].$
(S32)

The corrected image is then estimated by the computation

$I^{c}(\mathbf{x}_{o})=\left|\int
T^{c}(\mathbf{x}_{o},\mathbf{x}_{s})\,d^{2}\mathbf{x}_{s}\right|^{2}.$ (S33)

The phase correction pattern $\phi_{r}(\mathbf{x}_{s})$ calculated with our
modified DASH algorithm is defined by

$\phi_{r}(\mathbf{x}_{s})=\mathrm{angle}(C_{I,N}),$ (S34)

where $C_{I,N}$ is the final correction mask after the iterative optimization
and $N$ is the total number of recorded holograms, (i.e., the total number of
$\mathbf{x}_{s}$ positions). A modulated wavefront mode, $M_{n}$, a phase
step, $\phi_{p}$, and a reference field $C_{i,n}$ contribute to constructing
the phase pattern

$\Phi_{i,n,p}=\mathrm{angle}\left(\sqrt{1-f}\,\frac{C_{i,n}}{|C_{i,n}|}+\sqrt{f}\,e^{i(M_{n}+\phi_{p})}\right),$
(S35)

where the subscript $i$ is the iteration number, $n$ is the mode number, $p$
is the phase step which is set to be $\phi_{p}=p(2\pi/P)$. We chose $P=5$, as
recommended in the original DASH work, such that the
$p=[0,2\pi/5,4\pi/5,6\pi/5,8\pi/5]$. The value of $f=0.3$ is chosen under
recommendation of the authors of the original DASH work. The reference field
for the first iteration is set to zero. The phase pattern is then applied the
transmission matrix to maximize the cost function

$J=\int I^{c}(\mathbf{x}_{o})\,d^{2}\mathbf{x}_{o},$ (S36)

where

$I^{c}(\mathbf{x}_{o})=\left|\int
T(\mathbf{x}_{o},\mathbf{x}_{s})\,\exp\left(i\,\Phi_{i,n,p}\right)\,d^{2}\mathbf{x}_{s}\right|^{2},$
(S37)

The algorithm undergoes iterations to continually update the reference field

$C_{i,n}+_{1}=C_{i,n}+a_{n}\exp[i(M_{n})]$ (S38)

determined by an amplitude weighting function $a_{n}$. The weighting function
is defined as

$a_{n}=\sum_{p}\frac{\sqrt{J}\exp(i\phi_{p})}{P}$ (S39)

such that the phase correction mask is built from the basis set of modulated
wavefront modes $M_{n}$ and the weighting function gives modes that increase
the cost function amplitude, and thus, more contribution to the final phase
correction mask $C_{I,N}$. After several iterations (typically 7 iterations
for N=676 modes) the reference field converges to the final field $C_{I,N}$,
which maximizes $J$.

The result is a phase correction pattern that, when applied to the
transmission matrix, maximizes the intensity of the transmission matrix, and
thus corrects for the hologram-to-hologram phase distortions,
$\phi_{r}(\mathbf{x}_{s})$. The spatial overlap of each hologram provides the
constructive/destructive interference necessary to optimize. Without spatial
overlap, the phase of the holograms become independent of one another. The
modified DASH algorithm provides a fast-converging optimization ($<$10
minutes) for moderately sized data ($\sim$ 4GB).

Our modified DASH algorithm is used here for correcting phase differences
between the signal and reference field measured for each $\mathbf{x}_{s}$
position. However, the pupil plane distortions still exist. We correct pupil
plane phase aberrations through the distortion operator (DO) approach we have
adapted for the THG transmission matrix.