Effect of Scatterering on Coherent Anti-Stokes Raman Scattering (CARS) signals

We develop a computational framework to examine the factors responsible for scattering-induced distortions of coherent anti-Stokes Raman scattering (CARS) signals in turbid samples. We apply the Huygens-Fresnel Wave-based Electric Field Superposition (HF-WEFS) method combined with the radiating dipole approximation to compute the effects of scattering-induced distortions of focal excitation fields on the far-field CARS signal. We analyze the effect of spherical scatterers, placed in the vicinity of the focal volume, on the CARS signal emitted by different objects (2{\mu}m diameter solid sphere, 2{\mu}m diameter myelin cylinder and 2{\mu}m diameter myelin tube). We find that distortions in the CARS signals arise not only from attenuation of the focal field but also from scattering-induced changes in the spatial phase that modifies the angular distribution of the CARS emission. Our simulations further show that CARS signal attenuation can be minimized by using a high numerical aperture condenser. Moreover, unlike the CARS intensity image, CARS images formed by taking the ratio of CARS signals obtained using x- and y-polarized input fields is relatively insensitive to the effects of spherical scatterers. Our computational framework provide a mechanistic approach to characterizing scattering-induced distortions in coherent imaging of turbid media and may inspire bottom-up approaches for adaptive optical methods for image correction.


Introduction
Coherent anti-Stokes Raman scattering (CARS) microscopy is a nonlinear, label-free imaging technique that has matured into a reliable tool for visualizing lipids, proteins and other endogenous compounds in biological tissues and cells based on their spatially-dependent third-order polarization [1][2][3]. In the CARS process, a pair of incoming beams (named "pump" and "Stokes") are exploited to coherently and resonantly excite selected vibrational levels of a population of molecules. To this end, the beam frequencies are chosen so that their difference matches a vibrational frequency of the oscillating dipoles of interest. As a consequence of the interaction of the vibrationally excited molecule with a third photon, the nonlinear polarization radiates through emission of a fourth photon: the CARS signal [4].
CARS microscopy is most commonly executed in the point illumination mode, in which the signal is generated in the focal volume of a high numerical aperture lens. Similar to all forms of microscopy that rely on the formation of a tight focal spot, the CARS signal is sensitive to the characteristics of the three-dimensional focal volume. Distorted or aberrated focal fields generally compromise CARS signal generation and degrade the CARS signal [1,3]. In contrast to other nonlinear optical microscopy techniques, such as two-photon excited fluorescence (TPEF) [5], CARS microscopy relies on the spatial phase of the excitation field and is particularly sensitive to wavefront distortions. As a result, the CARS emission is dictated by both the amplitude and the phase of the focal fields. Moreover, since the pump beam and the Stokes beam have different wavelengths, their focal fields may exhibit different aberration characteristics.
The heterogeneity of biological samples, which results from structures of variable size and effective refractive index, modifies the propagation of focused optical wavefronts resulting in distorted focal volumes [6]. The scattering-induced modification of the focal volume distribution is the primary factor for the deterioration of CARS signals at greater sample depths in turbid samples and results in attenuated signals, reduced contrast, and degraded resolution [7]. While these effects may be less pronounced in thin samples such as cell cultures, refractive index variations still affect the focal volume and can alter the CARS radiation profiles, leading to signal loss or unaccounted image artifacts.
The deleterious effects of light scattering in coherent imaging methods can be mitigated by shaping the excitation optical wavefront to compensate for the anticipated scattering-induced wavefront distortions [8][9][10]. Such adaptive optics approaches offer the possibility to restore signal levels and retrieve high resolution images in turbid media [6]. In the context of linear optical microscopy, wavefront shaping techniques have been used to almost completely counteract the effects of light scattering, or to leverage scattering in the medium to achieve image resolution surpassing that obtained in non-scattering samples [6]. In recently published work, Judkewitz and co-workers [11] used scattered electric field point spread function as a guidance to compensate the effect of scattering. Such adaptive optics method may not work when transmission signals acquired from reference and scattered beams lack sufficient correlation. Adaptive optics approaches have also been applied to CARS microscopy, by using the maximization of the CARS intensity as an objective function to optimize the shaping of the excitation wavefront [12].
Virtually all adaptive optics approaches are based on empirical optimization of experimentally accessible parameters, such as the signal intensity [8][9][10][11][12]. In this approach, the sample is considered a black box, which can be characterized by an effective transmission matrix that does not require a detailed understanding of the physical origin of the wavefront distortions. In many cases, such a strategy has proven to work well for counteracting scattering effects in linear optical microscopy applications. However, in nonlinear optical microscopy, there is evidence that maximizing signal intensity may not represent an appropriate optimization metric, resulting in the convergence to local extrema that correspond to focal shapes and positions that are markedly different from those obtained under non-scattering conditions [13]. This possibility underscores that a general strategy to manage the deleterious effect of light scattering effects must go beyond empirical optimization of signal intensities. This notion is particularly pertinent to CARS microscopy, where subtle amplitude and phase effects can have dramatic consequences for the observed signal intensities [14]. Instead of tackling the problem through an empirical black box approach, a fundamental understanding of the physics that gives rise to scattering artifacts in CARS is imperative. In this regard, a bottom-up, computational approach, that considers how wavefront aberrations affect CARS imaging, may provide the insights necessary to devise experimental approaches for recording CARS images devoid of scattering artifacts.
Such a detailed, fundamental understanding of linear scattering effects in coherent Raman scattering does not currently exist. Several model-based approaches have been used to investigate the effect of light scattering on the generation of coherent Raman signals in scattering media [15,16]. These include the use of Monte Carlo methods to simulate Raman scattering in turbid samples [16]. However, Monte Carlo simulations are unable to rigorously model diffraction or properly account for the amplitude and phase characteristics of propagating fields. These deficiencies prevent Monte Carlo simulations from accurately modeling spatial coherence, which is a critical determinant for the generation of coherent nonlinear optical signals. While full-field simulations can be conducted using finite-difference time domain (FDTD) methods to study the effect of scatterer size and orientation on near-field CARS signals [17][18][19], they are impractical for extensive parametric studies due to the substantial computational cost [20].
In this work, we aim to take several important steps toward building a fundamental, real-space picture of how linear scattering affects experimental observables in CARS microscopy. Recently we introduced a new efficient method to compute focal field distortions produced by scattering particles using Huygens-Fresnel wavelet propagation [21] and field superposition methods [22]. This Huygens-Fresnel Wave-based Electric-Field Superposition (HF-WEFS) approach provides accurate focal field predictions in the presence of single or multiple scatterers with arbitrary size, spatial configuration, density and orientation. Here, we apply a computational framework that employs HF-WEFS to examine CARS signal generation and far-field detection in the presence of scattering. Our framework first employs the HF-WEFS method to compute scattering-induced focal volume distortions of both the pump and Stokes beams. Next, we determine the CARS signal generation by computing the spatially dependent third-order dielectric polarization density produced by the pump and Stokes fields. Finally, we use the radiating dipole approximation [1,23] to compute the CARS signal as measured by a far-field detector. This approach enables the simulation of the far-field CARS signal with pump and Stokes beams of arbitrary polarization state, spatial distribution, illumination and detection numerical aperture, scatterer configuration, and scatterer shape.
As test samples, we simulate an isotropic solid sphere, a myelin cylinder and a myelin tubular structure. Myelin is a biological structure that envelopes a subgroup of nerve fibers in the gnathostomata and functions to increase nerve impulse conduction efficiency. We chose to simulate myelin due to its biological relevance, morphology and molecular characteristics that make it suitable for CARS imaging. CARS microscopy is frequently employed in myelin imaging, thanks to the strong CARS signal obtained by targeting its extremely abundant CH 2 bonds. Consequently, myelin has been studied using CARS imaging under normal physiologic [24,25] and pathologic [26][27][28][29] conditions.

Methods
Our framework deconstructs the process of CARS excitation, emission and detection into three sequential computations: (a) focused beam propagation in a scattering medium, (b) production of a nonlinear polarization field within the focal volume, and (c) far-field dipole radiation. A schematic of these components is shown in Fig. 1. We use the HF-WEFS method to rigorously model the focal fields generated by the propagation of pump and Stokes beams.
Once the focal fields have been computed, we compute the third-order dielectric polarization density, P (3) (r) produced by the incident pump and Stokes electric field distributions within the focal volume based on the nonlinear susceptibility of the medium. The emission that follows from P (3) (r) is then modeled as a collection of radiating dipoles in focus, which couple to, and is detected in, the far-field [1,23].
We detail each of the processes represented in Fig. 1 in the following subsections.

Focus Beam Propagation
We consider monochromatic pump and Stokes beams incident upon an aplanatic lens, and propagating independently towards the focal volume. In this study, we use the fundamental Hermite-Gaussian spatial mode (HG00) for both pump and Stokes beams. The electric field amplitude distribution of a Gaussian beam at the plane of an aplanatic lens can be expressed as [30]: where E 0 = 1 and ω 0 is the radius of the Gaussian beam at which the electric field amplitude falls to 1/e of the maximum axial value. The aplanatic lens system can be geometrically represented by a reference spherical surface that has a center at the origin [22,30]. The HF-WEFS method considers forward propagation of Huygens-Fresnel spherical waves from the reference spherical surface. We determine the propagation origin of each HF spherical wave at the lens surface by generating a set of uniformly distributed points on the reference surface [31]. Each spherical wave is represented by the summation of outward propagating Huygens-Fresnel plane wavelets (HF wavelets) [22,32]. In the absence of linear scattering in the space between the lens surface and the focal region, this Huygens-Fresnel description accurately reproduces the three-dimensional, diffracted-limited focal volume as predicted by diffraction theory [32]. The amplitude of an HF wavelet at each radiating point is given by |E inc (x, y)|. The parallel and perpendicular electric field components (E , E ⊥ ) of the HF wavelet at the spherical reference surface are given by [22]: where JV is the Jones vector that describes the polarization of light, and n inc and n are the refractive indices of the medium before and after the lens. φ and θ are azimuthal and polar angles of the HF plane wavelet with respect to the global coordinate system. The unscattered electric field components (E unscat , E unscat ⊥ ) at a distance d from the point of emission can be expressed as: where k is the wave number = 2π/λ. When scatterers are present in the medium, we consider each scatterer sequentially and account for all possible HF plane wavelets that may interact with it. In this study, we select spherical scattering particles, for which full-amplitude scattering matrices can be readily obtained using Lorenz-Mie theory [33]. For a scatterer located at point D, the parallel and perpendicular polarization components of the scattered electric field for a specific polar angle θ s and distance from the scatterer r s can be expressed as [22]: where E D and E ⊥D are the parallel and perpendicular incident electric field components at point D. The unscatterted and scattered fields can be superposed to obtain the total field at any location [33]. The parallel and perpendicular electric field components calculated in Eqs. 3 and 4 are transformed into x, y, and z components before superposition [22]. The components of the total electric field at a location r, E(r), can be computed as: Equations 1-5 are used to propagate the pump beam E p (r) and Stokes beam E S (r) in a scattering medium to obtain their x, y, and z components of the electric field in the focal volume.

Polarization Signal Computation
In CARS microscopy, the i th component of the spatially-dependent third-order dielectric polarization density induced at location r by the pump electric field and Stokes electric field is computed from: where i = (x, y, z) and χ (3) i jkl (r) is the third-order non-linear susceptibility tensor of the objects or media. E pj (r) and E pk (r) are electric field components of the pump beam and E * Sl (r) is conjugate electric field components of the Stokes beam at location r.
In this study, we consider the third-order nonlinear susceptibility tensor of spherical objects and cylindrical and tubular myelin structures placed within the focal volume. The nonlinear susceptibility is a tensor of rank 4, with 81 elements. The number of nonzero and independent elements depends on the spatial symmetry of the sample object. We assume the spherical objects to be uniform and isotropic, which results in 21 nonzero tensor elements, of which only four are independent ( [34]. For the cylindrical and tubular myelin structures, we employed published tensor element values that were experimentally determined for myelin sheaths [35]. Although different from the isotropic case, the nonlinear susceptibility of myelin sheaths is also described by 21 nonzero elements, with four independent tensor elements [35][36][37][38]. Because myelin layers are organized in concentric cylinders, their constituent molecules are rotated with respect to the laboratory frame depending on the location in the myelin structure. To model the measured response in the laboratory frame, the molecular nonlinear susceptibility is rotated with the proper Euler angles to find the overall CARS response of the system [38].

Far-field Dipole Radiation
Once the nonlinear polarization is determined within the focal volume, the resulting far-field CARS emission can be modeled using an ensemble of radiating dipoles [1,23]. For this purpose, each volume element in the vicinity of the focus is considered a point dipole. The magnitude of the dipole is given by Eq. 6. Each dipole radiates, and the resulting electric field is detected in the far field. The total amplitude of the electric field, E C (R), at a far field location R is the sum of the amplitude contributions from all point dipoles emanating from r [30,39]: where k C = 2π/λ C , λ C is the CARS wavelength in the medium, and V is the excitation volume.
To calculate angular resolved CARS radiation patterns shown in Figs. 3 and 4, we compute |E C (R)| 2 by making use of Eq. (7). The total CARS signal intensity I C captured by the far-field lens with an acceptance angle of α max can be written as [14] I C ∝ To obtain the CARS intensity as a function of the y-z grid (Fig. 5) and to simulate CARS images (Figs. 6 and 7), we compute the total CARS intensity using Eq. 8.

Numerical Simulation
In this study, the wavelengths of pump and Stokes beams are selected as λ = 800 nm and 1064 nm, respectively. We consider HG00 beams with filling factors ( = ω 0 / f NA ex ) equal to unity [30], where f is the focal length of the lens. We consider (n/n inc ) = 1 and compute the excitation within a 3 µm × 3 µm × 6 µm volume centered about the focal point. This volume is subdivided into a three-dimensional grid with 50 nm cubic voxels. We compute the distorted pump and Stokes electric fields at each grid point separately using Eqs. 1-5.
We consider the CARS imaging of three separate objects: 2 µm diameter sphere, 2 µm diameter myelin cylinder and 2 µm diameter myelin tube. The myelin tube has wall thickness of 250 nm and is centered or offset from the optical axis. The refractive indices of the medium and the scatterers are 1.33 and 1.49, respectively. Even while the CARS active objects have different refractive indices, we assume them to be index matched when modeling light propagation. The χ (3) of each object is considered as non-resonant, i.e., we ignore tentative phase effects due to the presence of spectral resonances. The values of the nonlinear susceptibility tensor elements of the objects are obtained as described above. The χ (3) i jkl of the surrounding medium, including the empty center portion of the tube, is set to zero. The x, y, z components of P (3) (r) are computed using Eq. 6, with χ (3) i jkl and the electric field distribution of pump and Stokes beams as inputs. After calculating P (3) (r) in the volume element of each grid point, the far-field amplitude is computed using Eq. 7. Computation of the CARS far-field emission is accomplished by placing a hemispherical detector in the far-field. The total CARS intensity is computed by integrating the far-field CARS radiation pattern over the detector acceptance angle, as in Eq. 8. We consider detection with acceptance angles of 71.8 • (NA det = 0.95) and 33.4 • (NA det = 0.55). We examine the CARS signal under scattering and non-scattering conditions for different excitation numerical aperture (NA ex = 0.825 and 0.55). Figure 2 depicts the various simulation geometries. In Fig. 2(a) we depict the effect of scatterer locations within the y-z grid on the far-field CARS intensity. The y-z grid has an overall dimension of (y, z) = 10 µm×8 µm with 0.5 µm spacing. In Fig. 2(b) we depict the generation of CARS images using point illumination without scattering as references. In Figs. 2(c) and (d), we depict two cases used to examine the effects of a discrete scatterer on CARS imaging. Figure 2(c) considers the effect of a 2 µm diameter scatterer placed along the optical axis 5 µm below the focal plane. Figure 2(d) considers the same scatterer placed at the same depth but offset 1.5 µm to the right of the optical axis. We consider CARS images generated using x-polarized and y-polarized light for both pump and Stokes beams separately.

CARS radiation profiles with a scatterer
We first consider the CARS radiation profiles resulting from pump and Stokes beams in the absence and presence of scattering objects. Figure 3(a) provides far-field CARS radiation patterns for different objects located at the focus in a medium without scattering. In these computations, the excitation fields are focused by a microscope objective with NA ex = 0.825 and the far-field radiation is detected using a lens with NA det = 0.95. The CARS emission intensity is shown as a function of the collection angle in the far-field. For comparison, each radiation profile is multiplied by a factor shown in brackets to provide plots have the same maximum radiance relative to the non-scattering case of the 2 µm spherical object. The inset of each panel shows the amplitude and phase cross sections (y-z plane, x=0) of P (3) (r) in the focal volume. Displays of the phase cross-section are masked with the amplitude distribution to emphasize the regions of the focal volume that contribute most significantly to the CARS emission. Insets to the left of each radiation pattern show y-z cross-sections of the amplitude (upper) and phase (lower) of P (3) (r). Insets in rows (c) and (f) show the amplitude (left) and phase (right) differences of (b) and (e) relative to the corresponding non-scattering cases, (a) and (d), respectively. Each inset spans 2µm × 2µm. Each radiation profile was multiplied by the number in the bracket to provide same maximum radiance. The percentages in (b) and (e) indicate the CARS intensity relative to the corresponding non-scattering case. Detection numerical aperture is fixed at NA det = 0.95.
In the non-scattering case, the CARS emission is highly forward directed and results from the phase-matching of the CARS radiation along the optical axis [14]. This situation changes when a scatterer is introduced in the vicinity of the focal excitation volume. In Fig. 3(b), we show CARS radiation profiles for these same objects in cases where a scatterer is positioned along the optical axis 5 µm below the focal plane. The insets show profiles of the amplitude and phase differences of P (3) (r) in the focal volume relative to the non-scattering case. The wavefront aberrations produced by the scattering object result in a nominal shift of the maximum amplitude of the polarization density to positions just below the focal plane. The scattering also distorts the phase profile of the induced polarization. Along the optical axis, scattering of the excitation fields introduces an extra phase shift in the nonlinear polarization approaching π/3 across the focal volume. This additional phase shift is responsible for the reduced intensity and modified angular distribution of the CARS radiation profiles. When the object is centered about the optical axis, the forward directed CARS signal is depleted significantly, whereas the off-axis radiation is more prominent. This is also observed in the shifted myelin tube (fourth column of Fig. 3(b)) where uneven amplitude and phase profiles within the tube contribute to an asymmetric CARS emission profile. These results clearly illustrate that the presence of a scattering particle not only modifies the overall amplitude of the nonlinear polarization, but also the spatial phase distribution. The lobed radiation pattern results from a scattering-induced phase shift along the optical axis, as can be inferred from the spatial phase profiles of the nonlinear polarization. This observation highlights the need to consider both the amplitude and phase of the excitation fields to properly account for the interference effects that occur within the focal volume.
In Fig. 3(c) we show how these CARS radiation patterns are altered when using excitation illumination with a reduced numerical aperture (NA ex = 0.55) in a non-scattering medium, while keeping the detection NA unchanged (NA det = 0.95). The smaller illumination NA ex introduces a narrower range of spatial frequencies into the sample and results in a broader and more elongated focal excitation volume. The longer interaction volume provides a more directional, phase-matched CARS signal along the z-axis. Figure 3(d) displays the CARS radiation profiles in the presence of the scatterer. In contrast to the results of Fig. 3(b), when using the smaller illumination NA the CARS signals from the solid sphere and cylinder remain highly directional along the optical axis with much smaller changes in the radiation pattern. This shows that the CARS emission of solid objects are less sensitive to phase aberrations carried by the smaller spatial frequencies associated with the lower NA of illumination. By contrast, CARS emission from the hollow myelin tubes remain sensitive to scattering-induced phase changes carried by the lower spatial frequencies of the excitation light, resulting in more CARS radiation profiles that remain distorted. Also noteworthy is the lack of attenuation of the CARS signal in the case when the edge of the myelin tube is at focus, as compared to the centered case.
The percentage values shown in Fig. 3(b) and Fig. 3(d) provide the total CARS signal intensity relative to the non-scattering case for each NA ex . As expected, the presence of the scatterer attenuates the CARS intensity for all objects. The attenuation relative to the non-scattering case is larger when illuminating the sample with the lower NA. However, the angular distribution of the CARS radiation remains very directional.
In Fig. 4, we show the effect of lateral particle position on the CARS radiation patterns for single scatterers placed at different positions along the y-axis, 5 µm below the focal plane (z = −5 µm). We consider scatterer locations of y = −2, −1, 0, 1 and 2 µm. The scatterer diameter is varied from 1 to 4 µm (Figs. 4(a)-(d)). For a 1µm scatterer diameter (Fig. 4(a)), the distortion and attenuation of the far-field radiation profile is minimal. Larger scattering particles provide more substantial amplitude attenuation and phase distortion resulting in more pronounced variations in the CARS radiation profiles and overall signal attenuation. The largest attenuation and distorted radiation profiles are seen for the 4 µm diameter scatterer (Fig. 4(d)) because the scattering induced phase shift in the nonlinear polarization along the optical axis approaches 2π (Fig. 4(d)). The peak intensity along the optical axis is greatly affected for scatterer locations directly under the spherical object. Importantly, highly asymmetric radiation profiles are produced when the Fig. 4. The far-field CARS radiation patterns from a 2 µm diameter solid sphere (blue) located at the focal point in a medium with a single scatterer (gold) placed 5µm below the optical plane at y locations of y = −2, −1, 0, 1, 2 µm as shown. The effect of scatterer size is shown for diameters of (a) 1 µm, (b) 2 µm, (c) 3 µm, and (d) 4 µm. Each radiation profile was multiplied by the number in the bracket to provide same maximum radiance. Insets to the bottom of each radiation pattern show y-z cross-sections of amplitude difference (left) and phase difference (right). Amplitude/phase differences are calculated by subtracting amplitude/phase of P (3) (r) induced in a non scattering medium. Excitation and detection numerical apertures are fixed at NA ex = 0.825 and NA det = 0.95. particle is displaced laterally from the optical axis. These asymmetric profiles result from spatial phase distortions carried by the nonlinear polarization in the focal volume and directly impact the CARS signal detection. We provide computed amplitude differences and phase differences of the nonlinear polarization relative to the non-scattering case. Figure 4 demonstrates the impact of lateral scatterer position on the angular profile of the CARS emission.

Effect of scatterer position and detection numerical aperture on total CARS intensity
We now examine how the y-z position of a single scattering particle impacts the total CARS intensity with NA det = 0.95. Figure 5 shows the integrated CARS signal from several objects as a function of position of a 2 µm diameter spherical scatterer. The scatterer positions correspond to those in the y-z grid shown in Fig. 2(a). We start from z = −4 µm to avoid overlap between the scattering particle and the focal volume under consideration. Figure 5(a) shows that for the case of imaging a 2 µm diameter solid sphere, myelin cylinder and myelin tube, the attenuation of the CARS signal is more prominent for scatterer positions more proximal to the focal volume. Fig. 5. The far-field CARS intensity as a function of the y-z particle location grid. The object that is placed at the focus is 2 µm sphere, 2 µm myelin cylinder, and 2 µm myelin tube (centered and shifted by 0.875 µm left of the optical axis). Scatterer diameter is 2 µm. NA det is (a) 0.95 and (b) 0.55. (c) Intensity ratio after dividing (b) by (a). The center value (white color) of the ratio color bar represents the ratio obtained for the non scattering medium. Excitation numerical aperture is fixed at NA ex = 0.825. Figure 5(b) provides these same results for a reduced detection numerical aperture of NA det = 0.55. We observe a similar trend for the variation in the total CARS intensity although the overall CARS intensity is reduced due to the lower collection angle. The effect of the detection NA is emphasized in Fig. 5(c), which displays the ratio of the CARS intensity obtained with NA det = 0.55 divided by that obtained using NA det = 0.95. The color code in the ratio images has been scaled relative to the non-scattering case, where the white color corresponds to the ratio obtained when no scatterer is present. Blue regions indicate positions where the CARS signal ratio using these two detection NA's is smaller as compared to the non-scattering case. For majority of scatterer locations, the relative amount of signal loss due to scattering is greater for the lower detection of 0.55. Omission of large angles that have more signal contribution as shown in Fig. 4(b)) increases the relative signal loss. Red areas, however, correspond to scatterer positions where higher ratios are observed. These latter regions tend to be in the shadow of the geometrical focus. The higher ratios occur due to scattering-induced redirection of excitation field density to areas that are otherwise depleted of excitation energy.

Effect of incident polarization on CARS imaging in scattering samples
In the previous Sections, we examined the variation of the CARS emission profiles and total signal as a function of scatterer position. These results show that both the angular distribution and the intensity of the CARS signal depend on the scatterer size and location. Here, we examine the effects of scatterers on CARS images, and how these images are affected by the polarization state of the excitation beams. We first consider CARS images in the absence of scattering particles. In Fig. 6 we show simulated images of three objects: the 2 µm diameter solid sphere, myelin cylinder and myelin tube considered previously. Figure 6(a) shows simulated CARS images obtained when both input beams are x-polarized, whereas Fig. 6(b) provides images obtained using y-polarized incident beams. The differences in the images obtained using these different polarizations is a direct consequence of the anisotropy of χ (3) . Figure 6(c) displays the ratio of the xand y-polarization images. The CARS ratio image of the sphere is uniform because an isotropic χ (3) was chosen. We also see a ratio of 1 in the middle portions of the myelin cylinder. The outer boundaries of the myelin cylinder and myelin tube have polarization ratios larger than 1. Recall that the cylinder and myelin tube are modeled as radially-ordered lipid membrane sheets, which exhibit a highly anisotropic χ (3) in the laboratory frame that changes with the orientation of the lipid. The non-uniform ratio images thus reflect the variation of the lipid orientation.  Figure 7 shows simulated CARS images of the same objects in the presence of a scatterer at different locations. The scattering scenarios considered here correspond to those shown in Figs. 2(c) and 2(d). Figs. 7(a) and 7(d) provide images obtained using x-polarized incident beams and Figs. 7(b) and 7(e) provide images obtained using y-polarized incident beams. Fig. 7(c) and 7(f) provide the polarization ratio images. In both scattering scenarios, the images are significantly distorted by the presence of the scatterer. For scatterers positioned along the optical axis, the center of the image is attenuated as a direct consequence of the scatterer position as seen in Fig. 4. This is clearly seen for the myelin cylinder and myelin tube, whereas the scattering induced uniform amplitude attenuation is observed for the sphere. Once the scattering particle is offset by 1.5 µm relative to the optical axis, the image is distorted in an asymmetric pattern. Comparison of the two scenarios clearly demonstrate that scattering-induced distortions in CARS images are strongly affected by scatterer position. Moreover, the x-polarized and y-polarized images are affected somewhat differently by the scatterer. This difference originates primarily from the anisotropy of the lipid sample, which was also evident in Fig. 6. As a result, the polarization ratio images do show sensitivity to the presence of the scatterer. The polarization ratio images shown in Figs. 7(c) and 7(f) differ by less than 5% from the ratio images simulated in an non-scattering medium (Fig. 6(c)). This is a consequence of the invariance of the scattered field from spherical particles relative to the linear polarization direction of the excitation field. Thus, while the CARS intensity image is affected significantly by the presence of scatterers, the polarization ratio image, which displays the anisotropy of the dipolar Raman scatterers, provides a much less distorted view of scattering samples.
Because of the complexity of the problem, we have avoided linear scattering by using an index-matched object in the focal volume. It is known that refractive index mismatches between objects and the surrounding medium within the focal volume give rise to additional effects [40]. In particular, the internal and external scattered fields from refractive index mismatched objects slightly alter the electric fields in the P (3) (r) calculation. In addition, the phase and amplitude differences between vibrationally resonant particles and the nonresonant medium can modify the CARS radiation in the far-field. [1,19,41,42] However, including such additional effects complicates the analysis of wavefront distortions introduced by scattering particles in out-of-focus regions. Therefore, to isolate the effects of linear light scattering away from the focal volume, we have chosen to only consider non-resonant targets in a refractive index-matched medium. In future studies, we will examine index-mismatched objects in focus, as well as the effects of linear scattering from particles with more complicated geometrical shapes that are more representative of actual tissue structures.

Conclusion
We have presented a computational framework to simulate the propagation, generation, and detection of CARS signals in a medium containing scattering particles at deterministic locations. We utilize the HF-WEFS method to simulate the pump and Stokes excitation fields, which are distorted by the presence of scatterers in the propagation path toward the focal region. Using the perturbed excitation fields, we calculate the spatially-dependent third-order polarization of objects in the focal volume, and apply the radiating dipole approximation to calculate the far-field CARS radiation. This framework is applied to examine the effects of scattering on the far-field CARS radiation pattern from three index-matched objects for various scatterer sizes and locations under different illumination and detection conditions.
Our results demonstrate that the presence of small scattering objects proximal to the focal volume results in an attenuation of the CARS signal intensity. The signal attenuation can be minimized by using lenses with increased NA for both excitation and detection. In addition to signal attenuation, we have also observed significant distortions to the angular distribution of the CARS radiation. This effect can be related to the scattering-induced phase shifts imprinted by the excitation fields on the nonlinear polarization in the focal volume. The attenuation and propagation direction of CARS radiation is highly dependent on the size and position of the particle, an observation that underlines that the effects of light-scattering in a coherent technique like CARS microscopy are complex and require a full view of amplitude and phase distortions. Our computations confirm that placement of a scattering object near focus produces noticeable artifacts in the CARS intensity image. However, our simulations also show that CARS anisotropy images are much less sensitive to the presence of spherical scatterers.
The framework and results presented in this work provide a platform for detailed mechanistic study of the effects of light scattering on the quality of CARS images. With subsequent improvements of the model, including the consideration of multiple scatterers and scatterers of varying shape and refractive properties, we expect that the bottom-up understanding gleaned from these simulations will foster the development of adaptive optics strategies for coherent nonlinear optical microscopy.