Synthetic aperture holographic third harmonic generation microscopy

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.


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. 45][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,18Particularly as suitable laser systems become more readily available, THG is increasingly playing an important role in biological imaging.
0][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 pointscanning 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.[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
]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 ref (x m ) and a signal field E sig (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 ∼ 7µm, or roughly 15× 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 TiO 2 window, where strong THG light is generated. 1,2,32The 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, ω 3 = 3 ω 1 that is three times the fundamental frequency ω 1 .These two fields interfere at the camera surface to form a holographic intensity pattern ).We numerically isolate the interferometric pseudo-intensity term Ĩh = E * ref (x m ) E sig (x m ) to estimate the complex THG image field using standard protocols. 23,24or 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 x (j) s produces a two-dimensional (2D) complex-valued THG field map, E (j) sig (x m ).An example set obtained by scanning x s is illustrated in Fig. 1.The spatial coordinates of the camera, defined by the pixel pitch, are mapped to measurement coordinates, 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 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) sig (x m ), in the data stack is flattened into a linear vector that is inserted into a column of a matrix, D(x m , x s ) as illustrated in Fig. 1.The column for each specimen scan position, x (j) s is set to match the lexicographical order of the 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). 31other 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) sig (x m ), for the j th sample position x (j) s with a mean phase φ 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.
The first correction in our algorithm corrects for the dominant phase distortion, φ r (x s ) that arises from the hologram measurement process.This phase is uniformly described over 2π, leading to speckle on an estimated image.The first step to correcting φ r (x s ) maps the spatial-spatial NDO to the THG transmission matrix T (x o , x s ).We built an optimization algorithm that maximizes the cost function, x o position of the THG field, a correction of φ r (x s ) maximizes the estimated object intensity at each spatial point by exploiting the fact that the spatial map is highly correlated along x o .This is possible because the data shares the spatial structure of the sample nonlinear susceptibility, χ (3) (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 φ r (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, φ o (u o ), by using the concept of the distortion operator 29 that is adapted for nonlinear imaging. 31To use the NDO for aberration correction, we map the "DASH corrected" spatial-spatial transmission matrix, T (x o , x s ), into the spatial frequency-spatial distortion operator D(u o , x s ).To do so, we first shift the spatial fields back from x o to x s to produce D(x m , x s ).Then each 2D output field for each x s is Fourier transformed to the spatial frequency space, so that the NDO provides a relationship between x s and the output spatial frequency coordinates 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(u o , 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 u o across the x s spatial positions.Then, an inverse Fourier transform of u o brings back the spatial-spatial NDO, D(x m , 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 (x o , 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 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. 23Focal propagation distances are estimated using a gradient sharpness criterion.
The first set of images shown in Fig. 2 are of a 6µm thick mouse tail cross section prepared on a glass slide (Triarch Incorporated).THG reconstructed images are of the epidural region   The final THG reconstructed image presented here is that of monolayer molybdenum disulfide (MoS 2 ) (Fig. 4).Monolayer MoS 2 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 MoS 2 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.

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 pointscanning microscopy except for a recent report on intensity widefield THG with a lowrepetition rate amplifier. 33The 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, χ (1) (x), and nonlinear, χ (3) (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 (χ (3) (x)).Because (χ (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. 34Conventional 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 χ (3) (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 2 0 ∼ 0.5 mm, where k 1 = n 1 ω 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, ∆n = n 3 − n 1 > 0 due to the dispersion of the RI.The collinear phase mismatch parameter is defined as ∆k = k 3 − 3k 1 = ∆n = n 3 − n 1 > 0.Here we have also written the THG wavenumber as k 3 = n 3 ω 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 χ (3) (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 (r), that is imprinted onto the THG field.In addition, spatial inhomogeneities in n 1 (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. 25cause the χ (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,38G phase variations can also arise from the physics of the third-order THG optical tensor.A complex-valued χ (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 MoS 2 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. Forinstance, previous spectroscopy measurements of THG scattering have revealed two-photon absorption resonances in oxy-and deoxy-hemoglobin.41 Because χ (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. 42In 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., Im{χ (3)   Meanwhile, after the polarizing beam splitter cube, the reference beam path is focused on the distal side of a TiO 2 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.

B. THG holography and field estimation
The THG hologram is formed through interference between a signal field, E sig (x m ), and a reference field, E ref (x m ), and the intensity of that interference, where uninteresting constants of proportionality are suppressed, is given by Here, we define the reference intensity as , the signal intensity as , and the pseudo intensity, from which we estimate the complex THG The signal and reference fields may be separated into the amplitude, A, and phase, φ, as The intensity pattern formed during the recording of the hologram is an interferogram with phase information encoded as intensity modulations in the phase difference, ∆φ = φ sig (x m ) − φ ref (x m ) + φ 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 sig (x m ) = A 0 is a constant and the phase is that of a tilted plane wave, Here k r is the wavevector of the reference wave relative to the optical axis, along ẑ, the z-direction, of the holographic imaging system.
We have included an additional random phase, φ r , that arises from relative phase differences accumulated in the holographic recording process.Over a set of holograms, one may see that φ r is a random variable uniformly distributed over [0, 2 π].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 ref (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 S23-S25 and is illustrated in Fig. S2.The estimation of E sig (x m ) begins by taking the Fast Fourier Transform (FFT) of the THG hologram data that is described mathematically as Eq.S1, i.e., F {I h (x m )}.Because we have implemented off-axis holography, the terms in the FFT of Eq.S1 separate into the first two terms, F {I ref (x m )} and F {I sig (x m )}, appearing centered at the zero, i.e., "dc", spatial frequency, and The standard algorithm is to multiply the sideband centered at k r by a filter to isolate that component and compute the inverse FFT (IFFT) to recover an estimate of E sig (x m ).
This field is corrupted by two factors: aberrations and amplitude variations in E ref (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 φ 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, x s .This process is repeated for many sample positions, x

C. 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 S31 for spatially scanned extended aperture field imaging for THG.
As noted above, we obtain a scaled estimate of the scattered THG signal field, E sig (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, , that is a pulsed laser beam oscillating at a center frequency ω 1 .The THG field, E 3 (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 (x), on the spatial coordinates of the object x.This fundamental field is produced by a focusing optic with an input fundamental field, E i (u i ), presented to the back focal plane of the lens with coordinates specified by the input spatial frequencies u i .The illuminating fundamental field is present in the back focal plane of the illumination lens and is written as Here we have included the possibility of aberrations in the illumination field that are modeled as a pupil phase, φ i (u i ), of the illumination lens with the pupil P i (u i ) = exp(i φ i (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 χ (3) (x) E 3 1 (x).For a sample with a thickness less than the phase mismatch length, λ 3 /(n 3 − n 1 ), the image of the radiated field produced in a 4-f microscope can be written as where x m represents the spatial coordinates of the plane on the camera detector surface.
The coherent spread function, H o (x), is the inverse Fourier transform of the complex-valued in the output isoplanatic imaging pupil of the THG microscope.This pupil consists of a support defined by a circular characteristic function, X o (u o ), centered on the origin of the output spatial frequency coordinate u o with a radius, NA/λ 3 , defined by the THG wavelength, λ 3 , and the imaging system numerical aperture, NA.For a thin object, the THG field is directly proportional to the product , and thus the desired object nonlinear susceptibility, χ (3) (x), may be obtained from a set of E sig (x m ) measurements that are extracted from the THG hologram.this coherent spread function for the THG 4-f imaging system is given by 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 x s , leading to the set of data E 3 (x m ; x s ) = E 3 (x − x s ).The shifted positions are a discrete set with the j th position denoted by s ).Adding the shifted data to Eq. S3, we obtain the expression Making the variable substitution x ′ = x − x s , we may write The fields are mapped to a synthetic spatial aperture by defining the coordinate s , so that we may write Unfortunately, the extracted THG signal fields, E sig (x m ), are corrupted by aberrations that manifest as phase distortions.These phase distortions arise from the measurement phase differences, φ r , and the phase distortions from the THG imaging system, φ o (u o ), and from the cube of the fundamental field, E 3 1 (x) that drives the THG source term in the sample.This fundamental field is also described with an isoplanatic model, E 1 (x) = P i (u i ) e −i 2 π u i •x i d 2 x i , with P i (u i ) = g(u i ) exp(iφ i (u i )), where g(u i ) is the input fundamental beam that is mildly focused to produce the quasi-widefield fundamental illumination beam and φ i (u i ) are pupil-plane aberrations imposed onto the illumination beam.
Estimation of these aberrations requires a set of THG fields, E (j) sig (x m ), that exhibit spatial overlap, and thus redundancy of the underlying χ (3) (x) that we wish to extract.These measurements are indexed by j in which the object is translated in the plane by x (j) s .

Transmission matrix
The transmission matrix, T (x o , x i ), is an operator used in linear optics S43 that provides a map from an input field, E i (x i ), with coordinates x i denoting the input plane where a field is incident on an optics system.The field exiting the system, E o (x o ), at an output plane with coordinates x o reads In a linear system, the transmission operator for a system with a planar scattering object, described by χ(x), is defined by The Green's function for illumination at the input plane, x i , and propagation to the scattering plane at x, is defined as H i (x, x i ).Similarly the exiting Green's function mapping the field from x to x o is denoted H o (x o , x).In the case of an isoplanatic (i.e., a shift invariant system), the transmission operator becomes Analogously to the case of second harmonic generation holography S31 , we define a transmission matrix for THG with an isoplanatic imaging system as Comparing this expression to Eq. S7, we see that they are identical if we make the following associations: ).Now, we may write the transmission operator as and where we note that we can exchange x i and 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 χ (3) (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.

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 3 1 (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 (x).Let us consider an idealized case where the fundamental field is a normally-incident, unity-amplitude plane wave, i.e., E 1 (x) = 1.Then, the THG transmission operator becomes We see that in this approximation, we obtain a low-pass filtered image of χ (3) (x).In the limiting case where H o (x) → δ(x), then we simply have However, because the illumination field has a finite spatial extent, the estimate of the THG field is obtained by integrating over the scan coordinate, x s , leading to Applying the definition of this estimator to Eq. S7, we obtain where we have defined the synthetic illumination function

D. 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.S43 We use the distortion operator approach inspired by Badon et.al to correct for pupil plane phase distortions.S29 This approach relies on the measurement of the distortion operator, D. Here, we show that our experimental setup directly measures D(x m , 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.S31 In the experiments, a set of THG field data are acquired, one for each specimen location, x s = x (j) s , as given in Eq.S5.This expression is the definition of the spatial-spatial THG distortion operator, D(x m , 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(u o , x s ).S29 To obtain this mixed-space operator, we transform the output spatial coordinate to the conjugate spatial frequency coordinate, u o , using the transform Using the definition of the spatial-spatial NDO in Eq.S5, we obtain Making the variable substitution x ′ = x − x s and using Identifying the quantity in the parentheses as the transmission operator, T (x o , x s ), we obtain the standard definition of the distortion operator S29,S31 where we have defined the spatial-frequency/spatial transmission operator as Applying our assumption of an isoplanatic system, we may write We may then write the distortion operator in the isoplanatic case as In order to understand how and when the aberrations can be estimated, we consider the correlation where and with the illumination correlation function defined as where we have identified ∆x = x ′ − x and ∆u = u o − u ′ o .The left singular vectors of the SVD of the distortion operator D will then yield the eigenvectors of the correlation DD * .
The form of Γ depends on two important characteristic lengths: the spatial support of the illumination, ℓ 1 , and a typical length for χ (3) (x), ℓ 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 ℓ 0 corresponds to the spatial support of χ (3) .Here, E 3 1 (x s ) is a function of x s /ℓ 1 and χ (3) (x) is a function of x/ℓ 0 .When ℓ 0 ≪ ℓ 1 , then the ∆x can be ignored in K and replaced by 0, so that and then 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 −1 o ≥ ℓ 1 (but still satisfying ℓ 0 ≪ ℓ 1 ), to extract the aberrations.In that case, DD * (u o , u ′ o ) is more extended around the diagonal and we have where we used the string of inequalities u In the opposite situation where ℓ 1 ≪ ℓ 0 , the analysis is similar as above with χ (3) and E 3 1 swapped.Indeed, x ′ is now localized around x in χ (3) * (x ′ ), which becomes χ (3) * (x), and changing variables as x ′ → x − x s , we find The first term on the right above localizes ∆u around 0 at the scale ℓ −1 0 , while the integral of K can be recast as Ê3 where Ê3 1 is the Fourier transform of E 3 1 .Based on the analysis for the case ℓ 1 ≫ ℓ 0 , the aberrations can then be reconstructed provided u −1 o ≥ ℓ 0 , namely the support of χ (3) is sufficiently small compared to u −1 o .The condition u m ≪ ℓ −1 1 is also necessary to estimate P o (u o ) for all u o since Ê3 1 (u o ) becomes small when the length of u o is greater than ℓ −1 1 .
When ℓ 0 corresponds to the scale of oscillation of χ (3) , we rewrite χ (3) as χ (3) (x) = ϕ(x/ℓ ϕ )e i2πuχ•x where ϕ is an envelope function and the length of u χ is ℓ −1 0 .Here ℓ 0 ≪ ℓ ϕ , and we assume ℓ ϕ ≫ ℓ 1 .In that case, a short calculation shows that Γ is as in (S31) with χ (3) replaced by ϕ and the integral of K becomes The correlation matrix is not too narrow when u −1 o ≥ ℓ ϕ .The term P o (u o ) has then to be separated from Ê3 1 (u o − u χ ), which is possible when u −1 o ≫ ℓ 1 , and P o (u o ) can be estimated for all u o provided ℓ −1 1 ≫ u m and the length of u χ , namely ℓ −1 0 , is less than u m .In the scenario where ℓ 0 and ℓ 1 are comparable, the expression of Γ does not simplify, but the kernel K still does localize u o around u ′ o at the scale ℓ −1 1 .Hence, as before, the aberrations can be reconstructed when the illumination is focused enough so that u −1 o ≥ ℓ 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(u o , x s ), which is the eigenvector of the correlation of the NDO, DD * (u o , u ′ o ).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.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 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 (∼ 4GB).
Our modified DASH algorithm is used here for correcting phase differences between the signal and reference field measured for each 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.

x s 2 ,
that described the power of the image, I c THG (x o ) = T (x o , x s ) d 2 estimated from a phase-corrected transmission matrix, T (x o , x s ) = T (x o , x s ) exp[−iφ r (x s )].Because aberration phase from the hologram phase difference induces a phase shift between otherwise identical complex vectors for each

FIG. 1 .
FIG. 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(x m , x s ), then input to our aberration correction algorithm.In the first step of the algorithm, the THG field images are shifted, mapping x m → x o + x m , which transforms the NDO into the transmission matrix, T (x o , x s ).THG field images synthesized with T (x o , x s ) using Ê3 (x m ) = T (x o , x s )d 2 x s from the raw data are highly distorted (no correction smiley face).The dominant relative phase aberration, φ r (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(u o , x s ), by estimating and correcting the pupil plane distortions φ o (u o ).The last step rescans the corrected spatial-spatial NDO, D(x m , x s ) to the coordinates of the corrected transmission matrix, T (x o , x s ).Spatial filtering and truncation of the singular value decomposition (SVD) of T (x o , 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.F{•} represents a fast Fourier transform operation that is being performed in the transformation to the D(u o , 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.

FIG. 2 .
FIG. 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 φ r (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 µm.

FIG. 3 .
FIG. 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 φ r (x s ).(Bottom) NDO pupil plane correction phase map, φ o (u o ).Scale bars are 10 µm.

FIG. 4 .
FIG. 4. MoS 2 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 µm.
(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.Supplemental Materials: Synthetic aperture holographic third harmonic generation microscopy A. 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 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).
FIG. S1.Experimental layout of the THG Holographic microscope for collection in a transmissiongeometry.The nonlinear crystal used for the reference is TiO 2 .The red line is the fundamental light and the purple line is THG.

s
's to build a synthetic spatial aperture, enabling extended FOV synthetic aperture microscopy.Sample scanning is done by translating the stage in steps of 2.78 µ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.
FIG. 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 (2NA/λ).The scale bar is 7µm.

)
This function gates the region of illumination of the sample.Under conditions of ideal convolution, the F 1 (x) → ½(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 3 1 (x) and H o (x o + x).With the estimated THG field, we then obtain the estimated intensity from I THG (x o ) = | Ê3 (x m )| 2 .

)
We need to introduce the scale of oscillation of the phase function φ o , denoted u o , to extract information from the equation above.We suppose then that φ o (u o ) is a function ofu o /u o .Since E 3 1 is localized at the scale ℓ 1 , K localizes u o around u ′ o at the scale ℓ −1 1 .Hence,for a nearly widefield imaging scenario when ℓ 1 is large compared to u−1 o , K(u ′ o − u o , 0) is narrow, and we can replace χ(3) (u ′ o ) * P * o (u ′ o ) by χ(3) (u o ) * P * o (u o ), resulting in DD * (u o , u ′ o ) = | χ(3) (u o )| 2 K(u ′ o − u o , 0),(S30)making DD * (u o , u ′ o ) nearly diagonal.

− 1 o
≥ ℓ 1 ≫ ℓ 0 and the change of variables u ′ o −u o → u ′ o .Since χ(3) (u o ) varies at the scale ℓ −1 1 , which is much greater than that of P o (u o ), the function χ(3) (u o ) is nearly constant compared to P o (u o ) and as a consequence P o (u o ) can be separated from | χ(3) (u o )| 2 and is approximately an eigenvector.Now, for P o (u o ) to be estimated for all u o , the largest wavenumber that can be measured by the experimental system, denoted u m , has to be much less than ℓ −1 0 since | χ(3) (u o )| 2 becomes small when the length of u o is greater than ℓ −1 0 .In that scenario, DD * is approximately of rank one and the leading eigenvector of DD * gives an estimation of P o (u o ).
FIG. 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.