1. 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].

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 [5–8]. 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 scattered power from each known focal point location. Because THG scattering is a coherent process, the total scattered power derives from the interference of THG light scattered through the focal volume.

In conventional laser-scanning THG imaging, no signal is observed in the far field when focusing in a homogeneous medium that is isotropic and exhibits normal optical dispersion. This lack of detectable THG occurs due to the axial symmetry of the Gouy phase shift that leads to destructive interference of the THG signal generated through the focus [9]. Any optical inhomogeneity on the order of the size of the beam disrupts this destructive interference, leading to detection of a strong THG signal [1,2]. Thus, the signal in THG imaging is normally generated by optical inhomogeneities that have a similar size to the beam focus [10].

THG microscopy can be used to study important cell and tissue properties by exploiting the physical origin of the THG signal strength that is driven by optical linear and nonlinear susceptibility heterogeneities and variation in refractive index (RI) mismatch that occur at the scale of the excitation focal volume. The biggest contributor to THG signals in biological samples is from lipids residing next to either water or cytoplasm. This lipid specificity, which has been validated by correlative imaging using fluorescent staining, leads to image contrast driven by cell membranes, lipid droplets, and lipid vesicles [10,11].

The ability of THG imaging to identify these features allows for anatomical identification of cell types across a variety of samples and myriad biomedical contexts. For example, THG has allowed imaging and identification of all cell layers in mouse retinas, including the identification of a new contrast source from pigment granules in retinal pigment epithelium cells [12]. THG can identify tissue and cellular morphological features considered classical hallmarks of cancer, such as cell irregularity and the ratio of nuclear to cytoplasm area in cells. This information is used to distinguish benign and cancerous breast tissue [13] and to identify infiltrative glioma in unstained, fresh human brain tissue [14]. 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 cancer subtypes [15]. Increased THG brightness and mean signal area from lipid droplets and membranes have been demonstrated as a unique marker used to identify HER2-positive breast cancer cells [16]. In addition, changes in THG signal intensity serve as a marker for microglial cell activation [17]. THG imaging of osteocytes and the complex porous network of bone can be probed through interface specificity revealed by changes of linear and nonlinear susceptibility, providing detailed anatomical images and informing our understanding of bone remodeling dynamics [18–20]. The size-dependent signal based on changes in the linear and nonlinear susceptibility also produces THG image signals for amyloid $\beta$ plaques, composed of densely packed amyloid $\beta$ proteins and lipids, and neurofibrillary tangles, composed of tau protein, allowing label-free differentiation of the plaques and tangles [21]. The same interface sensitivity mechanism for THG generation was used to observe both red and white blood cells, interfaces between cells and the extracellular matrix, peripheral nerves, perfused blood vessels, collagen bundles, adipocytes, striated muscle fibers, and microvesicles released from tumor cells [22], with application to monitoring invasion of metastatic melanoma cells into the mouse dermis with epi-collected THG signals [23].

It is worth noting that the brain and the central nervous system offer unique contrast for label-free morphological imaging due to a relative lack of THG contrast from cells. Here myelin sheaths are composed of multiple layers of lipid membranes, which produces a much stronger signal than a cell or nuclear membrane and thus dominates the THG signals from the brain. As a consequence, in brain imaging, cells appear as dark structures upon a bright background of THG signal [24–26]. Such images can be used to quantitatively assess various cell morphology biomarkers such as cell shape irregularity, density of cell nuclei, vascularization, and tumor-associated macrophages, allowing grading of gliomas [15].

Particularly as suitable laser systems become more readily available, laser-scanning THG imaging is increasingly playing an important role in biological imaging. However, laser scanning comes with several drawbacks. Traditional laser-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. Further, beyond these biological contexts, complex field THG microscopy offers new capabilities for imaging materials, including two-dimensional (2D) materials and metamaterials, as demonstrated in several studies [27–29] as well as in the present work.

In this paper, we introduce holographic imaging to directly record complex-valued THG fields, providing the first demonstration, to the best of our knowledge, of amplitude and phase-sensitive THG microscopy. This microscope is configured as a widefield, stage-scanning THG holography. Use of a lower local illumination intensity avoids photo damaging the sample [30–33]. 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. The configuration we describe enables large field of view (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 can in theory reduce the time needed for recording images. While this current work is performed with a sub-optimal laser source that cannot increase speed over laser-scanning THG microscopy, we demonstrate that our strategy in principle can improve speed, and we present the laser source parameters that can meet this requirement.

Finally, this work also leverages and develops adaptive optics approaches to correct THG field distortions arising from aberrations that lead to poor image quality of the illumination field. These aberrations originate from scattered light and wavefront phase distortions introduced by imperfections in the optics and inhomogeneities in the specimen. We have adapted a computational adaptive optics [34–38] 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. In conventional laser-scanning THG microscopy, previous work has demonstrated the use of active adaptive optics [39–41]. This work made use of measurements of a THG image, perturbations imparted mechanically in the optical line, and optimization of the perturbations to improve the THG image. Because we can directly measure the complex THG field, we can implement computational adaptive optic techniques for THG imaging [34]. We test our system on two pathology samples and a two-dimensional semiconductor, demonstrating its versatility for biological imaging and materials science applications. 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 [32] and to study the physical origin of complex-valued nonlinear susceptibility in both biological systems and materials [29].

2. 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 [30–32,38]. 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 [31] to boost weak THG signals.

 

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({{\textbf{x}}_m},{{\textbf{x}}_s})$, and then input to our aberration-correction algorithm. In the first step of the algorithm, the THG field images are shifted, mapping ${{\textbf{x}}_m} \to {{\textbf{x}}_o} + {{\textbf{x}}_m}$, which transforms the NDO into the transmission matrix, $T({{\textbf{x}}_o},{{\textbf{x}}_s})$. THG field images synthesized with $T({{\textbf{x}}_o},{{\textbf{x}}_s})$ using ${\hat E_3}({{\textbf{x}}_m}) = \int T({{\textbf{x}}_o},{{\textbf{x}}_s}){d^2}{{\textbf{x}}_s}$ from the raw data are highly distorted (no correction smiley face). The dominant relative phase aberration, ${\phi _r}({{\textbf{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({{\textbf{u}}_o},{{\textbf{x}}_s})$, by estimating and correcting the pupil plane distortions ${\phi _o}({{\textbf{u}}_o})$. The last step rescans the corrected spatial–spatial NDO, $\tilde D({{\textbf{x}}_m},{{\textbf{x}}_s})$, to the coordinates of the corrected transmission matrix, $\tilde T({{\textbf{x}}_o},{{\textbf{x}}_s})$. Spatial filtering and truncation of the singular value decomposition (SVD) of $\tilde T({{\textbf{x}}_o},{{\textbf{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. ${\cal F}\{\cdot \}$ represents a fast Fourier transform operation that is being performed in the transformation to the $D({{\textbf{u}}_o},{{\textbf{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.

Download Full Size | PDF

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}}}}({{\textbf{x}}_m})$ and a signal field ${E_{{\rm{sig}}}}({{\textbf{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 (${\sim}200\,\,{\rm{mW}}$ incident in the focusing optic) that illuminates the sample and a beam (${\sim}400\,\,{\rm{mW}}$ incident on the THG crystal) used to generate the third harmonic reference field. The fundamental field is gently focused (to a beam diameter of ${\sim}7\,\,\unicode{x00B5}{\rm m}$, approximately modeled as a Gaussian beam waist, which encompasses an illumination area roughly $15 \times$ larger than 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 {\text-} {\rm 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 ${\rm{Ti}}{{\rm{O}}_2}$ window, where strong THG light is generated [1,2,42]. The THG reference light is collimated and then combined with the THG signal field using a thin film polarizer positioned at $45{\rm deg}$. An additional half-wave plate is used to adjust the reference intensity to optimize hologram interference contrast, which is tuned for each specimen. 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}({{\textbf{x}}_m}) = |{E_{{\rm{ref}}}}({{\textbf{x}}_m}{)|^2} + |{E_{{\rm{sig}}}}({{\textbf{x}}_m}{)|^2} + E_{{\rm{ref}}}^*({{\textbf{x}}_m}) {E_{{\rm{sig}}}}({{\textbf{x}}_m}) + {E_{{\rm{ref}}}}({{\textbf{x}}_m}) E_{{\rm{sig}}}^*({{\textbf{x}}_m})$. We numerically isolate the interferometric pseudo-intensity term ${\tilde I_h} = E_{{\rm{ref}}}^*({{\textbf{x}}_m}) {E_{{\rm{sig}}}}({{\textbf{x}}_m})$ to estimate the complex THG image field using standard protocols [30,31]. 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 Supplement 1.

Each measured hologram at a sample position ${\textbf{x}}_s^{(j)}$ produces a 2D complex-valued THG field map, $E_{{\rm{sig}}}^{(j)}({{\textbf{x}}_m})$. An example set obtained by scanning ${{\textbf{x}}_s}$ is illustrated in Fig. 1. The spatial coordinates of the camera, defined by the pixel pitch, are mapped to measurement coordinates, ${{\textbf{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 ${{\textbf{x}}_s}$. Datasets of THG holograms are acquired by scanning the sample to both extend the spatial region that is imaged and 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_{{\rm{sig}}}^{(j)}({{\textbf{x}}_m})$, in the data stack is flattened into a linear vector that is inserted into a column of a matrix, $D({{\textbf{x}}_m},{{\textbf{x}}_s})$, as illustrated in Fig. 1. The column for each specimen scan position, ${\textbf{x}}_s^{(j)}$, is set to match the lexicographical order of the ${{\textbf{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) [38].

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 and 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 [38]; the averaged field extracted from the dominant singular vector constitutes the measured THG field, $E_{{\rm{sig}}}^{(j)}({{\textbf{x}}_m})$, for the $j$th sample position ${\textbf{x}}_s^{(j)}$ with a mean phase ${\phi _r}$.

Once we have obtained a set of THG field data, flattened the field, and placed them into an NDO matrix $D({{\textbf{x}}_m},{{\textbf{x}}_s})$, and then transformed the NDO matrix to the transmission matrix, the synthetic THG field is obtained using the estimator

The formula for this estimator given in Eq. (1) is derived in the supplements beginning with Eq. (S7), which is the model for the recorded THG image field for each ${\textbf{x}}_s^{(j)}$.

As noted, the measured data are intrinsically described by an isoplanatic model, where the object is translated through a stationary imaging system. Because the holographic measurement selects the polarization component of the THG field that is in the same direction as the reference field, we may describe the recorded field with a scalar formalism once we account for the reference field polarization and the incident fundamental field polarization. With this simplification the distortion operator describing the measured data in the space–space domain given in Eq. (S5) reads

The transmission operator is obtained by translating the measured THG fields by the shift imparted to the object, ${{\textbf{x}}_s}$, so that the fields are lined up in output space coordinates defined by ${{\textbf{x}}_o} = {{\textbf{x}}_m} + {{\textbf{x}}_s}$. The derivation of this mapping results in Eq. (S12) that defines the THG transmission operator, given by

Here, the coordinates in the shifted object space, ${{\textbf{x}}_s}$, are mapped to the space of the object, ${\textbf{x}}$, where THG scattering occurs, and ${\chi ^{(3)}}({\textbf{x}})$ is the third-order nonlinear optical susceptibility that provides the local amplitude and phase of the THG scattered light driven by the cube of the incident field [43]. The effective Green’s function for this mapping is given by the shifted cube of the incident fundamental field, $E_1^3({\textbf{x}})$, that serves as the driver of the local THG source. Mapping from the object source to the camera plane is facilitated by the Green’s function of the third harmonic 4-f imaging system, ${H_o}({{\textbf{x}}_m} + {\textbf{x}})$, in Eq. (3).

As shown in the supplements, application of the estimator, Eq. (1), to the THG transmission operator, Eq. (3), leads to the expression given by Eq. (S16) and reads

The estimator, Eq. (1), is used to produce a synthesized THG field, ${\hat E_3}({{\textbf{x}}_o})$. Then the image intensity is obtained from the synthesized field, i.e., ${\hat I_3}({{\textbf{x}}_o}) = |{\hat E_3}({{\textbf{x}}_o}{)|^2}$. For a discrete matrix, the integral is a sum of the discrete set of ${{\textbf{x}}_s}$. This sum is taken over complex-valued field components so that we can extract the amplitude and phase of the THG image field. For the synthesized field to be useful, the sum over ${{\textbf{x}}_s}$ must be largely devoid of distortions, particularly phase distortions. When successfully combined, the THG field provides a meaningful measurement of the spatial amplitude and phase variation of the measured THG field. However, to obtain such a field, an algorithm is needed to optimize the phase correction to produce an unaberrated spatial synthetic aperture THG image field.

3. THG HOLOGRAPHIC ABERRATION CORRECTION

Distortions to an image are bound to accumulate in any imaging process. In a complex-valued field imaging process as described here, phase distortions 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. Three distinct distortions to the final THG field images occur in our system. The first arises from the fact that any fluctuations in the relative phase between the signal and reference fields appear as relative phase distortions on the complex signal field extracted for each measured hologram. We denote this phase ${\phi _r}({{\textbf{x}}_s})$, and this is the first phase distortion that we correct with the algorithm. The second distortion that can arise in a system are specimen-induced aberrations. As the samples are mounted between a slide and a coverslip, the mounting can introduce aberrations. Furthermore, the tissue slices are ${\sim}6\,\,{\unicode{x00B5}}{\rm m}$ thick, and as the optical path represents ${\sim}6 \lambda$ and ${\sim}18 \lambda$ for the fundamental and harmonic fields, respectively, some limited specimen aberrations can accumulate in these samples. Finally, the ${{4 {\text -} f}}$ imaging system that magnifies the THG field on the camera for capturing holograms also introduces aberrations. These latter two aberration contributions are removed by the NDO. 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 handles the dominant phase distortion, ${\phi _r}({{\textbf{x}}_s})$, that arises from the hologram measurement process. This phase is random and uniformly distributed over $2\pi$, leading to speckle on an estimated image formed from uncorrected field data. The first step to correcting ${\phi _r}({{\textbf{x}}_s})$ maps the spatial–spatial NDO to the THG transmission matrix $T({{\textbf{x}}_o},{{\textbf{x}}_s})$. We built an optimization algorithm that maximizes the cost function, $J = \int I_{{\rm{THG}}}^c({{\textbf{x}}_o}) {d^2}{{\textbf{x}}_o}$, that describes the power of the image, $I_{{\rm{THG}}}^c({{\textbf{x}}_o}) = {| {\int \tilde T({{\textbf{x}}_o},{{\textbf{x}}_s}) {d^2}{{\textbf{x}}_s}} |^2}$, estimated from a phase-corrected transmission matrix, $\tilde T({{\textbf{x}}_o},{{\textbf{x}}_s}) = T({{\textbf{x}}_o},{{\textbf{x}}_s}) \exp [- i{\phi _r}({{\textbf{x}}_s})]$. The spatial overlap of adjacent holograms, as stated in the previous section, provides necessary redundancy in the data for our algorithm to converge. The regions of spatial overlap achieve maximum intensity when the adjacent fields are coherently combined in phase relative to one another. Thus, our algorithm maximizes the cost function to optimize the THG intensity estimated from the transmission matrix. This intensity optimization is one of many possible aberration-correction optimization metrics introduced by the astronomy community for adaptive optics [34,44,45]. These and other image quality metrics have been implemented by the microscopy community for linear complex field imaging [36,37,46–50] and fluorescent imaging [51,52] and form the basis for optimization in computational adaptive optics algorithms.

Because aberration phase from the hologram phase difference induces a phase shift between otherwise identical complex vectors for each ${{\textbf{x}}_o}$ position of the THG field, a correction of ${\phi _r}({{\textbf{x}}_s})$ maximizes the estimated object intensity at each spatial point by exploiting the fact that the spatial map is highly correlated along ${{\textbf{x}}_o}$. This is possible because the data shares the spatial structure of the sample nonlinear susceptibility, ${\chi ^{(3)}}({\textbf{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 [53]. Details of the optimization algorithm are provided in Supplement 1. With ${\phi _r}({{\textbf{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}({{\textbf{u}}_o})$, by using the concept of the distortion operator [36] that is adapted for nonlinear imaging [38]. To use the NDO for aberration correction, we map the “DASH corrected” spatial–spatial transmission matrix, $T({{\textbf{x}}_o},{{\textbf{x}}_s})$, into the spatial frequency–spatial distortion operator $D({{\textbf{u}}_o},{{\textbf{x}}_s})$. To do so, we first shift the spatial fields back from ${{\textbf{x}}_o}$ to ${{\textbf{x}}_s}$ to produce $D({{\textbf{x}}_m},{{\textbf{x}}_s})$. Then each 2D output field for each ${{\textbf{x}}_s}$ is Fourier transformed to the spatial frequency space, so that the NDO provides a relationship between ${{\textbf{x}}_s}$ and the output spatial frequency coordinates ${{\textbf{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({{\textbf{u}}_o},{{\textbf{x}}_s})$ [38]. Following a similar procedure for the linear DO [36], 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 ${{\textbf{u}}_o}$ across the ${{\textbf{x}}_s}$ spatial positions. Then, an inverse Fourier transform of ${{\textbf{u}}_o}$ brings back the spatial–spatial NDO, $\tilde D({{\textbf{x}}_m},{{\textbf{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({{\textbf{x}}_o},{{\textbf{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 ${{\textbf{x}}_s}$ location analogously to smart-OCT filtering [35] 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.

4. 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 [30]. Focal propagation distances are estimated using a gradient sharpness criterion [44].

It is interesting to consider the differences between the single brightfield images and THG amplitude and field images in Figs. 2–4 below. What is normally seen in conventional THG laser-scanning imaging is the total power collected within the NA of the objective. This includes interferences accumulated from propagation (phase matching effects) that is distorted by aberrations as the field propagates out of the sample. This conventional signal cannot be used to extract the phase of the THG field, and there is no way to correct for distortions of the signal imparted when the field propagates out of the sample and onto the detector. The brightfield image brings contrast from either absorption and scattering from the sample due to one-photon interactions with a UV LED. Given that the sample is thin, and we are measuring the intensity with the LED illumination, the contrast is dominated by linear absorption at the wavelength of the LED (365 nm). Moreover, the biological specimens are pathology slides stained with hematoxylin and eosin (H&E) dyes to highlight nucleic acids and proteins, respectively. We note that staining with hematoxylin has been shown to provide resonantly enhanced THG signal generation from cell nuclei, but not from the cytoplasm or the extracellular space or matrix [54–57]. Inspection of the images does indicate significant differences in the spatial amplitude information. More work will be required to understand the subtleties, but we can examine differences in these images.

 

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 ${\phi _r}({{\textbf{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\,\,\unicode{x00B5}{\rm m}$.

Download Full Size | PDF

 

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) First corrections applied in our algorithm to estimate and correct ${\phi _r}({{\textbf{x}}_s})$. (Bottom) NDO pupil plane correction phase map, ${\phi _o}({{\textbf{u}}_o})$. Scale bars are $10\,\,\unicode{x00B5}{\rm m}$.

Download Full Size | PDF

 

Fig. 4. ${\rm{Mo}}{{\rm{S}}_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\,\,\unicode{x00B5}{\rm m}$.

Download Full Size | PDF

The first set of images shown in Fig. 2 are of a $6\,\,\unicode{x00B5}{\rm m}$ thick epidural region of a mouse tail cross section prepared on a glass slide and with a cover slip (Triarch Incorporated). This is a prepared slide from Triarch with H&E staining and fixed with Bouin’s solution. THG reconstructed images are of the epidural region of the mouse tail cross section, with the LED ground truth image shown in Fig. 2(a). 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 Figs. 2(b) and 2(c). After first correcting the hologram phase distortion ${\phi _r}({{\textbf{x}}_s})$, the image quality improves significantly, as shown in Fig. 2(d). The fully corrected, aberration-free THG holographic images are shown in Figs. 2(e) and 2(f). Individual tissue strands in intensity image Fig. 2(e) 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. 2(f) can be identified in the phase map of the aberration-free image, reporting characteristics of sample morphology. In comparing Figs. 2(a) and 2(e), there are differences in morphology [dimmer regions in Fig. 2(a)] that should correlate with brighter regions in Fig. 2(e), but these are difficult to quantify without a detailed mathematical analysis. There are locations on the fibers where a $\pi$ phase jump is observed, possibly indicating a switch in collagen orientation; however, such speculation will require tensor analysis to be confirmed [58].

The next biological sample we demonstrate our THG holographic microscope on is a developing bone sample, approximately $6\,\,\unicode{x00B5}{\rm 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. 3(c) is of much higher quality compared to the intensity image shown in Fig. 3(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. 3(e), with the top panel’s phase map showing the phase corrections applied to the transmission matrix to correct hologram phase distortions, ${\phi _r}({{\textbf{x}}_s})$ and the bottom panel showing the pupil plane phase corrections acquired by the NDO step in our algorithm.

The final THG reconstructed image presented here is that of monolayer molybdenum disulfide (${\rm{Mo}}{{\rm{S}}_2}$) (Fig. 4). Monolayer ${\rm{Mo}}{{\rm{S}}_2}$ is a 2D semiconductor with a thickness of ${\lt}1\,\,{\rm{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 ${\rm{Mo}}{{\rm{S}}_2}$ monolayers were originally grown on a sapphire substrate via chemical vapor deposition, and then mechanically transferred to a glass microscope slide, after which a coverslip is affixed to keep the sample in place. 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 Figs. 2 and 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. Finally, there are intriguing features in the THG images that are not present in the brightfield LED image. In particular, there are at least two clear vortex phase structures in that phase with a corresponding null in the THG amplitude image. By contrast, the LED images are flat in these regions.

5. 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 laser-scanning microscopy except for a recent report on intensity widefield THG with a low-repetition rate amplifier [59]. 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)}}({\textbf{x}})$, and nonlinear THG, ${\chi ^{(3)}}({\textbf{x}})$, third-order 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. THG imaging signals for both conventional laser-scanning THG and widefield THG microscopes are determined by spatial heterogeneity of the spatially varying third-order nonlinear susceptibility, ${\chi ^{(3)}}({\textbf{x}})$. Heterogeneity of the third-order nonlinear optical susceptibility is defined by $\delta {\chi ^{(3)}}({\textbf{x}}) \equiv {\chi ^{(3)}}({\textbf{x}}) - {\langle {\chi ^{(3)}}({\textbf{x}})\rangle _{\textbf{x}}}$. Here we define the average, background, nonlinear susceptibility as a spatial average, ${\langle {\chi ^{(3)}}({\textbf{x}})\rangle _{\textbf{x}}}$. While both point-scanning and widefield THG imaging reveal the heterogeneity, $\delta {\chi ^{(3)}}({\textbf{x}})$, of the third-order susceptibility, the physical origin of the signal is based on distinctly different physics of propagation.

Despite the fact that linear susceptibility and third harmonic nonlinear susceptibility are both universally symmetry allowed, there are significant differences in the information revealed by the two modalities. In both cases, information of interest in imaging is carried by inhomogeneities in the optical susceptibility. In the case of linear susceptibility, the heterogeneity is defined analogously to the THG case as $\delta {\chi ^{(1)}}({\textbf{x}}) = {\chi ^{(1)}}({\textbf{x}}) - {\langle {\chi ^{(1)}}({\textbf{x}})\rangle _{\textbf{x}}}$. Again, the average linear susceptibility is defined as spatial average ${\langle {\chi ^{(1)}}({\textbf{x}})\rangle _{\textbf{x}}}$. We note that the refractive index is defined in terms of the linear susceptibility, ${n^2} = 1 + {\chi ^{(1)}}$, and a small refractive index perturbation at a boundary is expressed as $\delta n \approx \delta {\chi ^{(1)}}/(2 n)$. Below, we discuss the differences between conventional laser-scanning and holographic THG microscopy. We also discuss the new information brought by THG microscopy and phase-sensitive THG widefield imaging as presented here in particular. THG imaging provides additional information compared to linear imaging as the THG tensor, ${\chi ^{(3)}}({\textbf{x}})$, carries much more information in additional resonances and tensor elements than the linear tensor, ${\chi ^{(1)}}({\textbf{x}})$, which is a lower-rank tensor. Accessible information in the ${\chi^{(3)}}({\textbf{x}})$ tensor includes multiphoton resonances and organization of the molecular ultrastructure of macromolecules carrying harmonophores and spatial organization of materials.

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. We consider the RI at the fundamental and third harmonic wavelengths, ${n_1} = n({\omega _1})$ and ${n_3} = n({\omega _3})$. In normal materials, $\Delta n = {n_3} - {n_1} \gt 0$ due to the dispersion of the RI. The collinear phase mismatch parameter is defined in terms of the RI dispersion $\Delta k = {k_3} - 3 {k_1} = \Delta n {k_{3,0}}$. Here we have also written the THG wavenumber as ${k_3} = {n_3} {\omega _3}/c$ and ${k_{3,0}}$ as the THG wavenumber in free space. The effect of phase mismatch manifests differently in conventional laser-scanning THG and in widefield holographic THG described here.

The differences between THG imaging with laser 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)}}({\textbf{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 [9]. 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)}}({\textbf{x}})$ is broken by either an interface or a small inclusion along the direction of propagation of the beam [1,10]. The signal produced through disruption of THG destructive interference from throughout the focal volume is proportional to differences in optical properties at an inhomogeneity between two sample constituents, $a$ and $b$, which is expressed as ${| {\delta \alpha} |^2}$, where $\delta \alpha \approx [{\chi ^{(3)}}/{n_3} \Delta n] (\delta \chi _{\textit{ab}}^{(3)}/{\chi ^{(3)}} - \delta {n_3}/{n_3})$ [8], where $\delta \chi _{\textit{ab}}^{(3)} = \chi _b^{(3)} - \chi _a^{(3)}$ and $\delta {n_3} = {n_b}({\omega _3}) - {n_a}({\omega _3})$ are the nonlinear susceptibility and RI differences between the constituents. These differences are assumed to be small in this approximation, which is consistent with typical experimental conditions. The interference of the measured power on the detector is modulated by the variations in the RI within the focal volume of the nonlinear microscope. Images are formed in this configuration by recording some fraction of the scattered THG power, and such measurements lose all phase information. In short, laser-scanning THG imaging does not provide access to the phase of the THG field.

As seen above, conventional point-scanning THG microscopy is proportional to the square of the difference in fraction change in the third-order nonlinear susceptibility and the RI at the harmonic frequency. Thus, there is a dual dependence on the linear and nonlinear properties. By contrast, linear scattering-based imaging using laser-scanning microscopy can be implemented as either transmission or epi as a laser-scanning confocal microscopy. The signal model in transmission scales as $\propto {\rm{Im}}[\delta \chi _{\textit{ab}}^{(1)}]$, and only occurs for linear optical absorption and provides an image identical to brightfield microscopy. By contrast, reflectance confocal signals collected in the epi direction scale as ${| {\delta \chi _{\textit{ab}}^{(1)}} |^2} = ({\rm{Im}}[\delta \chi _{\textit{ab}}^{(1)}{])^2} + 4 {n^2} \delta {n^2}$, and thus depend on both absorption contrast and changes in the RI [60]. Here, we see, unlike the THG case, there is only a dependence on the changes in the RI of the sample. This combined with the rich detail about the sample that is encoded in the ${\chi ^{(3)}}$ tensor provides a much richer set of information that may be extracted with THG microscopy compared to linear scattering-based imaging.

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 the 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\,\,{\rm{mm}}$, where ${k_1} = {n_1} {\omega _1}/c$ is the wavenumber at the fundamental wavelength, and ${w_0}$ is the radius of a Gaussian beam focus model. Within the confocal parameter, the beam propagation can be reliably modeled as a plane wave. It is important to note that with planewave illumination, the brightness of the THG field will be moderated by phase mismatch, producing only very weak THG signal when $\Delta k L \gt \pi /2$, or for a thin sample when $\Delta k L \ll 1$.

For the case of planewave illumination in nonlinear scattering, very little axial information can be obtained in the forward-scattered geometry [30–32]. Nonlinear scattering is described as a homogeneous wave equation driven by a source term that does not depend on the harmonic field [32], and as such the model for nonlinear scattering is identical to that of linear scattering in the first Born approximation [60], but with a source term driven only by a power of the fundamental field [32]. In both linear scattering in the first Born approximation and nonlinear scattering with planewave illumination, the forward scattered field carries no information about axial structure in the object. Thus, in nonlinear scattering, even with weak focusing, the detailed structure of the axial structure of the object cannot be recovered; however, some three-dimensional information may be inferred by backpropagating the measured field when the fundamental field is gently focused into a sparse sample [30,31]. Three-dimensional spatial structure from nonlinear scattering can be obtained by spatial frequency diverse fundamental beam illumination coupled with a computational inversion of the 3D structure of the nonlinear susceptibility [32,33].

By contrast, backscattered fields are known to contain information about the axial structure of the object when collected over a broad spectral range. In linear scattering, such an experimental arrangement is generally implemented as optical coherence tomography [61,62]. In our recent work, we demonstrated the first backscattered coherent nonlinear scattered widefield optical sectioned measurement of a nonlinear field by resolving the SHG scattered in a synthetic aperture imaging holography system, where the depth and thickness of the recorded field were gated by the cross-coherence length of the holographic interference between short SHG signal and reference laser pulses [38]. By contrast, synthetic aperture transmission holography measurements in thin [38] and thick [63] samples do not provide axial information from the sample.

The THG widefield image is impacted by phase mismatch in a decidedly different way than with laser-scanning THG microscopy. Here we describe a nonlinear scattering problem where the amplitude of the scattered third harmonic field is much weaker than the fundamental field. As a result, we can invoke the undepleted pump approximation to describe the THG field scattering. This leads us to a model described by solution to the Helmholtz operator with a forcing term driven by the THG polarization density. The THG polarization density is proportional to ${\chi ^{(3)}}({\textbf{x}})$, and as the Helmholtz equation is linear and shift-invariant, we may consider the effects of the homogeneous contribution to the THG susceptibility, ${\langle {\chi ^{(3)}}({\textbf{x}})\rangle _{\textbf{x}}}$, separately from the heterogeneous fluctuations $\delta {\chi ^{(3)}}({\textbf{x}})$. We note that the spatial average term, ${\langle {\chi ^{(3)}}({\textbf{x}})\rangle _{\textbf{x}}}$, includes the average over the slide, and the coverslip and the heterogeneous component is the contribution from the object and thus contains the information that we seek in the THG image. The total solution to the THG field, ${E_{3t}} = {\bar E_3} + {E_3}$, is the sum of the THG field generated by the homogeneous driving term, ${\bar E_3}$, and the signal term given by ${E_3}$. Because the total length, $L$, of the sample (that includes the total thickness of the sample, slide, and coverslip) is long compared to the coherence length, ${\ell _c} = \pi /\Delta k$, the third harmonic field amplitude contribution for this part of the nonlinear susceptibility, ${\bar E_3} \propto {\rm{sinc}}(\Delta k L) \ll {E_3}$, can be neglected. As a result, the total solution is dominated by the THG scattering from the inhomogeneities, $\delta {\chi ^{(3)}}({\textbf{x}})$. This means that ${E_{3t}} \approx {E_3}$. It is this component that must the thin enough to avoid needing to solve the full 3D inverse scattering problem [32].

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}({\textbf{r}})$, that is imprinted onto the THG field. In addition, spatial inhomogeneities in ${n_1}({\textbf{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 susceptibilities. 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 [32].

The tensor structure of nonlinear susceptibility used to model THG generation contains rich detail about the specimen under investigation. Information in the THG signal can be encoded in the amplitude of the susceptibility that can be enhanced by one-, two-, and three-photon resonances [64,65]. In many samples, the harmonophores that generate the nonlinear signals for THG are isotopically 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 [66], as can structured illumination [67]. When the harmonophores are not isotopically oriented, such as in a birefringent crystal, THG images can be recorded with circularly polarized incident light [68].

In experimental imaging studies of the cornea, polarization-dependent THG scattering was used to select between specific tissue regions [69,70]. Isotropic tissue regions, as exemplified by the corneal epithelium, contain the usual interface specificity that reveals cell walls and nuclear membranes when illuminated with linearly polarized fundamental light, but negligible signal was produced with circularly polarized illumination. By contrast, the stroma consists of alternating thin layers of organized collagen, with the direction of orientation alternating between layers. The strong birefringence of this tissue produces weak THG light with linearly polarized excitation, but strong THG with circularly polarized excitation. Theoretical analysis indicated that the scattered THG field is elliptically polarized, meaning that the measured THG field can take on any phase value, and that the phase of the field could provide significant information about the organization of the stromal layers.

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. Third-order nonlinear susceptibility tensors contain many elements, and while many of the tensor elements vanish due to symmetry, the non-zero elements can be probed to reveal details not accessible to linear measurements, such as linear optical scattering. By analyzing the dependence of the imaged THG field on incident illumination polarization and the THG polarization, the relative amplitude of the ${\chi ^{(3)}}$ tensor components can be studied, which allows for the determination of the organization and orientation of molecular species generating the THG field [58,71–73]. While such tensor analysis is possible with SHG holography [74], the relaxed symmetry selection rules for THG permit application of this measurement strategy to a much wider range of structures in a specimen.

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 by 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 ${\rm{Mo}}{{\rm{S}}_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 [75,76], 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 [64,65]. We also observe phase vortices in Fig. 4, which suggests spatially resolving chiral THG scattering that has previously been inferred from temporal measurements that are necessarily spatially averaged [28]. THG field imaging will allow for spatially resolved characterization of nonlinear optical scattering from 2D materials [29] and metamaterials [27].

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, including optical damage and the laser source. 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 [77]. 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.

To realize this synthesis of a large FOV for THG field imaging, we have created an algorithm to both estimate and correct these random phase errors that inevitably occur in any experimental system. 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 [38]. 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.

To extract information from the NDO to correct aberrations, a separation of length scale of correlations of the object (nonlinear susceptibility) and the Green’s function represented in either the space or spatial frequency domain is needed. The NDO allows for correction of aberrations for an isoplanatic imaging system when in the spatial frequency-spatial space, $D({{\textbf{u}}_o},{{\textbf{x}}_s})$. Transformation from $D({{\textbf{x}}_m},{{\textbf{x}}_s})$ to $D({{\textbf{u}}_o},{{\textbf{x}}_s})$ is shown in Supplement 1, with the spatial frequency–spatial space [Eq. (S24)] reading

In Supplement 1, we define three correlation lengths: ${\ell _0}$, ${\ell _1}$, and ${u_o}$. ${\ell _0}$ quantifies the correlation length of the object (${\chi ^{(3)}}({\textbf{x}})$), ${\ell _1}$ quantifies the correlation length of the spatial extent of the spatial distribution of the THG scattering driving term, $E_1^3({\textbf{x}})$, and ${u_o}$ quantifies the correlation length of the output pupil.

For imaging, we are interested in resolving spatial structure in the object across each of the regions of illumination provided by ${E_1}({\textbf{x}})$. This means we are considering the scenario where ${\ell _0} \lt {\ell _1}$. To motivate the capability of the NDO for aberration correction, consider the case of a constant-valued susceptibility, where ${\ell _0} \ll {\ell _1}$. In this case, the approximate expression for the NDO reads

This is a rank one matrix, from which ${P_o}({{\textbf{u}}_o})$ and $E_1^3({{\textbf{x}}_s})$ may be estimated from the left and right singular vectors. Under general conditions, the information that may be extracted from the NDO can be analyzed by studying correlations of $D({{\textbf{u}}_o},{{\textbf{x}}_s})$ in ${{\textbf{u}}_o}$. In a detailed analysis in the supplements, we show that to extract correlations, we require that the fundamental field be moderately focused into the sample, with a beam size that still allows us to satisfy the condition ${\ell _1} \lt u_o^{- 1}$. This is the experimental scenario that we have implemented, and this allows for estimation and correction of aberrations in the THG holographic imaging system.

Analysis of the NDO demonstrates that the phase to be corrected must be smooth on a suitable length scale to be estimated with the SVD. In this analysis, we show that to be able to estimate any phase information, there must be a separation of scales. This requirement follows from the fact that the correlation lengths in the object susceptibility, illumination field, and phase correlation in the output pupil must differ suitably to provide structure in the NDO that carries information that can be extracted. This information is derived from eigenfunctions of the correlation matrix formed from the NDO. Such a requirement is not applicable to estimation of ${\phi _r}({{\textbf{x}}_s})$ because this phase is discontinuous from point to point in ${{\textbf{x}}_s}$, and as shown in Supplement 1 this phase is random and uniformly distributed in $2 \pi$. This is why we use an optimization algorithm that deploys an image-based metric to correct for ${\phi _r}({{\textbf{x}}_s})$. We tested several optimization algorithms and found that, in our implementations, DASH provided the best convergence to an optimum and found the solution more rapidly than other optimization algorithms. We tested several image metrics [44] to use in the cost functional for the optimization, and all the metrics that we tested exhibited similar performance.

Our present laser system is sub-optimal for improving imaging speed and reducing laser exposure. However, a few recent experiments have demonstrated the improved speed and reduced power exposure that is an advantage of widefield nonlinear microscopy [33,78,79]. We derive conditions for widefield nonlinear microscopy sources in a manner similar to those in the supplemental information of a paper describing a recent widefield coherent anti-Stokes Raman scattering experiment [33]. For illumination, we consider two different laser systems, with one for widefield illumination that has an illumination beam area of $A$ and a second laser system for point-scanning nonlinear microscopy, with an illumination beam area of $a$. Our calculations consider a harmonic order, $q$, of the nonlinearity for signal generation (here $q = 3$ is the THG). In this proof-of-principle experiment, each hologram exposure was set to 300 ms, with a set of 10 holograms recorded at each ${{\textbf{x}}_s}$ point. Reconstructing the set of 10 holograms, averaging the set using our SVD method to correct for fluctuating phase, and translation of the stage to each ${{\textbf{x}}_s}$ position takes approximately 7 s. Therefore, it takes approximately 10 s to reconstruct a single THG hologram at position ${{\textbf{x}}_s}$ and store it within the transmission matrix. As the entire FOV is constructed with 676 reconstructed THG holograms (${{26 \times 26}}$ holograms), the total data acquisition time is 6760 s. This imaging time could be significantly decreased using galvanometer mirrors to scan the illumination beam, an optimal laser source, and optimized signal averaging and beam overlap. The DASH correction for ${\phi _r}({{\textbf{x}}_s})$ takes approximately 10 min, and the NDO aberration correction for ${\phi _o}({{\textbf{u}}_o})$ followed by the spatial filter and SVD denoising takes approximately 5 min.

To account for the potential for faster imaging, where the total image time of the widefield image, $\Delta {t_A} = \eta \Delta {t_a}$, is shorter than the point-scanning image acquisition time, $\Delta {t_a}$, by the factor $\eta \lt 1$. As damage for sub-picosecond pulses is approximately independent of pulse fluence, we consider the case of similar pulse fluence, with ${J_A} = \alpha {J_a}$, where if $\alpha \lt 1$, we consider a lower fluence (and less optical damage and disruption) for the widefield case compared to the laser-scanning case. Here, ${J_A}$ is the pulse fluence of the widefield source, whereas ${J_a}$ is the fluence for the laser-scanning source.

Taking these two conditions and considering images of equal quality, where the same total signal is acquired in both cases, we obtain the conditions for the repetition rate, average power, and pulse energy relationship (assuming identical pulses and beam shapes). The repetition rate of the widefield laser is $\nu _r^{(A)} = {\alpha ^{- q}} {\eta ^{- 1}} (a/A) \nu _r^{(A)}$, where $\nu _r^{(A)}$ is the repetition rate of the laser-scanning source. This corresponds to an average power of the widefield source as ${P_A} = {\alpha ^{- q}} {\eta ^{- 1}} {P_a}$, and ${P_a}$ is the laser-scanning source average power. Finally, the relationship between pulse energies is ${{\cal E}_A} = (A/a) {{\cal E}_a}$. For typical parameters for laser-scanning THG microscopy, where ${P_a} = 20\,\,{\rm{mW}}$, $\nu _v^{(a)} = 40\,\,{\rm{MHz}}$, and ${{\cal E}_a} = 0.5\,\,{\rm{nJ}}$, the laser source for an imaging scenario where the pulse fluence is reduced by a factor of two ($\alpha = 0.5$) and the imaging speed is increased by a factor of 10 ($\eta = 0.1$), the widefield laser parameters would be ${P_A} = 1.6\,\,{\rm{W}}$, $\nu _v^{(A)} = 0.64\,\,{\rm{MHz}}$, and ${{\cal E}_A} = 5\,\,\unicode{x00B5}{\rm{J}}$. Our current laser system does not meet these parameters, but an optimal source would significantly improve imaging speed and reduce laser exposure.

While we demonstrate THG holographic imaging with thin samples in this paper, future work will exploit harmonic optical tomography (HOT) for full 3D imaging by solving the inverse scattering problem [32]. This will provide an estimate of the spatial distribution of ${\chi ^{(3)}}$ relative to the mean susceptibility value. In contrast, laser-scanning THG primarily provides a map of interfaces and inclusions, with a nonlinear relationship to ${\chi ^{(3)}}$. Both methods are useful, and the physics describing the signal are rather different. Some applications will be better suited to conventional THG microscopy and some to THG holographic microscopy. The optimal method will be driven by the application one is interested in studying.

6. CONCLUSION

The results here show the promising capabilities of our aberration-free, synthetic spatial aperture THG holography for 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 using a new computational adaptive optics algorithm. 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 laser-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., ${\rm{Im}}\{{\chi ^{(3)}}({\textbf{x}})\}$. The THG phase information from the atomically flat ${\rm{Mo}}{{\rm{S}}_2}$ sample provides new insight to material physics by depicting spatial phase variation and vortices that are not observed with intensity measurements. THG phase maps of biological tissues, as in the mouse tail sample, reveal interesting phase changes along individual collagen strands, which can be further investigated with polarization experiments to uncover structural details of organized tissue layers. A broad understanding of the implications and interpretation of the THG phase provides a new direction of research for nonlinear microscopy, and further experimental and theoretical work is required to fully unravel the information contained within the THG phase. In addition, in future work, three-dimensional image reconstructions can be produced by adapting the HOT image reconstruction algorithm [32] to aberration-corrected THG image field data, while the heterodyne amplification of the THG signal opens the path for widefield THG imaging deep inside of scattering media.