Coupled and uncoupled dipole models of nonlinear scattering

: Dipole models are one of the simplest numerical models to understand nonlinear scattering. Existing dipole model for second harmonic generation, third harmonic generation and coherent anti-Stokes Raman scattering assume that the dipoles which make up a scatterer do not interact with one another. Thus, this dipole model can be called the uncoupled dipole model. This dipole model is not sufficient to describe the effects of refractive index of a scatterer or to describe scattering at the edges of a scatterer. Taking into account the interaction between dipoles overcomes these short comings of the uncoupled dipole model. Coupled dipole model has been primarily used for linear scattering studies but it can be extended to predict nonlinear scattering. The coupled and uncoupled dipole models have been compared to highlight their differences. Results of nonlinear scattering predicted by coupled dipole model agree well with previously reported experimental results. to have a maximum value of unity. This helps us to show the distribution of induced second order polarizations but not their relative magnitudes.


Introduction
Laser scanning nonlinear microscopy has become an important tool for imaging biological specimen. The principle of laser scanning nonlinear microscopy was first demonstrated in 1978 [1] and further developed into its current form in 1990 [2]. There are two forms of contrast in nonlinear optical microscopes -absorption followed by fluorescence, and scattering. In nonlinear fluorescence, two or more photons from the excitation laser are used up in exciting a fluorophore. On the other hand, in nonlinear scattering two or more photons from an excitation laser combine with the assistance of the sample to give a scattered photon. Here we are concerned mainly with nonlinear scattering. The three different types of nonlinear scattering which have been commonly used in imaging biological specimen are second harmonic generation (SHG), third harmonic generation (THG) and coherent anti-Stokes Raman scattering (CARS). In SHG, two photons of same frequency combine together with the help of the sample to give a single scattered photon of double the frequency. THG is very similar to SHG but it involves three excitation photons and a scattered photon of triple the frequency. CARS involves three excitation photons one each from the pump, probe and Stokes beams. The frequency difference between the pump and the Stokes photons overlaps with a vibrational frequency in the sample and hence there is resonance. The third excitation photon probes this resonance. All these nonlinear scatterings require samples with specific geometries for the process to be efficient. For example, SHG is strong from material lacking inversion symmetry like type I collagen fibers. SHG microscopy has been used for structural imaging of specially stained cell membranes [3], extracellular matrix [4], myosin [5] and collagen [6][7][8] in tissues. THG occurs primarily where there is a change of refractive index but it tends to cancel out in a focused beam if the sample is uniform. THG microscopy has been used to study the morphology of biological specimens [9,10]. CARS targets specific bonds in a sample. CARS microscopy has been used to image live cells [11] and lipid droplets [12] within cells.
While nonlinear scattering is a useful source of inherent contrast for imaging biological samples, interpretation of the images is not always simple. The intensity of the scattered light is obtained from the vectorial sum of electric field scattered by small regions within the focal volume of the excitation laser. Therefore the intensity of the scattered light is proportional to the scatterer's concentration raised to an exponent whose value can lie anywhere between -∞ to 2. Heterogeneities in the spatial distribution of scatterers in a sample make this process more complicated. Directionality of the scattered light also depends on the distribution of the scatterers. There are a few numerical methods like finite-difference time-domain (FDTD) method [13], finite element method (FEM) [14] and boundary element method (BEM) [15] which can be used to study nonlinear scattering from different samples. Dipole models are one of the simpler and intuitive numerical models for understanding complex nonlinear scattering processes. In these models, a sample is assumed to be made up of dipoles. The scattered light is calculated as a vectorial sum of radiations scattered by each of these dipoles. Mertz and associates proposed a dipole model for SHG and they were able to explain the SHG from different samples like giant unilamellar vesicles (GUVs) [16]. This model has also been extended to describe THG [17] and CARS [18]. This dipole model however has certain limitations. The dipoles in a scatterer interact with one another and with the incident wave. The final distribution of the excitation field is a result of this interaction. Similarly when nonlinear dipoles are induced, they interact among themselves and redistribute the net nonlinear field. This interaction between the dipoles is not captured by the dipole model we have discussed so far. Therefore this model can be aptly referred to as the uncoupled dipole model (UDM). Coupled dipole model (CDM) which takes into account interaction between dipoles was proposed by Purcell and Pennypacker [19], and it was further developed by Draine and associates [20][21][22]. CDM was developed for linear scattering but it can be extended to describe nonlinear scattering as well [23]. Although CDM is computationally expensive as compared to UDM, it gives results which are closer to experimental observations. Here we have compared the two models and shown how the results from the two models diverge with changes in the sample.

Theory
The approach of both the models, UDM and CDM, is the same. Any scatterer can be approximated as a collection of dipoles. An external electric field (E inc ) will induce dipole moments in all these dipoles. The dipole moment of each dipole can be calculated from the net local electric field. UDM is a simplified version of CDM and therefore we can start with CDM and reduce the results to get the UDM. In the case of CDM, if E inc,i is the incident electric field of angular frequency ω at the i th dipole and P i (1) is the resulting induced linear dipole moment, then these quantities are related by Eq. (1).
In the above set of equations, α ω,i is the linear polarizability of the i th dipole, N is total number of dipoles in the scatterer, k ( = ω/c) is the wave number of the incident electric field and r ij is the displacement between dipoles i and j. Electric field at the i th dipole due to the j th dipole is given by Eq. (2) and this accounts for dipole-dipole interaction. Equation (1) can be reduced to a simplified form (Eq. (3)) and solved for induced dipole moments. Equations (1)-(3) describe linear interaction between light and matter. Superscripts with parenthesis in the equations denote the order of interaction.
The linear local field (E loc,i ) at each dipole drives the nonlinear interactions and linear dipole moment of a dipole gives a good estimate of this linear local field (Eq. (4). The driving fields for various nonlinear interactions like SHG (Eq. (5)), THG (Eq. (6)) and CARS (Eq. (7)) were calculated. But these driving fields alone do not lead to the final nonlinear dipole moments. Interaction among all the nonlinear dipoles has to be taken into consideration for calculating the induced nonlinear dipole moments (Eq. (8)). β i and γ i are the first and second hyperpolarizabilities respectively of the i th dipole. α 2ω,i , α 3ω,i and α ω cars,i are polarizabilities of the i th dipole at the second harmonic, third harmonic and CARS wavelengths respectively. The term A ij (n) in Eq. (8) is defined in the same way as it is in Eq.
(2) but this definition will be based on the final scattering wavelength and not on the incident wavelength. Equations (3) and (8) involve large matrices and hence it is difficult to solve these equations by matrix inversion. These equations were therefore solved using conjugate gradient method. On removing the dipole interaction terms, CDM is reduced to UDM. Without dipole interactions, all the induced dipole moments can be calculated from the incident electric field (Eqs. (9)-(12)). (1) , (2) , , , Once the dipole moments are known, nonlinear scattered light in any direction can be calculated by vectorial summation of contributions from each dipole. In the case of microscopy, the signal is integrated over forward or backward collection aperture. All these numerical calculations were implemented in a MALAB code [23].

Results and discussion
The refractive index of a material is a measure of the interaction between the material and light. The greater the magnitude of refractive index of a material, the stronger is the interaction of light with that material. UDM takes into account only the shape of the scatterer but ignores its linear optical response. CDM, on the other hand, takes into account the shape as well as the optical properties of a scatterer. Therefore it was hypothesized that the results predicted by UDM and CDM should differ with an increase in the refractive index of a material. To test the hypothesis, second harmonic generation from hypothetical actin sample was analyzed. Second harmonic polarization is induced due to cross polarization terms in a thin layer of actin fiber bundle when it is excited by a focused beam of light [24]. When linearly polarized light, say x-polarized, is focused through a high numerical aperture (NA) lens, the polarization of light at the focus is not entirely along the x-axis and there is significant power in y-and z-polarized components as well [25]. In the hypothetical experimental set-up (Fig. 1), x-polarised light is focussed onto a thin layer of actin fiber bundles oriented along x-axis. The wavelength of the excitation source is 800 nm and the NA of the objective is 0.87 ( = sin60°). The vectorial distribution of electric field at the focus of a high NA lens was calculated [25]. In terms of magnitude, the strongest component of excitation field is E x , followed by E z and E y . The driving field for second order dipoles in the actin fibers was calculated (Eq. (13)) from the nonlinear susceptibility (χ (2) ) of actin fiber bundle [26]. The dipole size was set to 13 nm for both UDM and CDM simulations.
The second order polarization induced within the focal region of the lens can be calculated using UDM or CDM. In the case of UDM, induced second order polarization can be calculated without knowing the refractive index of the sample. But in the case of CDM, refractive index of actin fibers must be known. For the CDM calculations, refractive index of the sample was set to either 1.42 or 1.60. Such large values of refractive index were chosen to highlight difference in results predicted by the two models. Distribution of x-component of the second order polarization (P x (2) ) in the focal plane shows three different lobes (Fig. 2). The two side lobes are contributions from the E z component of the local field. The central lobe comes from the E x component of the local field. E y being the weakest component, its contribution to P x (2) is very small. Results calculated using UDM (Fig. 2(a)) show that the central lobe is very weak as compared to the side lobes. On the other hand, results predicted using CDM (Fig. 2(b) and 2(c)) show the central lobe becomes stronger with increase in refractive index of the actin fibers. The results show that with increase in refractive index, the E x component of the local field becomes stronger as compared to the E y and E z components. Any change in refractive index of the sample should affect the scattering from the sample and CDM predicts how it happens in this case.
The other components of the induced second order polarization, P y (2) and P z (2) depend on the products E x E y and E x E z respectively. Even if the E x component of local field becomes stronger with increasing magnitude of refractive index, the spatial distribution of the products E x E y and E x E z does not change significantly (Fig. 3). Here it should be noted that each plot in Figs. 2 and 3 has been normalized to have a maximum value of unity. This helps us to show the distribution of induced second order polarizations but not their relative magnitudes. Like UDM, CDM can be extended to other types of nonlinear scattering like THG and CARS. Scattering from beads was used as an example to demonstrate how the two models differ in their predictions. THG from a 1.5 µm diameter polystyrene bead was taken as an example for this study. Refractive index of polystyrene was taken to be 1.59 [27]. The beads were assumed to be surrounded by water. Excitation wavelength was 1200 nm and the numerical aperture of the objective was set to 1. The size of dipoles was set to 25 nm for both UDM and CMD calculations. Third order susceptibility (χ (3) ) of polystyrene was assumed to be the same as that of an isotropic material [28]. Forward scattered THG was calculated using UDM and CDM as the bead was scanned axially through the focus of the excitation beam. Since THG is strong at interfaces, axial scanning should give strong signals at the front and back surface of the bead. From the center of the sphere, the surface in the direction of incident light is referred to as the front surface. In the case of UDM, both these surfaces are symmetric and so is the signal generated by them (Fig. 4). In reality, the bead interacts with the excitation field and therefore the field distribution is different when light is focused at the front and back surfaces of the bead. CDM is able to capture this effect because it takes into account the presence of sample dipoles when calculating the field distributions. Therefore signal from the two interfaces is different (Fig. 4). CDM results agree well with experimental observations reported in literature [29]. Forward and backward CARS from a polystyrene bead of diameter 1.5 µm were calculated using CDM and UDM. In practice, the probe wavelength is the same as the pump wavelength. The pump/probe wavelength (λ p ), the Stokes wavelength (λ s ) and the CARS wavelength (λ c ) in this study are 750 nm, 852 nm and 670 nm respectively. This corresponds to a Raman shift of 1600 cm −1 in polystyrene beads [30]. The beads were assumed to be surrounded by water. The size of dipoles was set to 25 nm for CDM calculations. The γ of polystyrene should have resonant and a non-resonant component whereas the γ of water should have only the non-resonant component. The non-resonant component of γ was taken to be of the same form as that in THG calculations. The γ of water, which contributes to background here, is assumed to be 60% in magnitude of the non-resonant γ of polystyrene. For the resonant γ, the ratios of non-zero elements are γ xyyx / γ xxxx = 3/4, γ xxyy / γ xxxx = 1/8 and γ xyxy / γ xxxx = 1/8. Furthermore, the resonance of γ was taken to be 2.5 times the magnitude of non-resonant γ . Being a third order scattering process, CARS is present in all materials. The effect of the Gouy phase shift at the focus is partially compensated by the interaction between the pump field and the conjugate Stokes field [31]. Therefore CARS is not restricted to interfaces. CDM results for CARS during axial scan of the polystyrene bead show that the locations of peaks in forward and backward CARS do not overlap (Fig. 5(a)). These results match very well with experimental results reported earlier and this effect has been attributed to refractive index mismatch between the bead and its surroundings [32]. UDM results for CARS during axial scan of the bead show that the forward and backward CARS are of completely different nature (Fig. 5(b)). Forward CARS shows a single peak at the position of the bead whereas backward CARS shows two peaks at the front and back surfaces of the bead. These UDM results do not fit in with experimental observations [32].

Conclusions
A comparison between UDM and CDM shows that the latter gives better results when the refractive index of the scatter is large in magnitude. CDM is also able to predict the field enhancement at interfaces where there is a change of refractive index. The main difference between the two models is the interaction of different parts of the scatterer among themselves. This interaction becomes stronger with increase in the magnitude of refractive index. However, for scatterers with low magnitudes of refractive index, like biological specimen, the far field scattering results as predicted by UDM and CDM do not differ significantly. Under such circumstances UDM is a better choice because of its lower computational cost as compared to CDM. In experiments involving strong scatterers like metal nanoparticles or in experiments where surface effects are important, CDM gives results which are more close to experimental observations. Speed of CDM for linear scattering has been reported extensively in literature [22]. CDM for SHG and THG takes twice the time required by CDM for linear scattering. For CARS, CDM needs thrice the amount of time required for linear scattering. CDM will be very useful in this era of nanophotonics where optical properties of small particles have assumed great importance.