1.

Introduction

Imaging in thick samples is key to understanding more complex samples, such as developing embryos, organoids, or tissues. Optical imaging is a powerful tool for quantitative and functional imaging but has traditionally struggled with thicker samples where scattering dominates. For example, imaging of zebrafish embryos using lightsheet microscopy is very successful due in large part to the high transparency and low scattering in those particular embryos.^1,2 Achieving a greater depth of focus can be achieved using more sophisticated illumination geometries, such as Airy or Bessel beam illumination in lightsheet microscopy.³ However, imaging of samples such as chick embyros, which constitutes a more physiologically relevant model for human embryonic development, is greatly hindered by greater levels of scattering in the embryonic tissue.⁴

The standard approach to modeling light propagation in scattering media, from interstellar space to human tissue, is to consider photon transport via a Monte Carlo radiation transfer (MCRT) algorithm.^5–7 However, MCRT does not usually preserve the phase state of a photon with each scattering event and as such does not allow for the modeling of coherent scattering or interference effects. Several authors modified the MCRT algorithm in an attempt to model these effects.^8–11 However, most of these methods either inaccurately model Gaussian beams, are complex to implement, or carry a heavier computational burden. Most three-dimensional (3D) optical imaging approaches use a shaped illumination field, such as a Gaussian beam, a lightsheet Bessel beam, or an Airy beam. All of these require a spatially coherent illumination source, with the shape of the illumination arising from the interference of, for example, spherical, conical, or cubic phase distributions. We present a phase encoded MCRT ($\phi \mathrm{MC}$) that preserves the phase of the photons as they scatter in a model tissue, allowing the modeling of arbitrary beam types. This opens up the possibility to make direct and “fair” comparisons between different illumination geometries that are otherwise very challenging to directly compare experimentally.

The use of simple high quality Gaussian beam illumination of a sample, either as a scanning microscopy or in a digitally scanned lightsheet microscopy geometry, has been successful in imaging in a range of different samples.^2,4,12 However, scattering eventually limits the depth to which contrast and resolution can be maintained in an image. Hence, a different approach is needed for these more challenging tissues. Many different approaches to overcoming this challenge have been proposed and demonstrated. These include the direct physical manipulation of the scattering properties via optical clearing of the tissue,^13,14 the precompensation of the illumination wavefront to obtain “aberration free” illumination,¹⁵ and the use of beam types with different wavefronts such as the simple spherical wavefront of a Gaussian beam or more complex shaped light beams such as Airy or Bessel.^3,16 Often these approaches are combined, such as the use of deconvolution in tandem with Bessel beam illumination.¹⁷ What is not always clear is which of these approaches is the most likely to succeed for a sample, or even if these approaches necessarily provide more information on a sample than a simple approach performed well, especially as they inevitably present much greater optical or computational demands.

After first validating on a Young’s slit arrangement and the propagation of Gaussian beam in free space, we apply it to comparing the depth penetration of Bessel and Gaussian beams in homogeneously scattering model tissue. Bessel beams have what are often referred to as “self-healing” properties.¹⁸ This property arises because the beam, as observed subsequent to an obstruction, is formed by off axis rays that might not have been obstructed themselves. This self-healing is particularly powerful for maintaining optical intensity in the illumination in the presence of not only localized absorption but also localized scattering. However, in a homogeneously scattering medium, in which the complete wavefront of the illuminating beam is affected by the scattering, it is less clear whether or not a conical (Bessel) or spherical (Gaussian) wavefront will be more robust.

Because the energy is distributed so differently in a Bessel beam (between the central maximum and rings) than it is in a diffraction-limited Gaussian beam (all energy in a single spot), it is hard to make direct comparisons between the two. For example, should one compare equal total power in the beam, equal intensity in the Bessel beam’s central spot versus at the focus of the Gaussian beam, or something more physiological and sample specific such as equivalent phototoxicity? This range of approaches makes an experimental comparison very challenging. With such comparisons are made much simpler and subsequently fairer.

2.

Methods

3.

Validation

To verify the additions to the MCRT method, we model the double-slit experiment, diffraction through a square aperture in the near and far fields, and show that Gaussian beams and Bessel beams are correctly modeled through free space.

3.1.

Double-Slit Experiment

The first test of $\phi \mathrm{MC}$ is to compare our algorithm to the double-slit experiment. In the double-slit experiment, a monochromatic plane wave of light is incident on two slits distance apart $d$, width $2b$, and an interference pattern is observed on a screen a distance $L$ away from the slits. The experiment is usually carried out with the detector screen in the far field, the so-called Fraunhofer regime. The intensity pattern on the detector screen is as in Eq. (4):²³

Eq. (4)

$$I(x)\propto {\cos }^2\left(\frac{kdx}{2\sqrt{L^2+x^2}}\right)\sin c^2\left(\frac{kbx}{\sqrt{L^2+x^2}}\right),$$

where $k$ is the wavevector, $k=\frac{2\pi }{\lambda }$, and $x$ is the horizontal position on the detector screen. The simulation was carried out for a wavelength, $\lambda$, of 488 nm, a slit width of $10\lambda$, slit separation of 80 $\lambda$, and the detector screen positioned 10,000 $\lambda$ away from the slits. Using the Huygens–Fresnel principle, each slit is a source of Huygens wavelets. The detector screen has dimensions, $1{\mathrm{mm}}^2$ and there are ${2051}^2$ bins, giving a bin an effective size: $\approx 488\mathrm{nm}$ or $\approx \lambda$. The initial position of the photon packets is sampled uniformly from the slit area, after randomly choosing one of the slits to emit from. A random direction is then chosen to ensure that the packets will hit the detector screen. This is achieved by uniformly sampling a position on the detector screen and calculating the direction vector. The simulation was run with ${10}^9$ packets and the code was parallelized with MPI,²⁴ which took $\approx 10\min$ to run on an 8 core Intel Xeon machine. This gave an accurate match to the theoretical expression, as shown in Fig. 1.

Fig. 1

Comparison of theory and simulation for the double-slit experiment. (a) A slice through the computed image and the expected profile from theory. For clarity only every fifth MCRT data point is plotted. (b) The computed image.

3.2.

Diffraction by Square Aperture

We also validate the ability of $\phi \mathrm{MC}$ to simulate wave phenomena in the near field, the so-called Fresnel regime. To test this, we compare results from the simulation against theoretical predictions for diffraction through a square aperture, for both the Fresnel (using the paraxial Fresnel approximation) and Fraunhofer regimes. The Fresnel regime is defined when the Fresnel number, Eq. (5), is greater than 1.0, and the Fraunhofer occurs when the Fresnel number is $<1.0$. We define the Fresnel number as in Eq. (5),¹⁰ where $l$ is the slit width, $\lambda$ is the wavelength, and $r_0$ is the distance to the detector screen. As with Young’s slit experiment, the slit is a source of Huygens wavelets. However, we vary the distance to the detector from the square slit such that the Fresnel number changes to validate in the both diffraction regimes:

Eq. (5)

$$F=l\sqrt{\frac{2}{\lambda r_0}}.$$

We used a wavelength of 351 nm and a slit width/height of 100 nm. For the Fresnel regime, we used 300 bins covering $3600\mu \mathrm{m}$, and for the Fraunhofer, we used 100 bins covering $6000\mu \mathrm{m}$. More bins and photons were required to resolve the complex pattern for higher Fresnel number patterns. Figure 2 shows that matches the theory for both the Fresnel and Fraunhofer regimes.

Fig. 2

Comparison of theory and simulation for diffraction through a square aperture in the Fresnel and Fraunhofer regimes, for a variety of Fresnel numbers.

3.3.

Gaussian Beams

To show that $\phi \mathrm{MC}$ can model Gaussian beams, we first model a plano-convex lens and focus a Gaussian beam to its focal point and beyond. The plano-convex lens is based upon ThorLabs LA4249 UV-fused silica lens,²⁵ which has a radius of 5 mm, thickness of 2.2 mm, and radius of curvature 4.6 mm. A Gaussian beam with $\frac{1}{e^2}$ width²³ 0.5 mm and wavelength 488 nm is incident on the lens. This is propagated through the lens using Snell’s law. When the packet reaches the far surface of the lens, the Huygens–Fresnel principle is used to sample the packet onto the medium’s surface uniformly. The packet is then transported through free space with the usual MCRT method. Figure 3 shows the comparison of theory and in-silico experiment, with excellent agreement between the two.

Fig. 3

Comparison of Gausssian beam simulation and theory. (a) The Beam waist as a function of distance from its minimum. (b) The width of the Gaussian beam at its minimum.

4.

Bessel Beams

With a little adjustment to the simulation setup, we can also model $n$’th-order Bessel beams. We change the lens to that of an axicon lens for the 0’th-order beam and to an axicon with a helical delay for a higher order Bessel beams. The axicon used here has an opening angle of 5 deg ($\alpha$) and radius 12.7 mm. A Gaussian profile beam is incident on the axicon with $1/e^2$ width 1 mm. The beam is propagated through free space as described before to a detector screen 10 mm from the tip of the axicon lens. The detector screen has a size of $40\mu \mathrm{m}\times 40\mu \mathrm{m}$ with a bin resolution of $1\mu \mathrm{m}$. $8\times {10}^{10}$ photon packets were simulated taking $\sim 1\mathrm{h}$ on an 8-core Intel Xeon 3.5 Ghz machine. Equation (6) gives the equation of a theoretical Bessel beam at a depth $z\mathit{max}$ normalized. $J_0$ is the zeroth-order Bessel function, $r$ is the radial distance, and $k_r$ is the radial wavevector. Equation (6) is plotted against the simulation data, with the simulation normalized to the maximum intensity of the image generated. Figure 4 shows this comparison:

Fig. 4

Comparison of theoretical and MCRT simulation of a Bessel beams, with intensity normalized. The results from $\phi \mathrm{MC}$ show good agreement with the theory.

Fig. 5

Bessel beam in the far field. Bessel beams in the far field becomes a ring beam. Image shows a slice of intensity through the medium.

4.1.

Comparing Simulation to Experiment

To ensure our algorithm gets the correct results in turbid media, an experiment was carried in an intralipid scattering agent. The experiments consisted of imaging a Bessel beam after its propagation through an Intralipid medium of varying turbidity. The laser used to create the Bessel beam has a Gaussian profile, and a wavelength of 488 nm, which is incident on an axicon lens with an opening angle 5 deg. The laser beam had a $\frac{1}{e^2}$ diameter of 2 mm at the surface of the lens. The Bessel beam was propagated through 10 cm of air before encountering the Intralipid media. The Intralipid was kept in a cuvette of side 2 mm. The cuvette was filled with $500\mu \mathrm{L}$ of water, and various volumes of the scattering agent are added. Figure 6 shows the experimental set-up.

Fig. 6

Experimental setup for propagating a Bessel beam through a cuvette filled with varying concentrations of Intralipid 20%. Bessel beam is imaged by an $20\times$ objective lens and a Grasshopper 3 camera.

The scattering agent used is Intralipid 20% (Sigma-Aldrich), which is diluted as shown in Table 1. Figure 7 shows the optical properties of Intralipid 20%. Dilutions of Intralipid are kept below 2% scattering particle concentration, so the scattering exhibited by the Intralipid is in the independent scattering regime. The independent scattering regime is where $g$ (the anisotropy factor, which is a measure of scattering direction after a scattering event) is dependent on the size, shape, and material properties of the scattering particle, and the material properties of the bulk material but not the number of scattering particles.^27,28 This allows the linear scaling of the optical properties by concentration.^27,29,30 Images of the Bessel beam as it emerges from the cuvette are taken for comparison with our algorithm.

Table 1

Intralipid solutions used for experiment, see also Fig. 7.

Volume (μL)		Intralipid concentration		Optical properties
Intralipid	H2O	Volume (%)	Scattering particle (%)	Scattering coefficient (m−1)
0	500	0.00	0.00	0.00
2	500	0.39841	0.0908	557.14
4	500	0.79365	0.1816	1114.28
6	500	1.18577	0.2724	1671.42
8	500	1.57480	0.3632	2228.56

Fig. 7

Scattering properties of 20% Intralipid.²⁶

To model within $\phi \mathrm{MC}$, we simplify the experimental setup considerably. The simulation models the propagation of photon packets through the axicon to its conical surface. On the conical surface, the Huygens–Fresnel principle is invoked, and the packet is sampled onto the surface of the medium (cuvette). The sampling of the photon onto the surface of the medium speeds the algorithm up, as it does not need to simulate the photons that would “miss” the medium. From there, the usual MCRT method propagates the packet through the medium while tracking its phase, and scattering the packet until it leaves the medium. If the packet leaves the medium to any side other than the far side of the cuvette (e.g., any side of the cuvette not facing the objective lens), then it is discarded. If the packet leaves the medium on the objective lens facing side, then the packet is recorded by its phase onto an area element. For each Intralipid concentration $6.4\times {10}^{10}$ photons are run over 64 cores, taking $\sim 3\mathrm{h}$ for the $12\text{-}\mu \mathrm{L}$ Intralipid volume. Once all the packets have been run, the phase is converted into intensity, as in Eq. (3), but in two-dimensional.

Figures 8 and 9 show the results from the experiment and simulation. To allow for the comparison between the experimental and in-silico data, we normalized each image to its brightest pixel. The simulations show good agreement with the experimental data within experimental and simulation uncertainty.

Fig. 8

Comparison of experimental and simulation data for propagation of a Bessel beam produced by an axicon, through mediums of various turbidity. (a)–(g) The data from $\phi MC$, and (h)–(n) the experimental data. Volumes along the top are the volume of Intralipid in each solution as in Table 1. All images are cropped so they are the same size and normalized to the maximum value in each image.

Fig. 9

Line graph plots of slices taken through the generated and experimental images as shown in Fig. 8.

5.

Comparison of Bessel and Gaussian Beams

As Bessel and Gaussian beams are radically different from one another, it is hard to directly compare the two beams. Gaussian beams carry all their power in the “central core” of the beam, whereas in a Bessel beam, it carries the same amount of power in each ring. Bessel beams also have a much larger depth of focus than Gaussian beams. We attempt to give a fair comparison of the two beams, to predict which beam performs better in a heavily scattering medium using $\phi \mathrm{MC}$. Bessel beams are expected to perform better than Gaussian beams, due to their self-healing properties and nondiffractive core, this section aims to quantify how this property may or may not help penetration through a highly scattering medium.

As mentioned, Bessel beams and Gaussian beams are not alike, so to ensure a fair comparison the Bessel beams central core width [Eq. (7)] is set to that of the Gaussian beam’s waist:

The first simulation comparisons carried out between the Bessel and Gaussian beams is to use the same power to generate both beams. The beams are then propagated through mediums of varying degrees of Intralipid solution. Volumes of 0.0, 26, 52, 78, and $104\mu \mathrm{L}$ are used of Intralipid in $500\mu \mathrm{L}$ of water. The medium has a volume of $0.1\mathrm{mm}\times 0.1\mathrm{mm}\times 0.2\mathrm{mm}$, and voxel resolution of $1\mu \mathrm{m}$. For both beams, a wavelength of 488 nm and a power of 1 mW is used. One hundred million packets are simulated for each simulation. The results of this are shown in Figs. 10 and 11.

Fig. 10

First comparison of Bessel and Gaussian beams with equal power used to generate both beams. Plots taken at the Gaussian beams focus. The maxima at the sides of the Gaussian beam in the $0.0\mu \mathrm{L}$ plot are due to simulation effects, mainly the small size of the medium not allowing photons from further off the optical axis to interfere destructively.

Fig. 11

First comparison of Bessel and Gaussian beams, with equal power used to generate both beams. Plots taken at the bottom of the simulated medium. Medium has a 2 mm thickness.

The results show that for the same power, Gaussian beams propagate deeper into the medium compared with Bessel beams. This is to be expected as in a Gaussian beam all the power is in its “central core,” whereas the power is evenly distributed between all the Bessel beam’s rings. Therefore, for a second comparison, the power given to the Bessel beam is such that the central core maximum matches that of the Gaussian beam at its focus for the case where there is no scattering. To achieve this, the Bessel beam was given $\sim 15\times$ the power given to the Gaussian beam. The results of this comparison are illustrated in Fig. 12.

Fig. 12

Second comparison of Bessel and Gaussian beams for the case where the power given to each beam yields the same maximum at the Gaussian beams focus. These plots are taken from the Gaussian beams focus. Medium has a 2-mm thickness.

These results show as expected that the Bessel beam now performs comparably with the Gaussian beam in lower scattering media, with a drop off in performance in the higher scattering media.

6.

Discussion

For equal power beams in the previous section, Gaussian beams perform “better” in the highly scattering media, though this is expected as the power in a Bessel beam is spread evenly over its rings. Thus, the power in the central lobe of a Bessel beam is much less than that of the Gaussian beam.

To give a slightly fairer comparison of intensity in the central lobes of the beams, the Bessel beam was given $15\times$ more power than the Gaussian beam. This allows a better comparison between the Gaussian beam’s core and the Bessel beam’s core and gives a more comparable intensity between the beam’s at the location of the Gaussian beams focus. In this case, the Bessel beam appears to perform better in a highly scattering medium, as shown in Fig. 13. The Bessel beam shows comparable intensity with the Gaussian beam in the first three mediums, though the Gaussian beam out performs the Bessel beam in the higher scattering media. It would appear that the Bessel beams self-healing property does not help a Bessel beam propagate through a highly scattering medium. As photons propagate through the medium, they interfere with one another constructively and destructively to form a Bessel beam. However, if enough photons are scattered, then the Bessel beam becomes degraded and thus no longer is a Bessel beam, as these photons are no longer coherent with the rest of the beam, so they act as a negative factor in the beams formation. Another reason that the “self-healing” property of the Bessel beam does not “save” the beam from scattering is that the “self-healing” is not self-healing. The self-healing in reality is just photons from further off the optical axis forming the Bessel beam further down the optical axis, e.g., the photons that are impeded by the blockage are stopped, but the photons that are not impeded form a Bessel beam as expected. If you placed a blockage in front of the Bessel beam larger than the width of the input beam, the Bessel beam would not form at all.

Fig. 13

Comparisons of unequal powered beams at the bottom of scattering medium. Medium has a 2-mm thickness.

Bessel beams do have their positives, their self-healing property does help “reform” the beam past small blockages, and their depth of field is superior to an equivalent Gaussian beam, as their central core is “nondiffractive.”

This technique is fairly computationally heavy due to the stochastic nature of the simulations. However, with the advent of fast graphical processing unit, this computational cost going forward is not the barrier it may have once posed.

7.

Conclusion

This work presents an adaption to the Monte Carlo radiative transfer algorithm to allow the simulation of the wave properties of light. We showed that MCRT can be easily modified to allow the simulation of interference of light. These additions to the MCRT algorithm allowed a fair comparison between Gaussian and Bessel beams, something that is not easily achievable in a lab setting. We showed that despite Bessel beams self-healing property, we found that it does not outperform a comparable Gaussian beam when it comes to imaging at depth. However, the Bessel beams’ depth of focus is still vastly better than that of the Gaussian beam.