Abstract
A rigorous forward model solver for conventional coherent microscope is presented. The forward model is derived from Maxwell’s equations and models the wave behaviour of light matter interaction. Vectorial waves and multiple-scattering effect are considered in this model. Scattered field can be calculated with given distribution of the refractive index of the biological sample. Bright field images can be obtained by combining the scattered field and reflected illumination, and experimental validation is included. Insights into the utility of the full-wave multi-scattering (FWMS) solver and comparison with the conventional Born approximation based solver are presented. The model is also generalizable to the other forms of label-free coherent microscopes, such as quantitative phase microscope and dark-field microscope.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
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]
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
which involve operations of convolution and correlation. Thus, fast Fourier transform (FFT) can be employed to save memory and accelerate computation. Details of FFT implementaion for half-space dyadic Green’s function can be found in [36]. To simplify the calculation of gradient, we use the contrast current $\mathbf {J}(\mathbf {r})=\chi (\mathbf {r})\mathbf {E}_t(\mathbf {r})$ [37], instead of the total field, as the unknown variable in the equation. It can be obtained by solving the least-square problem
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
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Funding
HORIZON EUROPE European Research Council (804233); Norges Forskningsråd (288082); Universitetet i Tromsø.
Acknowledgments
We acknowledge the generous help of Prof. Zicheng Liu for sharing the codes of the work [13] and other technical discussions.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
Supplemental document
See Supplement 1 for supporting content.
References
1. J. W. Lichtman and J.-A. Conchello, “Fluorescence microscopy,” Nat. Methods 2(12), 910–919 (2005). [CrossRef]
2. P. F. Gao, G. Lei, and C. Z. Huang, “Dark-field microscopy: recent advances in accurate analysis and emerging applications,” Anal. Chem. 93(11), 4707–4726 (2021). [CrossRef]
3. G. S. Verebes, M. Melchiorre, A. Garcia-Leis, C. Ferreri, C. Marzetti, and A. Torreggiani, “Hyperspectral enhanced dark field microscopy for imaging blood cells,” J. Biophotonics 6(11-12), 960–967 (2013). [CrossRef]
4. Z. Wang, L. Millet, M. Mir, H. Ding, S. Unarunotai, J. Rogers, M. U. Gillette, and G. Popescu, “Spatial light interference microscopy (slim),” Opt. Express 19(2), 1016–1026 (2011). [CrossRef]
5. C. J. Mann, L. Yu, C.-M. Lo, and M. K. Kim, “High-resolution quantitative phase-contrast microscopy by digital holography,” Opt. Express 13(22), 8693–8698 (2005). [CrossRef]
6. E. Cuche, F. Bevilacqua, and C. Depeursinge, “Digital holography for quantitative phase-contrast imaging,” Opt. Lett. 24(5), 291–293 (1999). [CrossRef]
7. I. Yamaguchi, J.-i. Kato, S. Ohta, and J. Mizuno, “Image formation in phase-shifting digital holography and applications to microscopy,” Appl. Opt. 40(34), 6177–6186 (2001). [CrossRef]
8. K. Kim, A. S. Somkuwar, Y. Park, and R. K. Singh, “Imaging through scattering media using digital holography,” Opt. Commun. 439, 218–223 (2019). [CrossRef]
9. Y. Park, C. Depeursinge, and G. Popescu, “Quantitative phase imaging in biomedicine,” Nat. Photonics 12(10), 578–589 (2018). [CrossRef]
10. T. H. Nguyen, M. E. Kandel, M. Rubessa, M. B. Wheeler, and G. Popescu, “Gradient light interference microscopy for 3d imaging of unlabeled specimens,” Nat. Commun. 8(1), 210 (2017). [CrossRef]
11. Y. Cotte, F. M. Toy, C. Arfire, S. S. Kou, D. Boss, I. Bergoënd, and C. Depeursinge, “Realistic 3d coherent transfer function inverse filtering of complex fields,” Biomed. Opt. Express 2(8), 2216–2230 (2011). [CrossRef]
12. L. Hu, R. Chen, K. Agarwal, C. J. R. Sheppard, J. C. Phang, and X. Chen, “Dyadic green’s function for aplanatic solid immersion lens based sub-surface microscopy,” Opt. Express 19(20), 19280–19295 (2011). [CrossRef]
13. Z. Liu and K. Agarwal, “Silicon substrate significantly alters dipole-dipole resolution in coherent microscope,” Opt. Express 28(26), 39713–39726 (2020). [CrossRef]
14. E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. 1(4), 153–156 (1969). [CrossRef]
15. A. Devaney, “Inverse-scattering theory within the rytov approximation,” Opt. Lett. 6(8), 374–376 (1981). [CrossRef]
16. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University, 1999).
17. T. Chang, S. Shin, M. Lee, and Y. Park, “Computational approach to dark-field optical diffraction tomography,” APL Photonics 5(4), 040804 (2020). [CrossRef]
18. X. Chen, M. E. Kandel, and G. Popescu, “Spatial light interference microscopy: principle and applications to biomedicine,” Adv. Opt. Photonics 13(2), 353–425 (2021). [CrossRef]
19. L. Tian and L. Waller, “3d intensity and phase imaging from light field measurements in an led array microscope,” Optica 2(2), 104–111 (2015). [CrossRef]
20. A. M. Maiden, M. J. Humphry, and J. M. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29(8), 1606–1614 (2012). [CrossRef]
21. S. Chowdhury, M. Chen, R. Eckert, D. Ren, F. Wu, N. Repina, and L. Waller, “High-resolution 3d refractive index microscopy of multiple-scattering samples from intensity images,” Optica 6(9), 1211–1219 (2019). [CrossRef]
22. M. Chen, D. Ren, H.-Y. Liu, S. Chowdhury, and L. Waller, “Multi-layer born multiple-scattering model for 3d phase microscopy,” Optica 7(5), 394–403 (2020). [CrossRef]
23. K. S. Kunz and R. J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics (CRC, 1993).
24. M. Totzeck and H. Tiziani, “Interference microscopy of sub-λ structures: A rigorous computation method and measurements,” Opt. Commun. 136(1-2), 61–74 (1997). [CrossRef]
25. M. Thomas, R. Su, N. Nikolaev, J. Coupland, and R. Leach, “Modeling of interference microscopy beyond the linear regime,” Opt. Eng. 59(03), 1–034110 (2020). [CrossRef]
26. T. Pahl, S. Hagemeier, M. Künne, D. Yang, and P. Lehmann, “3d modeling of coherence scanning interferometry on 2d surfaces using fem,” Opt. Express 28(26), 39807–39826 (2020). [CrossRef]
27. T. Pahl, S. Hagemeier, J. Bischoff, E. Manske, and P. Lehmann, “Rigorous 3d modeling of confocal microscopy on 2d surface topographies,” Meas. Sci. Technol. 32(9), 094010 (2021). [CrossRef]
28. M. Moharam, E. B. Grann, D. A. Pommet, and T. Gaylord, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” J. Opt. Soc. Am. A 12(5), 1068–1076 (1995). [CrossRef]
29. M. Totzeck, “Numerical simulation of high-na quantitative polarization microscopy and corresponding near-fields,” Optik 112(9), 399–406 (2001). [CrossRef]
30. P.-I. Schneider, P. Manley, J. Krüger, L. Zschiedrich, R. Köning, B. Bodermann, and S. Burger, “Reconstructing phase aberrations for high-precision dimensional microscopy,” in Optics and Photonics for Advanced Dimensional Metrology II, vol. 12137 (SPIE, 2022), pp. 120–127.
31. W. C. Gibson, The Method of Moments in Electromagnetics (Chapman and Hall/CRC, 2021).
32. M. A. Yurkin and A. G. Hoekstra, “The discrete dipole approximation: an overview and recent developments,” J. Quant. Spectrosc. Radiat. Transfer 106(1-3), 558–589 (2007). [CrossRef]
33. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University, 2012).
34. W. C. Chew, Waves and Fields in Inhomogeneous Media (Wiley-IEEE, 1995).
35. G. Gao, C. Torres-Verdin, and T. Habashy, “Analytical techniques to evaluate the integrals of 3d and 2d spatial dyadic green’s functions,” Prog. Electromagn. Res. 52, 47–80 (2005). [CrossRef]
36. T. J. Cui and W. C. Chew, “Fast algorithm for electromagnetic scattering by buried 3-d dielectric objects of large size,” IEEE Trans. Geosci. Remote Sensing 37(5), 2597–2608 (1999). [CrossRef]
37. P. M. Van Den Berg and R. E. Kleinman, “A contrast source inversion method,” Inverse Problems 13(6), 1607–1620 (1997). [CrossRef]
38. M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving,” Journal of Research of the National Bureau of Standards 49(6), 409 (1952). [CrossRef]
39. H. A. Van der Vorst, “Bi-cgstab: A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems,” SIAM J. Sci. and Stat. Comput. 13(2), 631–644 (1992). [CrossRef]
40. A. Ahmad, N. Jayakumar, and B. S. Ahluwalia, “Demystifying speckle field interference microscopy,” Sci. Rep. 12(1), 10869 (2022). [CrossRef]
41. A. Butola, D. Popova, D. K. Prasad, A. Ahmad, A. Habib, J. C. Tinguely, P. Basnet, G. Acharya, P. Senthilkumaran, D. S. Mehta, and S. A. Balpreet, “High spatially sensitive quantitative phase imaging assisted with deep neural network for classification of human spermatozoa under stressed condition,” Sci. Rep. 10(1), 13118 (2020). [CrossRef]
42. I. J. Cox and C. J. R. Sheppard, “Information capacity and resolution in an optical system,” J. Opt. Soc. Am. A 3(8), 1152–1158 (1986). [CrossRef]
43. V. S. Hatamzadeh, M. M. Naser, and S. R. Sadeghzadeh, “Numerical method for analysis of radiation from thin wire dipole antenna,” International Journal of Industrial Mathematics 3, 135–142 (2011).
44. B. B. Baker and E. T. Copson, The Mathematical Theory of Huygens’ Principle (American Mathematical Soc.,2003).
45. J. M. Jin, The Finite Element Method in Electromagnetics (John Wiley & Sons, 2015).