1. Introduction

Visualization of biological samples is an important task in biologic research. Fluorescence microscopy [1] is widely used to image biological samples where fluorescence markers or dyes are used to make the specimens visible. Although it provides detailed structures of the biological sample with high specificity, it suffers from some limitations such as photo-toxicity and photo-bleaching.

Label-free imaging solutions, which rely on the scattering due to the intrinsic contrast of the sample, are therefore preferred in the situations where long-duration imaging, live cell friendliness, not perturbing the biological system through chemical additions, etc. are preferred. Label-free microscopy has been developed for several decades, including bright-field microscopy, dark-field microscopy [2,3] and interferometric microscopy [4,5]. Among them, arguably, bright-field and quantitative phase imaging microscopy are the most popular qualitative and quantitative label-free imaging solutions respectively.

One of the challenges with label-free imaging solutions is that their images are more difficult to interpret than fluorescence microscopy images. Besides the lack of specificity, another important reason is that multiple scattering of light in the sample region contributes to a non-linear map of the sample in the image region, which appears artefacts while it is not. Solutions such as incoherent or partially-coherent illumination systems have been used to suppress the effect of multiple scattering in the image plane. However, the interest in coherent imaging solutions still remains because the evidence of multiple scattering in the image plane indicates potential of being able to decode the near field interaction and thereby seek super-resolution. This potential however remains untapped and the artefact-prone nature of label-free images remain uninterpreted because of lack of computational tools in microscopy domain that correctly simulate and explain how and how much of the near field multiple scattering interaction actually reaches the imaging plane. There is some preliminary work done. For example, scalar field propagation based models have been proposed, especially used in holography [6–8] and other quantitative phase imaging techniques [4,9,10]. Analysis of imaging system is mainly established by first-order diffraction theory and scalar point spread function [11]. 3D dyadic Green’s functions (or vectorial point spread functions) have been presented [12,13] which gives the mapping of dipoles induced in the sample with various orientations to a vectorial field and can be considered foundational for coherent imaging. Further, there have been some simulators that use Born [14] or Rytov approximation [15] for calculation [9,16–18]. Explicit solutions can be obtained for these linearized cases and ease the analysis. Several forward models considering z-propagation of multiple scattering for thick samples have been proposed and used for inversion process. Multi-slice is a popular method where the sample is decomposed of thin planes along optical axis and the field is calculated layer by layer with some linear approximation, such as multi-slice beam propagation methods [19–21] and multi-layer Born methods [22]. Such methods are definitely more sophisticated than Born approximation. However, multiple-scattering effect is not accurately modeled in these work.

An accurate forward model is essential for 3D reconstruction of the sample, as the reconstruction is usually relied on decreasing the discrepancy between the measured data and the simulated data obtained by the forward model. Due to the wave property of the light, the interaction between light and sample be can be precisely described by Maxwell’s equations. Based on the equations, forward modeling can be achieved by solving time-dependent differential equations in time domain with finite difference time domain (FDTD) method [23], or by solving integral or partial differential equations in frequency domain with numerical solvers such as finite element method (FEM), boundary element method (BEM) and so on [24–27]. Yet, these methods are limited to 2D cases where the object is assumed infinite and invariant along one axis. Rigorous coupled-wave analysis (RCWA) [28,29] can give an exact solution of Maxwell’ equations to simulate the electromagnetic diffraction by grating structures. However, its application is also limited to period strutures. A rigorous model is used in [30] for 3D objects with scattered field in near field solved by FEM and mapping each plane-wave component of the scattered field to the camera to get the intensity map. In this work, we present a rigorous computational model based on the integral method with Green’s function which describes the field radiated by a dipole as sum of plane-wave components for a full-wave three-dimensional case with temporally coherent illumination and it can be applicable to arbitrary 3D objects. Simulation of light propagation and field captured on camera is given for a conventional microscope setup with support for reflection mode bright field microscopy.

Here, we briefly discuss the physics being simulated. The light is considered as wave. When monochromatic light is used to illuminate the sample, the atoms in samples are polarized and become dipoles which are oscillating at the same frequency with the applied field. These dipoles also radiate electric field and such field is described by dyadic Green’s function. The field produced by one dipole generates a field which changes the total field elsewhere and further influences changes in the polarization of other dipoles. Therefore, the total field at one point in the sample is influenced by and also influences the field at all other points. Such phenomenon, called multiple-scattering effect, is described by Lippmann-Schwinger integral equation. The fields produced by all dipoles in the sample are added up to form an image on the camera. Another Dyadic Green’s function is derived for this specific microscope setup which gives the field produced by a dipole at the far field. Bright-field image and interferogram are formed as the intensity map of combination of the complex scattered field, reflected incident field and reference field. This model gives a full description of the observed image. However, it is impossible to find an explicit solution for random samples, thus, numerical techniques are employed to simulate this process. The method of moments (MoM) [31] and discrete dipole approximation (DDA) [32] are popular tools to get numerical solutions. Here MoM based on a uniform grid with pulse-basis function and point-matching method is used to ease computational cost. This method allows the calculation of far field at the camera under different configurations and investigation of near field which helps understand the effect of samples with different physical parameters. This simulation is also valid for strong-scattering objects and also shows the potential to develop new imaging algorithms based on inverse methods which can break the diffraction limit due to the integration of multiple-scattering effect in the forward model.

In Sec. 2, we demonstrate the forward problem describing the image field on the camera given the distribution of refractive index (RI) of the sample and how the problem is solved. In Sec. 3 we show some simulation results that are validated by experiments. In Sec. 4, we present further insight into the role of multiple scattering by comparing the simulation results with Born Approximation and also showing the near-field results inside domain of interest. In Sec. 5, we discuss the scope of generalization of our model to different microscopes and samples. Conclusion is given in Sec. 6.

2. 3D Full wave multiple scattering solver

The present forward solver solves the problem in frequency domain thus temporally coherent illumination is required. Here we take conventional bright-field microscopy with reflection measurement and spatially coherent illumination as an example. However, the model can be employed with various kinds of illumination. See Supplement 1 (Sec. S3) for more information about the adaptation of our model for other illuminations.

2.1 Problem description

The microscope setup under investigation exploits reflection measurement for bright-field imaging, as shown in Fig. 1(a). Collimated coherent light with angular frequency $\omega$ ($\omega =2\pi c/\lambda$ with $c$ the light speed and $\lambda$ the wave length in vacuum) is employed to illuminate the sample. Time dependence $e^{-i\omega t}$ is discarded for simplicity in the following formulation. The sample, immersed in a medium with refractive index $n_{obj}$, is placed on the substrate with a refractive index $n_{sub}$. The scattered light, together with the reflected incident light, propagates through the objective lens, tube lens and forms the bright-field image on the camera. Here we define two local coordinate systems with the focal points of objective lens and tube lens being the original points and $z$-axis is chosen as the optical axis as shown in Fig. 1(b). The Domain of Interest (DoI) is a cubic with size of $X\times Y\times Z$ in the upper half-space (object region) with all the samples enclosed. With the distribution of RI $n(\mathbf {r})$ at point $\mathbf {r}=(x,y,z)$ inside DoI known, the forward problem can be described by equations [34]

 

Fig. 1. (a) Schematic diagram of the microscope setup. For simulation, laser source is used to generate temporally and spatially coherent illumination and bright-field image is simulated with reflection measurement mode. Partially spatially coherent unit: in experiment, incoherent illumination is also employed (as described in Sec. 3) to suppress unexpected fringes. Add-on module for interference microscopy: generalization to off-axis interferometric phase microscopy can be achieved by including the reference arm (described in Supplement 1, Sec. S2). Note that modules in dashed blocks are not considered in the following modelling. MO: microscopic objective lens, BS: beam splitter, L: lens, RD: rotating diffuser. (b) Coordinate systems defined for simulation. Objective and tube lenses are represented by Gaussian reference spheres [33]. There origins are chosen as the origins of the local coordinate systems with the optical axis being the $z$-axis. DoI is a cubic region on the substrate inside which RI is a function of location based on the distribution of the samples.

Download Full Size | PDF

2.2 Numerical solver

To solve the problem numerically, the equations are discretized with the method of moments. The whole domain of interest is discretized into $N=N_x\times N_y\times N_z$ small cells with side $h$ and all the physical variables are considered uniform in each cell, represented by the value at the center of each cell $\mathbf {r}_m,m=1,2,\ldots,N$. Considering the singularity, the half-space dyadic Green’s function $\overline {\overline {G}}_d(\mathbf {r},\mathbf {r}')$ is split into two parts as $\overline {\overline {G}}_d(\mathbf {r},\mathbf {r}')=\overline {\overline {G}}_h(\mathbf {r},\mathbf {r}')+\overline {\overline {G}}_r(\mathbf {r},\mathbf {r}')$ with $\overline {\overline {G}}_h(\mathbf {r},\mathbf {r}')$ representing the primary field and $\overline {\overline {G}}_r(\mathbf {r},\mathbf {r}')$ the reflected field from the substrate. The state equation is discretized as follows

The discrete equation can be simplified as

We can rewrite it as matrix-vector product and the total field can be obtained by matrix inversion given the illumination and distribution of contrast. However the system matrix has a size of $3N\times 3N$ which leads to difficulty in large-scale problem. In our simulation, we use an optimization method to get the solution. Exploring the expressions of Green’s function, we can find that

Similarly, we discretize the image plane into $M=M_x\times M_y$ pixels centered at $\mathbf {r}_{cam,m},m=1,2,\ldots,M$. Once the contrast source $\mathbf {J}(\mathbf {r})$ is obtained, the scattered field at each point can be calculated as

2.3 Illumination and collected light

The illumination is assumed to be collimated, represented as plane wave

Given reflection coefficient $R=\dfrac {n_{obj}-n_{sub}}{n_{obj}+n_{sub}}e^{-2ik_{obj}z_{sub}}$ for normal incidence with $z_{sub}$ representing $z$-coordinate of the substrate, the incident field is the sum of the illumination and the reflected light as

The reflected light will propagate through objective and tube lens and is captured by the camera as

Here $f_{obj}$ and $f_{cam}$ are the focal lengths of objective and tube lens respectively. $k_{cam}$ is the wavenumebr in the image plane. $M$ is the magnification of the microscope. Note that we have neglected the intensity loss when light passes through the objective and tube lens and we also assume that the light intensity is uniformly distributed on camera after refraction by the lens. The field generated from the sample arm on the camera is $\mathbf {E}_{sample}=\mathbf {E}_{s}+\mathbf {E}_{r}$ and the bright-field image is

 

Fig. 2. Illustration of our 3D full-wave multi-scattering (FWMS) forward solver process. Example is shown with a 1 micron-diameter bead ($n=1.6$) placed on a substrate with $n_{sub}=3.5$ and in immersion medium ($n_{obj}=1$) illuminated by $x$-polarized plane wave ($\lambda$=660 nm) with normal incident angle. Far-field DGF is derived to simulate a 40 $\times$ 0.65 NA objective lens. With known incident field and distribution of the contrast, contrast source is determined first which indicates how the dipoles are induced in the DoI, as $J_x$, $J_y$ and $J_z$ shown in three orthogonal planes in a 3D cubic. Each contrast source component will radiate vectorial scattered field on camera, which gives $3\times 3$ 2D data maps with the row index showing contrast source component and column index indicating the component of the scattered field. The final scattered field is obtained by adding up the data maps by column which compose the same component of the field. The square of the amplitude of each component is computed as the intensity component and the final camera image of scattered field on camera is the sum of these three components. In the shown example, the $x$-component of the scattered field produced by $J_x$ (1st row, 1st col. in the $3\times 3$ 2D maps) has the highest amplitude and keeps the shape of the object. The final image on the camera is dominated by this component.

Download Full Size | PDF

3. Examples and experimental validation

Some simulation examples and comparison with experiments are shown in this section.

3.1 Microscope setup

The setup of the microscope is summarized as follows. One setting uses microscope system with 20 $\times$ objective with NA$=0.4$ and another setting uses microscope system with 40 $\times$ objective with NA$=0.65$. Both microscopes can be used in (a) incoherent illumination mode which is not relevant to the current modeling but is used to find the regions of interest, and (b) the coherent illumination model which is modeled in this paper. The coherent illumination has a wavelength of 660 nm in vacuum. The pixel size of camera is 6.5 µm. Since the camera is in air, $n_{cam}=1$ is set for image plane.

3.2 Sample details

The samples are placed on silicon wafer substrate ($n_{sub}=3.8$) and immersed in air ($n_{obj}=1$). 1 µm-diameter polystyrene beads ($n=1.6$) are used as samples.

3.3 Simulation configuration

For simulation, the DoI has a size of 4 µm$\times$4 µm$\times$1.1 µm. The image planes for 40 $\times$ and 20 $\times$ objectives are set as 156 µm$\times$156 µm and 104 µm$\times$104 µm, respectively. The DoI is discretized into small cells with side length of 50 nm, leading to $80\times 80\times 22$ voxels. The grid size in image plane is also 6.5 µm, thus the images have $24\times 24$ and $16\times 16$ pixels for two objectives respectively. Circularly-polarized illumination is obtained by setting $E_{x0}=1$ and $E_{y0}=e^{i\pi /2}$.

All the simulations are run with Matlab R2020a on a desktop computer (Intel Xeon W-2125 CPU @ 4.00GHz) with 128 GB RAM. The optimization problem is solved by built-in function of Bi-CGSTAB algorithm.

In order to make the experimental and simulated results comparable, the brightness/ contrast of the experimental results has been adjusted via ImageJ for the best visual match to the simulated images. This is done in consideration of the presents of noise in the experimental images and limited control in the experimental settings of the selection of $z$-plane precisely. In order to accommodate for imprecise knowledge of the experimental $z$-plane, we simulated three different locations of the substrate, i.e $z_{sub} =$ 0 nm, −500 nm, −1000 nm, with 0 being the position of the focal plane. Among the three simulations, we chose the one that presented the best visual match with the experimental results. We also note that no noise or speckle has been added to the simulated results, whereas our experimental present both noise and speckle. In the future works, realistic microscopy noise such as shot noise, camera’s dark current, and photoelectric noise can be included in the simulated images as desired.

3.4 Utility of incoherent illumination mode

There is a presence of undesirable circular fringes in coherent microscopy images, see for example Fig. 3(b) and Fig. 4(b). These circular fringes are referred to as speckle pattern or spurious fringes in conventional microscopy parlance [40,41] and is present due to high temporal and spatial coherent nature of the light source. The source of these spurious fringes is stray reflections from different surfaces of optical components that end up coherently interfering and finally reaching the camera.

 

Fig. 3. Experimental and simulation result of scattered field produced by two 1-micron beads with 20$\times$0.4 NA objective. Experimental results of (a) bright field image with incoherent illumination and (b) bright field image with coherent illumination are shown, respectively. (c) and (d) give the corresponding zoom-in images of the region of interest to be simulated. (e) gives the estimated positions of the two beads: (-960, 1020) nm and (1412.5, -862.5) nm. (f) shows the simulated bright field image. (g-i) show the normalized intensity vs. distances from the round points along the blue, pink and green line sections respectively drawn in (c). FWHW is estimated with linear interpolation for line sections crossing the first bead in the simulated image. The intensity profile in (h) has a larger FWHM than the profile in (g) which shows the ellipticity of the simulated bead.

Download Full Size | PDF

 

Fig. 4. Experimental and simulation result of scattered field produced by two 1-micron beads with 40$\times$0.65 NA objective. Experimental results of bright field with incoherent illumination and coherent illumination are shown in (a), (b) (original results) and (c), (d) (zoom-in images). (e) gives the estimated position of the beads: (-568.75, -81.25) nm and (893.75, 406.25) nm. (f) shows the simulation result of bright-field image.

Download Full Size | PDF

Due to the presence of speckle and spurious fringes in coherent illumination, it is difficult to spot small bodies such as beads in high contrast and less afflicted by the fringes. So, we also record bright field image with incoherent illumination and estimate the locations of the beads before acquiring the coherent illumination images. The result shown below is chosen from these locations which gives the best match with the experimental results.

3.5 Results for 20 $\times$ objective lens and two beads

For 20 $\times$ objective, the experimental results and the selected region to image is shown in Fig. 3. The locations of the two beads estimated from the incoherent image are also shown in the figure. The center to center distance between the two beads is about 3.03 microns. With incoherent illumination (Fig. 3(a)), the images of the two beads are regular circles in the experiment (Fig. 3(c)). However, the shapes become irregular in the coherent image (Fig. 3(b),(d)) due to speckle pattern.

Beside this effect, the images of the beads themselves illustrate important effects that are explained using our full field simulation. One observation is that the images of the two beads appear individually radially symmetric or circular in the incoherent image while they appear elliptic in both the coherent and simulated coherent images. In order to understand this better, we compared the normalized intensity profiles across a section connecting the two beads (Fig. 3(g)) and also individual beads (Fig. 3(h),(i)). It is seen that the widths of the curves across the individual beads is somewhat similar for incoherent, coherent, and coherent simulated images. However, the dip of intensity in the region between the two beads is less pronounced for the incoherent image. This implies that incoherent image presents lesser contrast than the coherent and coherent simulated images. The improved contrast can be easily explained using the coherent simulated images as the result of destructive interference of the scattered electric fields resulting from the induced sources on the beads. This also indicates the potential of better contrast and better resolution in coherent imaging when destructive interference is observed, which is in line with the previous known theories in coherent imaging [42]. We also evaluated FWHM (full width at half minimum) to show the elliplicity of the simulated beads. The FWHMs of the intensity profiles along two approximately perpendicular lines in the first bead (located at (−960 nm, 1020 nm)) are estimated as 1771.4 nm (blue line in Fig. 3) and 2018.0 nm (pink line) which means the shape of the bead is squeezed along the line connecting the centers of the two beads. For the second bead, FWHMs are evaluated as 1716.8 nm (blue line) and 2016.9 nm (green line) which also demonstrates the elliptic shape of the simulated bead.

Though the normalized intensity profiles along the selected lines show good agreement between the incoherent, coherent and simulated results, there is visual mismatch between the simulated camera image and the experimental coherent image due to the unavoidable noise which leads to an obvious deformation of the beads in the coherent image and the estimation error in the location of the beads.

3.6 Results for 40 $\times$ objective lens and two beads

With 40 $\times$ objective, the selected imaging region and locations of two beads are shown in Fig. 4. The estimated distance between the beads is about 1.54 microns (Fig. 4(e)) computed using the incoherent image (Fig. 4(c)) of the beads identified as a region of interest in the original full field-of-view on incoherent image (Fig. 4(a)). The undesirable but unavoidable speckle pattern in the form of circular fringes (Fig. 4(b)) is present as expected in the experimental results with coherent illumination, which contribute to the distortion of the images of the beads (Fig. 4(d)). However, it is interesting to note that the image of the beads in the coherent image appears as connected at the junction of the two (see the yellow arrow), which gives the image a dumbbell-like appearance. This effect is also observed in the simulated coherent images (Fig. 4(f)) and is explained using the constructive interference of the fields created due to the induced currents on the beads. This example further shows that not only our simulated results match the coherent imaging results with sufficient fidelity to observe important optical effects, but it can also explain the observations. The bright-field coherent imaging result is the interference result of the scattered fields generated by the two beads and the reflected illumination. Since interference inherently involves superposition of the vector form of electric fields, it is important to take vectorial property of the fields and the polarization effect of the lens into consideration. The present solver has considered these effects: expression of fields as vector, dyadic Green’ function concerning orientation of dipole, TE ($s$-) and TM ($p$-) polarized wave decomposition in the Green’ function for the calculation of the reflection coefficient and focusing of the lens. These considerations have helped to give a more accurate model for simulation.

4. Further insights from FWMS solver

In this section, we use our forward solver to show some multiple-scattering effects inside one object and between two objects. Here we only focus on the scattered fields on the image plane, but we also show the physical variables, contrast source $J$ inside DoI which helps us analyze the results (corresponding total fields are shown in Supplement 1, Sec. S1). The results are compared to those obtained with conventional linear method where Born approximation (BA) is used and that the total field is simply replaced by incident field in DoI.

In this section, 40 $\times$ 0.65NA objective is used and the focal plane is set aligned with the center of the bead in all experiments. To explore the effect of polarization, the illumination is $x$-polarized by setting $E_{x0}=1$ and $E_{y0}=0$. The size of DoI is 3 µm$\times$3 µm$\times$2.1 µm. Other configurations are kept unchanged with the experiments in previous section.

4.1 Case 1: one object

In the first study we show the multiple-scattering effect even while simulating the image of only one object. First we show the simulated image of a 1-micron bead with $x$-polarized illumination, and then we investigate how the size and refractive index of the object affect the image. Our simulation result of the scattered field intensity produced by a 1-micron bead is shown in Fig. 5. For comparison, the result calculated with conventional linear BA and the relative difference is also shown in this figure.

 

Fig. 5. Simulation result of scattered field intensity produced by a 1-micron bead with $n=1.6$ calculated by FWMS and BA. The relative difference of the intensity at each pixel is also shown in the figure.

Download Full Size | PDF

Comparing the magnitude of the results, the Born Approximation gives a much higher intensity than ours, which means BA overestimates the scattered field, as also observed in [22]. The intensity in the BA result is focused at the center and decreases quickly so it gives a smaller size than the result given by FWMS. Though the bead is a sphere, an elliptic profile is observed in our simulation result. In contrast, BA generates a circular profile. We can further explore the result by visualizing the contrast source $J$ inside the object simulated by FWMS and BA. Here we show the $x$, $y$ and $z$ components of the magnitude and phase of induced current at the plane $z=$−25 nm in Fig. 6. For a one-to-one comparison, the top half of the plotted magnitude and phase maps shows the entities computed by FWMS and the bottom half shows the same entities computed by BA.

 

Fig. 6. Amplitude $|J_c|$ and phase $\angle{J_c}$ (degrees) of contrast source component $J_c$ at plane $z={-25}\,\text{nm}$ simulated by FWMS and BA method, separated by dashed line. Red circles represent the size of the bead on this cross section. The phase of $y$ component of the contrast source is symmetric about the circle center while others are symmetric about the dashed line. With BA, only $x$ component of contrast source is induced and shown under $x$-polarized illumination. In contrast, all the three components $J_x$, $J_y$ and $J_z$ are generated and shown.

Download Full Size | PDF

Born Approximation simply replaces the total field by incident field. As the illumination is $x$-polarized and the bead is homogeneous, the induced current only has $x$-component and distributes uniformly in the bead, as shown in the figure. $J_x$ has the same phase with the incidence, about 5.6 degrees at the selected plane. However, if we consider multiple scattering such as in FWMS, these induced $x$-polarized dipoles radiate and generate electric fields with $x$-, $y$- and $z$- components. These radiated fields change the fields that the dipoles sense and further change the polarization of the dipoles. Since the contrast of the bead is $\chi =(n^2-n^2_{obj})/n^2_{obj}=1.56$ and size of the bead is about $1.52\lambda$, the object is not weakly scattering. Such process should not be neglected and BA cannot give the correct result.

Our forward solver takes this multiple-scattering effect into consideration. As formulated by the state equation, the total field at each point is affected by all the fields inside the object which appear in the integrand in the right-hand side. In our simulation result, all the $x$-, $y$- and $z$-polarized dipoles are induced. Since the illumination is $x$-polarized, the $x$-component of the contrast source has a higher amplitude than other components. Besides, the distribution of the contrast current is not uniform in the bead, both in terms of amplitude and phase. The scattered field on the camera is mainly determined by $J_x$ and the $x$-component is dominant in the total intensity. Since $J_x$ is not uniformly distributed and has two regions with a higher amplitude distributed vertically, thus the final shape of the imaging result is elongated along $y$-axis.

Next we change the size and the contrast of the sample to investigate the influence of these parameters on the multiple-scattering effect.

First we show the simulation result of a 200 nm-bead whose refractive index is kept unchanged. Since the bead has a diameter of 200 nm, we use a smaller DoI and a finer discretization for simulation. The grid size is 10 nm and simulation results are shown in Fig. 7. Even though the bead is 200 nm large, it occupies more than 6 pixels in both results obtained by our solver and BA, which is consistent with diffraction theory. The deformation in the image of a 200 nm-bead is hardly seen and the shape is quite similar with the result given by BA. Since this sample has a smaller size, the estimated scattered field intensity becomes lower. The discrepancy between the amplitude of the intensity images is smaller but BA still overestimates the intensity at a percentage of 50 approximately. The multiple scattering effect in the single object has decreased with a smaller size.

 

Fig. 7. Simulation result of scattered field intensity of FWMS and BA and the relative difference produced by a 200 nm-bead with n=1.6.

Download Full Size | PDF

We can also investigate the contrast source inside the object as shown in Fig. 8. BA estimates a uniformly distributed contrast source in the object which is proportional to the incident field at a specific $xy$-plane. The $x$-component of the contrast source calculated by FWMS also is more homogeneous in both amplitude and phase. The magnitude of the contrast source is quite similar with the result of 1-micron bead and with BA, which decreases slightly from the center to the edge. The phase difference between the contrast source and the incidence also decreases slightly from the center to the edge.

 

Fig. 8. Amplitude $|J_c|$ and phase $\angle{J_c}$ of contrast source component $J_c$ at plane $z=-{5}\,\textrm{nm}$ simulated by FWMS and BA method. Red circles represent the size of the bead on this cross section. Since the DoI has $25\times 25$ pixels on this plane, the dashed line is placed on $y=-{5}\,\textrm{nm}$ with the center of the beads as the original point. The dominant component $J_x$ still has a high magnitude but becomes homogeneous. Phase $\angle{J_x}$ calculated by FWMS is not totally but nearly homogeneous while thr result given by BA is totally homogeneous as the the phase of incident field.

Download Full Size | PDF

Now we consider the refractive index of the object. In this experiment, we use a 1 µm-diameter bead with RI $n=1.05$, thus the contrast $\chi$ is $\sim$0.1 compared with 1.56 when RI is $n=1.6$. The intensity of scattered field and the induced current are shown in Fig. 9 and 10 respectively. The images produced by BA and FWMS are most similar to each other in both shape and amplitude in comparison to 200 nm small bead with contrast 1.56 and 1-micron bead with contrast 1.56. The discrepancy between BA and FWMS is quite small to be visibly obvious and multiple scattering is not obviously affecting the simulation result of FWMS. If we investigate the contrast source, we observe that BA still gives a uniformly distributed contrast current inside the bead but at a smaller magnitude due to a small contrast $\chi$, compared with Fig. 6. The electric field produced by these dipoles are small and they can affect the incident field in a limited way thus the distribution of the dipoles and the total field do not change a lot. Therefore unsurprisingly, in FWHM the $x$-component of the contrast source is more uniform and has a similar amplitude with BA result. $y$- and $z$- components with smaller magnitude are induced. The phase difference between the contrast source and incidence is also smaller compared with $n=1.6$ bead. Still, at the center of the bead, the phase difference is highest.

 

Fig. 9. Simulation result of scattered field intensity of FWMS and BA and the relative difference produced by a 1-micron bead with $n=1.05$.

Download Full Size | PDF

 

Fig. 10. Amplitude $|J_c|$ and phase $\angle{J_c}$ of contrast source component $J_c$ at plane $z=$25 nm simulated by FWMS and BA method. Red circles represent the size of the bead on this cross section. The amplitudes of all the components become smaller compared with $n=1.6$ bead. The main $x$ component becomes more homogeneous and similar with the BA simulation result.

Download Full Size | PDF

From the above discussion, we see that our forward solver gives a similar result with Born Approximation when the size or the refractive index of the object is small, which is the case when the BA is valid. On the other hand if the sample is optically large (order of wavelength) and has a high contrast compared to the background, Born approximation cannot give the correct result and FWMS is more accurate since multiple-scattering effect is considered.

4.2 Case 2: two objects

In this section we discuss the multiple-scattering effects between two objects. Since the importance of incorporating multiple scattering is obvious in this case, we only discuss the FWHM results. Distance and orientation of two samples are set as variables under investigation. First, two beads with $n=1.6$ are placed along $x$-axis, and the distance between their centers are $d=$ 1000, 1500 and 2000 nm. The scattered field intensity of these three cases simulated using FWMS is shown in Fig. 11.

 

Fig. 11. Simulation result of total scattered field intensity produced by two 1-micron beads placed on $x$-axis with distances $d =$1000 nm, 1500 nm and 2000 nm (from left to right) under $x$-polarized illumination.

Download Full Size | PDF

When the two beads are placed on $x$-axis, two peaks representing the two beads can be always identified even if the two beads are touching each other (when the center-to-center distance is 1000 nm). As the distance between the two beads increases, it is easier to distinguish them from the intensity map.

The contrast source induced in the two beads are shown in Fig. 12. Compared with Fig. 6 which gives the induced current in one bead, we can find that the distribution of the induced current has changed when the two beads are close to each other. As the two beads separate further, the contrast source distribution become similar to the single bead.

 

Fig. 12. Simulation result of induced current $J$ at $z={-25}\,\textrm{nm}$ with two beads placed on $y$-axis. From left to right: bead distance $d = {1000}\,\textrm{nm}$, ${1500}\,\textrm{nm}$ and ${2000}\,\textrm{nm}$. From top to bottom: components of induced current $|J_x|$, $|J_y|$ and $|J_z|$.

Download Full Size | PDF

Since the illumination light is $x$-polarized, we change the position of the beads to be located on $y$-axis to see if the orientation of the illumination has an influence on the beads placed along different directions. The distances between the two beads are still 1000 nm, 1500 nm and 2000 nm. The results are shown in Fig. 13 and 14.

 

Fig. 13. Simulation result of total scattered field intensity produced by two 1-micron beads placed on $y$-axis with distances $d =$1000 nm, 1500 nm and 2000 nm (from left to right) under $x$-polarized illumination.

Download Full Size | PDF

 

Fig. 14. Simulation result of induced current $J$ at $z={-25}\,\textrm{nm}$ with two beads placed on $y$-axis. From left to right: bead distance $d = {1000}\,\textrm{nm}$, ${1500}\,\textrm{nm}$ and ${2000}\,\textrm{nm}$. From top to bottom: components of induced current $|J_x|$, $|J_y|$ and $|J_z|$.

Download Full Size | PDF

It is difficult to identify the two beads when they touch each other, which is different from the previous case. If we explore the contrast source inside the beads, we can find that the component $J_x$ differ a lot with the result simulated in a single bead as in Fig. 6 even when the two beads are not so close to each other. This can be explained by the radiation pattern of a dipole, whose envelope for points of equal radiation intensity has a doughnut type shape, with the axis of the dipole passing through the hole in the centre of the doughnut [43]. For example, the field radiated by a dipole decays much faster along $x$-axis when the dipole is $x$-polarized and also locates on $x$-axis (as the weak radiation region in gray in Fig. 15). Thus, with $x$-polarized illumination, $x$-polarized contrast source dominates in the bead, but affects little the field on the $x$-axis. In contrast, when the two beads are vertically placed, the fields are influenced mutually and coupled with each other. The multiple scattering effect is stronger in this case which leads to difficulty in distinguishing these two beads in the scattered field. Figure 10 shows the influence between two dipoles with different locations.

 

Fig. 15. Illustration of the mutual influence between two dipoles. Black arrows indicate the polarization of the dipole. The points on the colorful surface have the same intensity of radiation. The radiation intensity of the dipole is low in the gray zones. In the left configuration, two dipoles are placed aligned with the orientation direction, i.e., in the weak-radiation of region of each other, thus they will be influenced weakly by each other. In contrast, the two dipoles are influenced by the radiation fields strongly when they are placed vertically (right configuration).

Download Full Size | PDF

5. Generalizability of FWMS model

5.1 Generalizability to other coherent microscopes

In this paper, we present the complete derivation of the model of bright-field microscopes. However, our model can easily be extended and generalized to most other coherent microscopes. We present the methods and examples of generalization to quantitative phase microscope (QPM) and dark-field microscope in the Supplement 1.

5.2 Generalizability to large realistic samples

It is of considerable interest to simulate large realistic samples such as 2D and 3D cell-cultures, but the current approach requires using voxel sizes small in comparison to the wavelength as well as the resolution of the system, which is at present a big compute problem even for the state-of-the-art computation hardware. Nonetheless, technically given sufficient computation resources, our model is generalizable to large realistic samples. Yet, it may be more useful to use surface based approaches which assume homogeneous interior of the objects being simulated and employs Huygen’s equivalence principle [44] for computing equivalent sources on the boundaries of the objects. In such a case, the adaptation of our model will be relatively straight forward, where a replacement of the multiple scattering model assuming voxel representation to multiple scattering model assuming boundary representation is used.

It may also be of interest to investigate adaptations of other modeling approaches such as finite elements methods [45] etc. for microscopy image simulation where versatile elements are generated to adapt to the complicated structures of biological samples. When simulating large scale structures, it might also be of interest to simulate multiple scattering of different orders depending on the vicinity of the elements for which multiple scattering is being simulated.

6. Conclusion

We demonstrate a 3D full wave multiple scattering (FWMS) forward model solver to simulate the image on camera for a coherent microscope. It simulates the wave behaviour of light in the sample and on the camera incorporating multiple scattering and reflections from surfaces, surpassing the limits of ray modeling, Fourier optics based modeling, and Born approximation based modeling of the coherent optical microscopes. This shall facilitate in understanding the challenges of imaging and identifying the opportunities for improvement in the microscope by checking the currently most rigorous simulation results generated by FWMS solver and physical variables inside region of interest and the microscope optical elements. Our forward model is generalizable to facilitate calculation of the scattered field and intensities, bright field images, as well as other conventional coherent microscopy setups. We validate the solver against the experimental brigh-field images and discuss how the solver can be used to understand the image characteristics. We also show that FWMS solver gives the similar results with Born Approximation when the scatterer is weak-scattering, in other words, when BA is valid. When the sample presents strong scattering such that BA is not valid, multiple scattering effect cannot be neglected and our FWMS solver fills the gap.

An important challenge of our FWMS solver is the computational cost. The computational cost mainly comes from calculation of Green’s function, where nine entries containing Sommerfeld integrals need to be calculated for mapping from one point to another. For example, it takes about 12 hours to calculate the far-field Green’s function used for 20 $\times$ objective by parallel computation in Matlab with four kernels, which has $9\times 140800(80\times 80\times 22)\times 256(16\times 16)\approx 324$ million entries. However, this can be decreased to about 3 hours if the symmetry of Green’s function is considered. Near field Green’s function $\overline {\overline {G}}_r(\mathbf {r},\mathbf {r}')$ can be used for all the configuration as long as the material of immersion and substrate are fixed. In contrast, far-field Green’s function $\overline {\overline {G}}_s(\mathbf {r}_{cam},\mathbf {r}')$ must be calculated again when the location of the substrate changes. Also, it needs to be recomputed for any change in the collection path of the microscope. On the other hand, if the collection path includes another interferometric component, such as in quantitative phase microscopy, then superposition of fields is sufficient and the far-field Green’s function need not be recomputed. Computational time for solving the forward problem varies from sample to sample, depends on how strong the sample is scattering. For example, in our simulation, it takes about 287 seconds to solve the scattering problem with one 1-micron bead with $n=1.6$ and 25 seconds for the bead with $n=1.05$. Implementation of FFT can greatly ease the computation burden, and this will be one of our priorities in the near future.

It is also worth mentioning that a full-wave forward model is important for designing imaging or reconstruction algorithms that can break the diffraction limit. Following the inverse theory, 3D distribution of physical parameters inside domain of interest may be retrieved by decreasing the discrepancy between the measured data and simulated data with present contrast iteratively, which is the focus of our future work.