3.1.

Synthetic Models

A synthetic model has been created with a geometry and optical properties mimicking those of skin. The model has cells, subcellular components such as mitochondria, melanin, and a nucleus, and a dermal–epidermal (DE) junction and is discussed in detail elsewhere.^2,7,12 An index of refraction ($n$) map can be found in Fig. 1(a). The indices of refraction were chosen to be representative for near-infrared light. The size and shape of the cells vary with depth. Cells near the DE junction have a circular shape with a radius of 10 μm and flatten out as they approach the skin surface, with their major axes linearly approaching 20 and minor axes approaching 5 μm. The mitochondria and melanin are randomly distributed within each cell, and the cell spacing is also randomly varied to ensure there is no periodicity introduced into the signal. Absorption is assumed to be negligible compared to scattering. The model is based on the work of Simon and DiMarzio⁷ who used indices of refraction collected from the literature by Dunn and Richards-Kortum.¹³

Fig. 1

(a), Geometry of a synthetic model illustrating cells, subcellular components, and a dermal–epidermal junction. The color bar indicates the index of refraction for each component. (b), Measurement geometry with the location of the ultrasound (US) beam is illustrated with a rectified beam in green. The source location is illustrated with a solid red line and the detector with a dashed yellow line.

The tissue model is illuminated with a plane wave. A plane wave with a wavelength of 632 nm in free space illuminates the surface of the simulated skin and is detected in a transmissive geometry as shown in Fig. 1(b). The input optical beam takes 4.17 fs to reach full intensity and has a beam width of 80  μm. The simulation space is divided into pixels of $63.2\times 63.2\mathrm{nm}$ and has a temporal step size of 0.12 fs. A US beam is utilized to modulate the pressure in the tissue; hence the index of refraction. The US wave is modeled by a 2-D Gaussian pressure field. The Gaussian wave is similar to the spherical Gaussian wave used in optics,¹⁴ except it is constant in the third dimension. The US wave is focused to a spot size of 5.9 μm, centered at 60 μm in depth and 40 μm from the left along the transverse axis.

A second model was created to simulate a highly scattering medium. The model simulates titanium oxide ($n=1.48$) suspended in gel ($n=1.33$), mimicking phantoms used recently in previous scattering experiments.¹⁵ The phantom is shown in Fig. 2. This is a 2-D simulation, so the scatterers are actually cylindrical and, for our numerical calculations, we assume they extend 1 m along the third dimension. The phantom measures $80\times 100\mu \text{m}$ and contains approximately 6400 randomly distributed cylinders. The cylinders have a diameter of 532 nm but are allowed to overlap, so they can make complicated shapes. Even though the particles are randomly distributed, their positions are fixed for all simulations. The phantom is embedded in a $100\times 100\mu \text{m}$, index-matched simulation space. The scattering coefficient of the medium is calculated to be $3165{\mathrm{cm}}^{-1}$ for a wavelength of 532 nm by computing the extinction of a plane wave in the simulation. While the scattering coefficient is an order of magnitude larger than expected for most tissue, it was chosen so that diffusion effects may be observed within a reasonable computational area. This phantom does not contain absorbers.

Fig. 2

Geometry of randomly distributed cylinders to create a random scattering medium. The US beam location is illustrated in green.

The model is discretized into pixels that measure $53.2\times 53.2\mathrm{nm}$ and a temporal step size of 0.11 fs. The US modulation is modeled as a cylindrical beam with a diameter of 3.7 μm centered at 50 and 30 μm along depth and transverse axes, respectively. This model would be difficult to implement physically, but it provides insight into the microscale effects of modulation. A plane wave with a 532-nm wavelength is turned on and plateaus at full intensity after 4.1 fs. The beam extends 80 μm along the transverse axis.

3.2.

Simulated Optical Detection

A simulated detector collects the optical signal for each time-step during an FDTD run. The analytical signal is computed for each of the five time-varying detector signals using Matlab’s Hilbert transform function. The analytical signal is then mixed with a reference optical beam, $\exp (j2\pi f_{\mathrm{opt}}t)$, signal to mix it down to baseband. The resulting sidebands are at the frequencies $\pm f_{\mathrm{us}}$. Each sideband may be moved to DC by further mixing the signal with a reference US signal.

The sidebands are conjugated to create a phase conjugate sideband (PCS) wave that will ideally deliver its power to the US beam waist. As mentioned earlier, phase conjugation will force the optical signal to retrace its path. Because we have isolated the light modulated by the US beam, it is expected to travel back to the US waist. It has been shown that phase conjugation can be achieved with only 0.02% of the incident optical power,¹⁶ and in our simulations, most of the incident optical power is collected (except for light lost at the boundary) because absorption is ignored. The PCS signal may be assigned an arbitrary amplitude that can be chosen to increase power delivery. The unmodulated carrier signal is also conjugated to create a phase conjugate carrier (PCC) wave. We expect this signal to retrace its path through the medium and converge back to the source without focusing in the medium. The amount of energy focused by the PCS and PCC delivery waves is compared to evaluate the performance of optical phase conjugation using a US guidestar.

4.1.

Tissue Phantom

An example simulation of a 632-nm wavelength optical probe beam propagating through tissue modulated by a US beam is shown in Fig. 3. The tissue model appears in gray, the source and detector positions are indicated by a solid red and dashed yellow line, respectively, and the rectified US and optical beams are illustrated with green and red, respectively. Because the signals are rectified, the visualization shows the waves having half their actual wavelengths.

Fig. 3

An illustration of the probe beam field in tissue. The tissue model appears in gray, the source and detector positions are indicated by solid red and dashed yellow lines, respectively, and the rectified US and optical beams are illustrated green and red, respectively. Because the signals are rectified, they appear to be oscillating at half their actual wavelengths.

It is interesting to note that the optical beam forms into channels within the medium, creating localized bands of increased energy. These optical channels appear to be coherent over long distances and can be as little as one wavelength in width. They are also primarily directed radially away from the source. The positions of the optical channels depend on the wavelength of light and geometry of the tissue. The presence of optical channels has been noted previously in FDTD simulations,¹⁷ although it was not seen in a different type of scattering medium.⁵ These models are all in 2-D, which raises the question of weather the channels exist in 3-D experiments. One published 3-D simulation seems to show similar channels.¹⁸

The relative irradiance is calculated by normalizing the detected power by the source power. An example result is shown in Fig. 4(a) for the steady-state case (US beam is off). The measured phase is unwrapped using Matlab’s unwrap function and is shown in Fig. 4(b). It is interesting to note that the phase delay is not uniform after being transmitted through the tissue phantom, even though the initial wave’s phase was. This is because waves scattered by several small particles will add coherently, like the Huygens wavelets introduced in optics courses. If the incident wave is a plane wave, the result will be, on average, a plane wave. The curvature of the wavefront from one particle would be canceled by that from others. However, in the present case, the incident wave has finite extent (80 μm in a computational area 100 μm in width). The curvature of the wavefronts of the illuminated scatterers are not completely canceled, and that gives rise to the measured field’s curvature.

Fig. 4

(a), Irradiance of the measured probe beam (with the US switched off) after propagating through the tissue model normalized by the source power. (b), The phase delay of the probe beam after propagating through the tissue. The data (blue) is fitted with a second-order polynomial (green) with a radius of curvature of 98.6 μm. (c), Difference between the detected phase and the second-order fit.

A second-order polynomial is fitted to the data utilizing least squares,¹⁹ which represents a spherical wavefront having a source near the origin and has a radius of curvature of 98.6 μm. As one would expect, the wavefront differs from this quadratic curvature due to scattering. The difference between the measured phase and polynomial is shown in Fig. 4(b).

The sidebands, discussed in Sec. 3, are then isolated by demodulating the measured signals, as discussed in Sec. 3.2. The relative demodulated probe power is shown in Fig. 5(a), and its phase is shown in Fig. 5(b) in blue. The radius of curvature is found to be 66.9 μm. The smaller radius of curvature indicates that the wavefront is centered in the medium closer to the US waist. The difference between the measured phase and the second-order polynomial varies by $>18$ radians and is illustrated in Fig. 5(c).

Fig. 5

(a), Irradiance of the demodulated optical signal normalized by the source power. (b), The phase of the demodulated optical signal. The phase (blue) is fitted with a second-order polynomial (green) with a radius of curvature of 66.9 μm. (c), The difference between the demodulated phase and the second-order fit.

Delivery waves are created from the phase conjugate carrier and sidebands and simulated in the synthetic tissue. The results [Fig. 6(a)] show that the PCS delivery wave exhibits strong focusing at the US guidestar position. The energy is highly localized in the horizontal, and less so in the vertical axis. These results agree with experimental results.⁴ The PCC [Fig. 6(b)] retraces its path through the medium and converges back to the probe beam source. This is because the signal is the phase conjugate of the probe beam. Phase conjugation forces it to retrace its path through the scattering medium.

Fig. 6

Phase conjugate carrier and sidebands in the synthetic tissue model. The model indices of refraction are illustrated in gray, the US beam is in green, and the optical field is in red.

To evaluate the performance of the guidestar, a measure of energy concentration is needed. We propose a metric that is independent of the source power. The time-averaged energy within the medium is normalized by the optical input power. The time-averaged energy is computed as

Fig. 7

The energy-to-input-power ratio for the focused (green solid) and unfocused (blue dashed) delivery waves on a log-log scale. The values are plotted for growing distances from the center of the US beam (shown as R).

To further investigate the focusing of the PCS beam, the time-averaged Poynting vector is calculated and illustrated in Fig. 8 for the probe beam [Fig. 8(a)], PCS [Fig. 8(b)], and PCC [Fig. 8(c)]. The direction of the Poynting vector is illustrated with color, where red indicates the vector is directed down along the depth axis and cyan indicates a vector in the opposite direction. The magnitude of the Poynting vector (irradiance) is utilized as an intensity envelope for the figure. Pixels with a higher intensity correspond to higher irradiance values. The pixels containing a magnitude below a threshold are set to white. The optical channeling phenomenon is visible. These graphs illustrate that the power is traveling along these channels with little divergence.

Fig. 8

The angle of the Poynting vector is utilized to create an image where cyan suggests the vector is directed down along the depth axis and red suggests the opposite direction. The magnitude is utilized as an intensity envelope for the figure. Pixels with a higher intensity correspond to higher magnitude values.

4.2.

Scattering Phantom

The same process is repeated for the titanium oxide phantom containing random scatterers. A 532-nm probe beam interrogates the medium as shown in Fig. 2. As before, the field is measured at the detector (position is indicated by the yellow line) for four US phases. The received signal is demodulated to isolate the sidebands, and Fig. 9 shows the resulting amplitude [Fig. 9(a)] and phase [Fig. 9(b)]. A second-order polynomial is also fitted to the unwrapped phase (green). The difference between actual phase and the fitted second-order function is shown in Fig. 9(c) and is larger than the previous example owing to the high scattering.

Fig. 9

(a), Irradiance of the US-demodulated optical signal normalized by the source power for an example titanium oxide phantom. (b), The phase of the demodulated signal. Note the phase (blue) is fitted with a second-order polynomial (green) with a radius of curvature of 57.5 μm. (c), Difference between the demodulated phase and the second-order fit.

It is well known that energy can be directed to some extent by focusing into tissue, as is done in confocal microscopy. To prove that the phase conjugation results here could not be achieved just by focusing, a delivery signal is created with a phase given by the second-order curve shown in Fig. 9(b), and the amplitude is unchanged. If we removed the scatterers, this beam, the phase conjugate fit (PCF), should focus to a point. The PCF wavefront is simulated as a delivery wave along with the phase conjugate sideband and carrier waves (Fig. 10). Recall, the phase conjugate carrier is the phase conjugate of the unmodulated probe beam. That signal [Fig. 10(a)] travels back to the source, as expected. The phase conjugate fit [Fig. 10(c)] exhibits some focusing in the medium, but the spot is very broad due to the scattering. However, the phase conjugate sideband [Fig. 10(b)] focuses to a tight spot within the US beam width despite the high scattering. This comparison demonstrates that the high irradiance of light delivery is the result of phase conjugation, and not merely focusing. It is also interesting to note that the optical channeling effect is still present even after the signal has traveled a mean free path in the medium.

Fig. 10

The phase conjugate carrier (a), phase conjugate sideband (b), and phase conjugate fit (c) delivery waves in an example titanium oxide phantom. The US beam is illustrated in green, the medium of random scatterers in gray, and the electric field in red.

The Poynting vectors are computed in in a similar manner as above for the phase conjugate carrier and sideband, and fitted fields are shown in Fig. 11. The results show that the PCF and PCC signals exhibit little focusing. The PCS exhibits a converging energy density at the beam waist that quickly diverges beyond it. Figure 10 illustrates the resulting EP metric for each signal on a log-log plot. The EP ratio is lowest for the PCC (green dotted line), but it is interesting to note that the phase conjugate fit (blue dashed line) has about the same ratio as the PCS signal near the center of the US beam. This appears to contradict the field and Poynting vector distributions in the medium. However, the ultrasound beam radius is 3.72 μm (indicated by the vertical cyan line in Fig. 12); the energy delivered to the US beam is 2.6 times higher for the PCS than the phase conjugate fit signal. This suggests that the heterogeneity of the optical channels is skewing the result for small distances. As discussed earlier, phase conjugation can achieve suboptical wavelength resolution through the diffusive lens effect.³ However, the US guidestar cannot select the location for this focusing to an accuracy less than the US resolution. Therefore, we shall use the radius of 3.72 μm to characterize the enhancement. Note that the EP ratio (energy for a given input power) for optical phase conjugation with a US guidestar drops as the opacity of the medium increases. This agrees with previously reported experimental results¹⁶ as well as other results obtained from other computational models.⁵ It is important to remember that these simulations allow for variation in two dimensions. Although we do not have 3-D simulation data, geometric arguments suggest that the enhancement would be squared in three dimensions. In our example, the enhancement by a factor of 2.7 in 2-D means the focused energy will be $>7$ times greater than the phase conjugate fit signal in three dimensions.

Fig. 11

Poynting vectors for the phase conjugate sideband, carrier, and fit sources in the medium containing randomly positioned scatterers.

Fig. 12

Log-log plot of the energy-to-input-power ratio of the phase conjugate fit (blue dashed), phase conjugate carrier (green hashed), and phase conjugate sideband (red solid) waves.