1. Introduction

Ultrasound optical tomography (UOT) is a nonintrusive medical imaging modality which since its emergence in the 1990s has seen steady development [1,2]. In this imaging scheme, tissue is illuminated by light while simultaneously being illuminated by ultrasound (US). In the interaction region between the acoustic and optical fields, part of the light is shifted in frequency and "tagged" due to the acousto-optic effect. As it is the overlap of the two fields which produces the signal, the high scattering of light allows for the much less scattered US to create a source of tagged light from the untagged light, which due to scattering reaches a large tissue volume. By steering the direction of the US pulses within this volume and timing optical pulses to probe different depths, it is possible to create images of the optical properties of the medium by recording the tagged signal strength for the varying US pulse positions.

This tagged light must, however, first be differentiated from the much stronger background of untagged light which has not interacted with the ultrasound. This separation is possible using different techniques, such as using: heterodyne signals of individual light speckles [3], contrasts in the output speckle pattern [3], photorefractive crystals imprinted with a reference beam [4] or, more recently, laser feedback signals from tagged light coupling into a laser cavity [5]. The tagged light intensity can also be extracted using spectral holeburning filters or slow light filters (SLFs) prepared in rare earth doped crystals [6–10], and the use of SLFs is the selected experimental filtration approach in this paper. If UOT can be realized to the optimized extent discussed in Ref. [8], it could allow for fast, non-intrusive medical imaging of optically deducible quantities, such as blood oxygenation, in human tissues at depths over 5cm in reflection mode or 10cm in transmission mode.

A key aspect for reaching this goal is knowing how strong the acousto-optic tagging effect is for the used acoustic scheme. This has been extensively studied for many schemes, such as focused continuous wave ultrasound [11–16] and pulsed wave ultrasound [17,18]. Even schemes using planar acoustic waves have been investigated, though use of a pulsed ultrasound has the advantage of yielding a higher spatial resolution [19]. However, in common for these previous studies has been the use of low acoustic pressures, possibly due to the issues arising from the necessary nonlinear description of acoustics for higher pressures. Although at least one study using high pressure ultrasound has been performed [20], this did not address the nonlinear frequency evolution experienced by the acoustic pulses at these pressures. As such, a model which can describe the tagging process by an arbitrary US field in pressure regimes up to and just below the medical safety limits is currently poorly covered by the literature. Furthermore, such a broad model would allow for prediction of the acousto-optical tagging (i.e. UOT signal) strength of any potential US pulse used in UOT.

In this paper we explore such a model by iterating on the work performed in Ref. [16], showing that the model presented there is also valid for high pressure pulsed ultrasound under certain limitations. Using the software made available by the authors of that paper as inspiration, we from the ground up re-implement the model in MATLAB in order to allow for simulations using arbitrary acoustic fields while also significantly increasing the computational speed. We continue to present experimental measurements of pulsed acousto-optic tagging efficiencies in tissue phantoms together with detailed measurements of the used acoustic pulses. The experimentally determined tagging efficiencies are then compared with simulations using the measured acoustic pulses. Finally, we discuss the US frequency’s and US pulse irregularities’ impact on the tagging process and their implications on UOT.

2. Theory

In the theory presented by Raman and Nath, as an electro-magnetic field propagates in a clear and transparent medium together with an acoustic wave, it accumulates a phase change varying in time with the same frequency as the acoustic field due to the acoustically induced change in refractive index [21]. This harmonic phase variation generates sideband frequencies $n f_{\textrm {US}}$ away from the initial optical carrier frequency $f_{\textrm {o}}$ where $n$ is the order of the sideband and $f_{\textrm {US}}$ the fundamental ultrasonic frequency.

As described in Refs. [12,13], the same frequency shift occurs for an electro-magnetic field propagating in an optically scattering medium with some key differences. In addition to the phase change variation caused by the change in refractive index, the optical field accumulates a phase shift from the movement of the scattering particles. This can further be divided into two effects: a variation in optical path length due to movement of the scattering centers and a change in the scattering rate ($\mu _s$) of the medium [18]. As discussed in Ref. [14], the latter of these effects is much weaker than the previous two and is neglected from here on. A useful quantity for describing the amount of frequency shifted light in each order is the tagged fraction $\eta _n$, defined as

In the following subsections, we discuss how a high pressure ultrasonic pulse propagates and can be described inside a liquid. We then continue with how the accumulated phase of a wave packet of light is affected due to the existence of such an acoustic pulse.

2.1 Acoustic wave propagation and particle displacement

The relation between a time dependent pressure field $P(\mathbf {r},t)$ and the particle displacement field $\mathbf {A}(\mathbf {r},t)$ can be described using fluid dynamics, assuming that the scatterers follow the movement of the surrounding fluid. We describe the fluid dynamics using Euler’s equations for an incompressible and non-dissipative fluid [22]

Taylor expansion of Eq. (5) yields that for propagation distances $|y| \ll D$, the transfer of power from the initial excitation frequency to its harmonics is negligible. This is the same linearization criteria fulfilled by low pressure excitations with the difference that in those cases, $D$ is much larger than the typical regions of interest due to small $\Gamma _0$.

We now turn back to Eq. (7) and treat that as an initial pressure excitation. We then approximate that each of the $N$ known harmonics behave as the $n=1$ harmonic in Eq. (5), i.e. they all display a negligible power transfer to even higher harmonics for $|y-y_0| \ll D$. We can therefore approximate $a_n$ as constant for $|y-y_0| \ll D$ which in turn allows us to describe our pulses as

2.2 Light phase accumulation and frequency shift in scattering tissue

We describe the light traveling inside the scattering media from our source to the detector as $M$ coherent wave packets, indexed with $m$, traveling along $J_m$ straight paths between in total $J_m-1$ scattering events. Given a homogeneous medium and assuming that the power is initially divided equally over all paths, we describe the total optical power spectrum at a detector, $I_{\textrm {det}}$, as

For a wave packet that has traveled through some scattering medium from a point source at $\mathbf {r}_{m,j=0}$, the accumulated phase only depends on the accumulated optical path length and takes the form

3. Experiments

Two experimental setups were used to measure the ultrasonic pressure field and the acousto-optic tagging efficiency at different ultrasound output powers. The ultrasound field was measured using a hydrophone fixed on a triaxial translation stage immersed in water. Using the hydrophone, five ultrasonic pulses of different lengths with either 1.6 or 3.5 MHz center frequencies were measured in the transverse plane of the transducer’s focii at $y=y_0=3.5$ cm. The ultrasound source was then moved to the optical setup where, using the same five ultrasound pulses, the tagging efficiency was measured in solid intralipid phantoms with the ultrasound again focused at a depth of 3.5 cm. The tagged signals were measured using an SLF prepared in a Pr:YSO crystal. Simplified versions of the setups are depicted in Fig. 1.

 

Fig. 1. Simplified illustration of the setups for the (a) acoustic field strength and (b) tagging efficiency measurements. i: motorized translation stage, ii: hydrophone probe iii: US transducer. The US field was measured in the plane of the US focus (fixed y position). iv: Beam path selector for filter creation or probing, v: intralipid phantom with US, vi: waveguide and slow light filter (SLF), vii: photomultiplier tube detector with shutter.

Download Full Size | PDF

3.1 Ultrasound measurement

The ultrasonic source was an EPIQ-7 ultrasound machine used with the transducers X5-1 and L12-3 (Philips Medical Systems (PMS), Bothell, WA, USA). The pulses were generated using the pulsed wave Doppler setting in the EPIQ-7 graphical user interface. The X5-1 matrix transducer was used to generate 1.6 MHz center frequency pulses and the L12-3 linear array transducer to generate 3.5 MHz center frequency pulses. The five specific pulses used are denominated by their generating transducer and the length of the pulse as given by the ultrasound device software. Each transducer was separately lowered into a water tank with room temperature water. A needle hydrophone (Precision Acoustics, 0.2 mm diameter aperture) was mounted on a triaxial translation stage powered by stepper motors. Using this stage, the hydrophone was moved to point directly at the transducer where it could be translated transversely to the ultrasound propagation. The pulse pressure in the focal plane of both transducers was then measured with a transverse resolution of 0.5 mm and a >50 MHz sampling frequency. The results from these measurements correspond to Eq. (7). The hydrophone was also positioned at the center of the focus (i.e. at max pressure) to measure the centerline wave form for pulses down to -20 dB of maximum transducer output. A fast Fourier transform of these centerline waveforms produced spectral peaks which relative heights were used to estimate the $a_n$ factors. The described setup and measurements can be seen illustrated in Fig. 1(a).

3.2 Tagging efficiency measurements

The tagging efficiencies were experimentally determined using the setup depicted in Fig. 1(b). First a $7.0\times 7.0\times 3.8\textrm {cm}^{3}$ large tissue phantom was created using a recipe based on the work of Ref. [27] with the mixture proportions 96.5% deionized water, 2.5% Intralipid-20% and 1% agar. As the phantom was mostly made of water, we assume that it is acoustically equivalent to water, i.e. Equations (3)–(4) for water are identical to those for the phantom. The shape and amplitude of the pulses in the phantoms are thus assumed to be the same as what is measured in the water tank. Similar to what was done in Ref. [8], the phantom’s optical properties were measured using photon time of flight spectroscopy. The phantom was determined to have the reduced scattering coefficient $\mu _{\textrm {s}}'={6.4}\textrm {cm}^{-1}$ and, although no absorber was added, the absorption coefficient $\mu _a={0.008}\textrm {cm}^{-1}$.

A Coherent CR-699 dye laser operating with rhodamine 6G tuned to 606nm was used as the light source. The laser was stabilized to a sub-kHz linewidth using the Pound-Drever-Hall technique to lock to a temperature stabilized ultra-low-expansion reference cavity. Using two acousto optic-modulators, we shaped the output of the laser into pulses at different frequencies. See Ref. [28] for a more detailed description of the laser system. The light source was used to create and probe SLFs in a 12 mm thick Pr:YSO crystal. An SLF is a sharp and narrow spectral hole in a highly absorbing material. Other than having high absorption for frequencies outside the hole, the SLF also slows the group velocity of a light pulse with the same frequency as the hole center frequency by several orders of magnitude. This combination allows for both absorptive filtering and time gate filtering of the signal. The prepared filters were 1 MHz wide with a transmission of 60%, a 30 dB contrast for diffuse light and a slow light retardation of ${6}\;\mu \textrm {s}$[8]. After the filter creation, the laser beam was diverted using a flip mirror, allowing for illumination of the tissue phantom. Using an acousto-optic modulators in the pulse shaping optics, the frequency of the laser light to the phantom was shifted to probe either $\mathcal {I}_0$ (no frequency shift), $\mathcal {I}_1$ (a frequency shift of $f_{\textrm {US}}$) or $\mathcal {I}_2$ (a frequency shift of $2f_{\textrm {US}}$). $\mathcal {I}_0$, $\mathcal {I}_1$ and $\mathcal {I}_2$ were measured for different peak US pulse pressures and averaged from 1000 shots at each such peak pressure setting. On the opposite side of the phantom, the light was collected using a liquid core light guide (1cm diameter, numerical $\textrm{aperture} = 0.59$). The collected light was then guided to the SLF for filtration. The filtered light was measured using a photomultiplier tube (Hamamatsu, R943-02).

Over the several hours the measurements were performed, the transmission through the measurement arm of Fig. 1(b) without any ultrasound present (i.e. $\mathcal {I}_{\textrm {ref}}$) slowly deteriorated and was compensated for in the data analysis. This loss in $\mathcal {I}_{\textrm {ref}}$ did not correlate to any loss in laser output power and $\mathcal {I}_{\textrm {ref}}$ was reset to previous values when readjusting the light guide’s contact with the phantom. We therefore believe this transmission loss to be caused by the tissue phantom slowly deforming over time, causing the optical coupling between the light guide and the phantom to deteriorate.

4. Ultrasound fit and simulations

The tagged intensities were simulated using a MATLAB (The MathWorks, Inc., Natick, MA, USA) implemented algorithm, which has been made freely available under the MIT license on our Bitbucket repository (Code 1, [29]). The simulation rests on two cornerstones: the calculation of $M$ optical paths and the US field description. The optical paths were generated using Monte-Carlo (MC) photon transport [30] modified to store scattering points. To match the experimental setup, the MC simulation domain was set to a $7.0\times 7.0\times 3.8\textrm {cm}^{3}$ large cuboid. As absorption was added later when calculating the spectral intensity using Eq. (11), $\mu _a$ was set to zero in the MC simulation. The scattering anisotropy factor $g=0.7$ for intralipid [31] was used to estimate the scattering coefficient $\mu _s=\mu _s'/(1-g) = {21.3}\textrm {cm}^{-1}$ from the measured $\mu _s'$ of the phantom. All packets were initialized in the point $x,y,z=0$ with direction $\hat {\mathbf {e}}_z$, see Fig. 1(b) for axes depiction. Detected packets are defined as all packets leaving the medium through the circle in the $z={3.8}$cm plane with center point $x,y=0$ and 1cm diameter.

In the tagging simulation, the ultrasound pressure is described as a fit of Eq. (8) to the measured ultrasound pulse, disregarding ultrasound harmonic components $n>5$. The following envelope template was deemed to adequately describe all investigated US pulses:

The measurements used to estimate the $a_n$ factors also showed that for the X5-1 4 mm pulse, $\sigma _y$ decreased with the peak pressure from 1.7 mm to 1.4 mm. The same was observed to a smaller extent in the other four pulses which could be neglected. The $\sigma _y$ used in the simulations for the X5-1 4 mm pulse was thus changed according to these measurements using linear interpolation (details in Appendix B). Other than the $a_n$ coefficients and this $\sigma _y$, all other fitting parameters were assumed to be constant in respect to the peak pulse pressure. The measured ultrasound for the X5-1 4mm pulse can be seen in Fig. 2 together with a fitted envelope template from Eq. (18) and $a'_n$ coefficients, where $a_n = a'_n/\sum a'_n$. The full corresponding fitting parameters for the X5-1 4 mm pulse and other pulses can be found in Appendix B.

 

Fig. 2. Characterisation of the X5-1 4mm ultrasound pulse, propagating toward positive $y$. (a) depicts the pulse in the $x = 0$ plane where the dashed lines correspond to the fitted envelope (Eq. (18)). (b) depicts the frontal profile of the pulse in the $xz$-plane where the dashed line corresponds to where the pulse is evaluated in (a), (c) and (d). (c) is (a) seen in the $P_1-z$ plane from $y = {-7}$mm, showing the ultrasound wings along the $z$ axis together with a graphical description of $R_1$ and $R_3$ from Eq. (18). (d) is a top-down view of the pulse in (a) where the dashed line corresponds to the pulse centerline. (e) shows the normalized Fourier transform of the pulse centerline with ${f}_1 = {1.6}\textrm {MHz}$. The peaks of this spectrum correspond to the values $a'_n$, from which $a_n = a'_n/\sum a'_n$ is calculated. In (f), this evaluation is shown for different pulse peak pressures together with the polynomial fit used to estimate the $a_n$ coefficients in the simulations.

Download Full Size | PDF

Once the fitting parameters $R_i$, $\sigma _i$ and $N_i$ ($i =$ placeholder index) are determined, the found envelope is tested so that it follows the requirements discussed in Appendix A after which the displacement field is estimated with Eq. (A.7). With both $P_1$ and $\mathbf {A}$ adequately described, the spectrum of the field at the detector is now simulated along the process depicted in Fig. 3, which details a numerical adaptation of Eqs. (10)–(17). As part of the simulation, time is discretized with a virtual sampling frequency $F_s$ such that $t\rightarrow t_s = (s-1-S/2)/F_s,\ s = 1,2,\ldots ,S$. $F_s= {40}$MHz was used to avoid aliasing for frequency shifts below 20 MHz, i.e. the five first sidebands generated by both US pulses. With $S=1024$, this $F_s$ results in a time span of ${26}\;\mu \textrm {s}$ which gives a frequency resolution of 40kHz and is more than sufficient to cover the ${1}\;\mu \textrm {s}$ optical pulse used in the experiment. A second discretization is performed for each free path which is split into ${Q}$ segments and ${Q}+1$ points such that $l\rightarrow l_q = qL_{m,j}(t_{s})/Q,\ q = 0,1,2,\ldots ,Q$. The drawback of discretizing each free path length in this manner is that since each length is randomly distributed, the spatial step length is not uniformly defined but instead has the median value $\ln (2)/(\mu _s Q)$. The advantages are that no unnecessary estimations of the pressure field is performed and that the memory structures containing the $l_q$ parameter are of fixed size. In the simulations $Q=20$ was used, giving a median resolution of ${16}\;\mu \textrm {m}$. In comparison, the highest frequency component of the 3.5 MHz US pulses has a wavelength of ${84}\mu \textrm {m}$.

 

Fig. 3. Simulation steps following the outline in Sec. 2. (i) Generate paths using Monte Carlo method. (ii) Extract the $M$ paths that arrive to detector and ($iii$) save their scattering positions. (iv) With time discretized into $S$ points, (v) estimate $\mathbf {A}(\mathbf {r},t)$ for all scattering and time points and calculate Eqs. (13), (14) and (17). ($vi$) Discretize the free path length between every scattering point into $Q + 1$ points and ($vii$) estimate Eq. (15). ($viii$) Calculate Eq. (12) with the integral numerically approximated using the trapezoidal method. ($ix$) Find the frequency shifted expected spectrum for each path using fast Fourier transform (FFT). ($x$) Calculate the expected spectrum at the detector and subsequently the tagged fractions.

Download Full Size | PDF

With $E_m(t)$ described numerically, the expectation value of Eq. (10) is given as [16]

In addition to the simulations with $P_1$ estimated from the fitting of Eqs. (8) and (18) to the measured US pulses, the tagged fractions were simulated with the US pulse wings omitted, i.e. $R_1=1,\ R_2 = R_3=0$. This was done in order to understand the effect of the wings on the tagging efficiency. Furthermore, tagged fraction simulations were performed for the measured pulses with either $\mathbf {A}\equiv 0$ or $\frac {\partial \mathfrak {n}}{\partial P}=0$ to find the relative contribution of these two effects. In addition to simulations with the fitted $P_1$, the tagging was, for comparison reasons, simulated with $P_1$ instead estimated directly from the measured US pulses using linear interpolation. Further calculating $\mathbf {A}$ from this interpolated $P_1$ led to the scattering points $\mathbf {r}_{m,j}$ not returning to their equilibrium points after the passage of the US pulse, which was deemed unrealistic. As such, only the case $\mathbf {A}\equiv 0$ is analyzed with the interpolated $P_1$. Lastly, the tagging from the simpler US pulse

5. Results

Figure 4 shows the simulated optical spectrum generated by the max pressure amplitude X5-1 4mm pulse using either the analytically fitted $P_1$ ($\mathbf {A}\not \equiv 0$) or the interpolated $P_1$ ($\mathbf {A}\equiv 0$). The source code for the simulation together with examples is available free

 

Fig. 4. Simulated optical spectrum for the X5-1 4mm pulse in Fig. 2 at the phantom output showing the first four US induced sidebands. $\Delta {f} = f-f_{\textrm {o}}$ is the shift in frequency from the initial optical center frequency. The blue full line is from the fitted $P_1$ expression with $\mathbf {A}\not \equiv 0$ and the red line is the spectrum from an interpolated $P_1$ with $\mathbf {A}\equiv 0$. The simulated sideband strengths $\mathcal {I}_n$ used in Eqs. (1)–(2) are estimated by integrating the simulated spectra from $f_n-{0.5}$MHz to $f_n+{0.5}$MHz.

Download Full Size | PDF

The simulated tagged fractions are compared to the corresponding experimental results in Fig. 5. Other than the predicted tagging from the fitted pulses, Fig. 5 also includes the predicted tagged fractions with the omission of the ultrasound pulse wings ($R_1=1,\ R_2 = R_3=0$).

 

Fig. 5. Tagged fractions as a function of peak pressure for the 5 investigated pulses. Full blue line is the simulated tagged fraction from the fitted US pulse. The dashed green line is the same but with US pulse wings omitted ($R_1=1$, $R_2 = R_3 = 0$). The shaded areas to each line depict the SEM equivalent of the simulations. The red dots depict the mean value and the bars the standard deviations of the experimental measurements.

Download Full Size | PDF

Figure 6 shows the tagged fractions when either the movement of the particles or the pressure induced refractive index change is omitted while Fig. 7 shows the tagged fractions from the simple pulse in Eq. (22) for different ultrasound frequencies.

 

Fig. 6. Tagged fractions as a function of peak pressure for the investigated 5 pulses. Full cyan line is the simulated tagged fraction from the fitted US pulse with $\mathbf {A}\equiv 0$. The dashed purple line are the same except with $\frac {\partial \mathfrak {n}}{\partial P}=0$ instead. The yellow capped marker is the simulated tagged fraction with $P_1$ estimated from interpolation and $\mathbf {A}\equiv 0$. The shaded areas to each line and the capped error bar depict the SEM equivalent of the simulations. The red dots depict the mean value and the bars the standard deviations of the experimental measurements.

Download Full Size | PDF

 

Fig. 7. Tagged fractions of the simple pulse from Eq. (22) as a function of carrier frequency of the ultrasound. The shaded area depict the SEM of the simulation.

Download Full Size | PDF

6. Discussion and outlook

In this section we discuss the results’ possible impact on the UOT technique and also what open questions still exist concerning the topics treated in this paper. Starting with the main result in Fig. 5, the simulated first and second order sideband strengths match well with the measured values while the difference between the simulated and measured $\eta _{\textrm {all}}$ gives rise to some questions. Similar to the results in Ref. [16], $\eta _{\textrm {all}}$ is consistent with the simulation for pressures <0.5 MPa - it is not until we reach higher pressures than those investigated in Ref. [16] that the discrepancy starts to arise. This difference could in part be explained by optical tagging by stray acoustic fields unintentionally output by the transducer or frequency noise in the ultrasound pulse which is not taken into account in the analytical pulse descriptions. From the acoustic measurements, these two factors are not expected to be so large as to account for the entire difference however. As such, we have no clear answer to why this discrepancy exists, though, if the model in this paper is to be used, it is a question that needs to be answered before $\eta _{\textrm {all}}$ can be used reliably in UOT.

In further agreement with Ref. [16] (which demonstrates a higher total tagged fraction due to the acoustic field encompassing a larger fraction of the scattering volume in their experiments), a nonzero amount of the light is shifted into higher order sidebands as $\eta _{\textrm {all}}$ does not equal the sum $2\eta _1 + 2\eta _2$ for neither the simulation nor the measurement results. A further discrepancy along this line is that only the 1st and 2nd order sidebands were strong enough to be directly measured with our setup. In hindsight it is possible that a 3rd order sideband signal could have been measured by averaging from many shots as Fig. 4 shows that the peak of the 3rd order sideband is ’only’ a factor 3 weaker than the 2nd. This was not done since no 3rd order sideband signal was visible at all for single shots at the maximum transducer output pressure. This discrepancy is therefore an avenue for future investigation, but the apparent strength of the higher order sidebands indicate that detection techniques in UOT where multiple sidebands are measured simultaneously may be superior due to the increased signal strength.

The US pulse irregularities (i.e. the wings) bring some insight to the sensitivity of the ultrasound tagging process. As seen in Fig. 5, disregarding the wings lead to a substantial reduction in predicted sideband strengths for the X5-1 transducer pulses while staying within the error margin for the L12-3 transducer pulses. In comparison, the wings are roughly 3 times as wide as the central lobe for both transducers. Comparing the pulse sizes when $R_1 = 1$, $R_2 = R_3 = 0$ with the fitted values (see Table 1), we find that for the pulse with the smallest wing volume (L12-3 2 mm), the wings still comprise 24% of the normalized pressure integral $\int \Gamma dV / \Gamma _0$. For the other pulses, the same ratios correspond to over 42% 49% and 79% for the X5-1 4, 2 and 0.5 mm pulses and 47% for the L12-3 4 mm pulse. This seems to roughly correspond to the signal drop in Fig. 5 when disregarding the wings for the X5-1 transducer pulses, though not for the L12-3 4 mm pulse. The UOT resolution is therefore expected to be worse due to the existence of these wings, though this needs to be validated with, for example, simulations of a medium with inhomogeneous absorption. Through a deeper study of the tagged fraction with respect to the spatial overlap between the US and each optical path, it may also be possible to deconvolute a local tagging rate $\eta (\mathbf {r},t)$ which could be used in the (compared to Monte Carlo) computationally fast light diffusion model. These investigations lie outside the scope of this paper however, even though the current treatment of the acousto-optic effect allows for it. Either way, the observed signal drop when disregarding the wings indicates that proper processing of an UOT image requires that the used US pulse shape is known at all points in the imaging sequence.

Similar to the results reached by Ref. [16], the two tagging effects do not add linearly and while the movement of scattering particles does have an independent effect, Fig. 6 shows that the tagging process can be well described with the change in refractive index alone. In fact, comparing the $\mathbf {A}\equiv 0$ results with the full simulation results in Fig. (5) show that they differ by less than 4%, which is well within the error bounds. This is in contrast to the findings presented in Ref. [16], which found that the tagging contribution by the movement of scatterers is stronger and not negligible. We attribute this difference to a possible error in their calculations as the transverse envelope is not employed for the acoustic amplitude $\textbf {A}$ in the example code made available by the authors of Ref. [16] on GitHub (Ref. [32]). The differences between our results disappear when we emulate the simulations performed in Ref. [16] with a transversely uniform acoustic amplitude. We therefore conclude that the movement of scatterers may be wholly ignored while still achieving accurate results which, if universally true, allows for speed up of the simulations. Again disregarding the particle displacement, the tagged fractions may also be simulated to a somewhat accurate level using $P_1$ interpolated from the measured pulses as seen in Figs. 4 and 6. Accepting this, the time consuming analytical fitting of the US pulses can be bypassed entirely, a task mainly performed in order to ensure that the scattering points return to their equilibrium points after passage of the US.

A final topic for discussion is how the choice of ultrasound frequency could impact an UOT measurement. A higher frequency allows for a higher pressure without exceeding medical safety limits [24], smaller minimum pulse lengths and a larger frequency separation of untagged and tagged light. This in turn would mean higher tagged fraction of light (if tagging efficiency is frequency independent), better spatial resolution and less demands on the SLF materials. The higher ultrasound frequency does have the drawback that ultrasonic absorption in tissue is dependent on frequency [23]. Raman-Nath theory predicts that the strength of the first order acousto-optic effect is unaffected by the ultrasound frequency [21], something that is not apparent from Fig. 5. The lower tagging performance of the higher frequency pulses can instead be explained by them occupying a smaller volume than their 1.6 MHz counterparts. Comparing the normalized pressure integrals show that the X5-1 2 mm and the L12-3 4 mm pulses are the two closest in size, where the latter is roughly 25 % larger than the former. Even so, the 1.6 MHz pulse exhibit a slightly higher tagging efficiency at the same peak pressure. Other than frequency, the major difference between the two pulses are their shapes where L12-3 4mm is longer and narrower while X5-1 2 mm is wider and shorter, something that could be equally impactful. We can therefore not draw any hard conclusions regarding the existence of a fundamental difference in tagging efficiency between the two frequencies from the presented experimental data alone. The results in Fig. 7 shed some light to this, with the caveat that the pulse used is unrealistic in the sense that it only has a single harmonic at those pressures. With that in mind, the results in Fig. 7 agree to an extent with Raman-Nath theory as $\eta _1$ is weakly dependent on the US frequency while simultaneously agreeing with Ref. [16], which found that higher frequencies show a lower $\eta _{\textrm {all}}$. The drop of $\eta _{\textrm {all}}$ for high frequencies is instead explained by more sidebands being nonzero for lower frequencies. As $\eta _{\textrm {all}}$ is the sum of all sidebands then if they decrease linearly with frequency, as Fig. 7 suggests, $\eta _{\textrm {all}}$ as a function of frequency will be curved as sideband intensity cannot be less than zero. That lower frequency pulses generate more sidebands is also consistent with the comparison between X5-1 2 mm and L12-3 4 mm, as the difference in $\eta _{\textrm {all}}$ between the two is larger than their difference in $\eta _1$ or $\eta _2$. As a final point, the increased tagging into the higher order sidebands for lower frequencies suggest that techniques where the UOT signal is measured from multiple sidebands may benefit from the usage of a lower US frequency.

7. Conclusions

In this paper, first a framework for how to locally characterize and describe high pressure pulses from measurements have been explored. Then using this framework, measured ultrasonic pulses have been analytically described in order to accurately predict detected individual optical sidebands generated by acousto-optic interaction in optically scattering media. Albeit that the individual sidebands were accurately predicted, a discrepancy was found between simulation and experiments for the total amount of light removed from the initial optical frequency for peak pressures above 0.5 MPa. We further found that the UOT signal strength is somewhat sensitive to aberrations of the US pulse and as such, full knowledge of the pulse shape is necessary to accurately predict the signal. Furthermore, the ultrasonically induced particle displacement’s effect on the tagging process in the investigated cases is negligible, thus allowing it to be ignored in these contexts. Lastly, we found that the frequency dependence of the tagging process of individual sidebands is weakly dependent on frequency with a stronger dependency for the total tagging efficiency. This is a result that is not fully consistent with Raman-Nath theory but aligns with more recent publications. The results in this paper are important for UOT design aspects, where a validated model describing pulsed acousto-optic interaction can be used as a virtual lab to investigate UOT imaging using arbitrary US pulses in arbitrary tissue types, enabling further development of the imaging technique outside an actual lab environment.

Appendix A: Approximate solutions to Euler’s equations for pulsed plane wave

The following appendix sets out to investigate and discuss an approximate solution to Eq. (4) for a known $P$. We start by rewriting Eq. (4) to

(A.1)$$\rho \left[\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v}\cdot \nabla) \mathbf{v}\right] ={-} \nabla P$$

where the left hand side describes the movement of the medium and the right hand side the force acting on the medium. As we presume we can describe $P= P_0 + P_1(\mathbf {r},t)$ with $P_1(\mathbf {r},t)$ from Eq. (8), we find ourselves with an inhomogeneous partial differential equation. We can from a physics standpoint disregard any general solutions to Eq. (A.1) as we do not expect the medium to move unprovoked. By applying the slow amplitude approximation to the ultrasound pulse, here quantitatively expressed as

(A.2)$$K\left|\frac{\partial^{i} \Gamma}{\partial r^{i}}\right| \gg \left|\frac{\partial^{(i+1)} \Gamma}{\partial r^{(i+1)}}\right|\ ,\ r = x,\ y\ \textrm{or}\ z$$

for any integer $i=0,1,2,\ldots$, we can simplify the right hand side of Eq. (A.1) to

(A.3)$$\begin{aligned}-\nabla P ={-}\Bigg[ \hat{\mathbf{e}}_x\frac{\partial \Gamma}{\partial x}\sum_n a_n \sin \zeta_n &+ \hat{\mathbf{e}}_y \left(\frac{\partial \Gamma}{\partial y}\sum_n a_n \sin \zeta_n - K\Gamma\sum na_n \cos \zeta_n \right) \\ &+ \hat{\mathbf{e}}_z\frac{\partial \Gamma}{\partial x}\sum_n a_n \sin \zeta_n\Bigg] \approx \hat{\mathbf{e}}_y K\Gamma\sum na_n \cos \zeta_n \end{aligned}$$

where we have introduced $\Gamma = \Gamma (x-x_0,Ct-y+y_0,z-z_0)$ and $\zeta _n = nK(Ct-y+y_0) + \varphi _n$ to shorten notation. In Eq. (A.3), in addition to the simplification of the $\hat {\mathbf {e}}_y$ term, we also assume that the force, and therefore also the movement, in $x$ and $z$ is negligible in comparison to that in $y$ by the same criteria (i.e. Eq. (A.2)), even though the movement in the prior two directions is not necessarily zero. If we after this simplification attempt to insert the low amplitude solution $\mathbf {v} = \hat {\mathbf {e}}_y(\rho C)^{-1}P_1$ as an ansatz into Eq. (A.1) and again ignore terms according to Eq. (A.2) (with harmonic sum factors $\approx 1$) we are left with the left hand side

(A.4)$$K\Gamma \sum n a_n\cos\zeta_n\left(1 - \frac{\Gamma}{\rho C^{2}}\sum_n a_n \sin \zeta_n \right)\ .$$

There is approximate equality between the left and right hand sides of Eq. (A.1) for this ansatz if we introduce the second criteria

(A.5)$$\frac{\Gamma}{\rho C^{2}} \ll 1\ .$$

The inclusion of higher order derivatives in Eq. (A.2) becomes apparent when we attempt to calculate the displacement field using Eq. (9). After applying integration by parts twice to Eq. (9) with $\mathbf {v} = \hat {\mathbf {e}}_y(\rho C)^{-1}P_1$, we have

(A.6)$$\begin{aligned}\mathbf{A}(\mathbf{r},t) = \hat{\mathbf{e}}_y &\Bigg(\Bigg[- \frac{1}{\rho KC^{2}}\Gamma^{*} \sum_n \frac{a_n}{n}\cos\zeta_n^{*} \Bigg]_{\tau ={-}\infty}^{\tau=t} + \Bigg[- \frac{1}{\rho K^{2}C^{2}}\frac{\partial \Gamma^{*}}{\partial y} \sum_n \frac{a_n}{n^{2}}\cos\zeta_n^{*} \Bigg]_{\tau ={-}\infty}^{\tau=t} \\ &+ \int_{-\infty}^{t}\frac{1}{\rho K^{2}C}\frac{\partial^{2} \Gamma^{*}}{\partial y^{2}} \sum_n \frac{a_n}{n^{2}}\cos\zeta_n^{*}d\tau\Bigg) \end{aligned}$$

where $\Gamma ^{*}$ and $\zeta _n^{*}$ are the same shorthands as before but with $t$ replaced with $\tau$. From Eq. (A.6) we can see that if the integral can be repeatedly expanded by integration by parts, each expansion $i$ yields a term $K^{-i}(\partial ^{(i+1)} \Gamma ^{*})/(\partial y^{(i+1)})$ times some decreasingly smaller sum of harmonics. Adding that $\Gamma (|\mathbf {r}|=\infty ) = 0$ (i.e. that the pulse is localized) we find from the criteria stated in Eq. (A.2) that

(A.7)$$\mathbf{A}(\mathbf{r},t)\approx \frac{\hat{\mathbf{e}}_y}{\rho KC^{2}}\Gamma(x-x_0,Ct-y+y_0,z-z_0) \sum_n \frac{a_n}{n}\sin\left[nK(Ct-y+y_0)+\varphi_n + \frac{3\pi}{2}\right]$$

if the pressure wave is described by Eq. (8). An important reminder is that this analysis rests on that $a_n$ can be seen as constant over the duration of interest as discussed in Section 2., i.e. the acousto-optic interaction time.

Appendix B: Additional information concerning fitted US pulses

In this appendix we detail the fitting parameters used in the analytical fitting of the US pulses.

In Table 1, the values of the fitting parameters for Eq. (18) are given for pulses at peak pressure. The 1.6 MHz pulses are generally larger in both the central lobe and the wings in comparison to the 3.5 MHz pulses. This is probably due to the two transducers used to generate the pulses are designed for different purposes. The X5-1 (1.6 MHz) transducer’s main application being heart imaging means that a trade off in regard to ultrasound image resolution (i.e. pulse size) and depth has been made in order to ensure that a sufficiently large echo returns from heart depths. Conversely, the L12-3 (3.5 MHz) is designed for imaging of the abdomen and the vascular system at shallower depths, which allows for smaller pulse sizes. The higher frequency of the L12-3 also allows for a higher peak pressure as the medical pressure limit is scaled with $f_{\textrm {US}}^{-1/2}$[24].

Table 1. Fitting parameters for the envelope template in Eq. (18) for the investigated pulses at maximum pressure.

View Table

Figure 8 depicts the change in $\sigma _y$ of the pulse centerline measurements for the X5-1 4 mm pulse. In the simulation of the pulses at lower than max peak pressure, $\sigma _y$ for the X5-1 was estimated from linear interpolation of these values with truncation to the lowest measured $\sigma _y$ for simulations at peak pressures lower than the ones measured in Fig. 8.

 

Fig. 8. Change of pulse length $\sigma _y$ from measurements of the X5-1 4mm pulse centerline at different pressures.

Download Full Size | PDF

Funding

Knut and Alice Wallenberg Foundation (KAW) (2016.0081); Lund University Faculty of Engineering; Wallenberg Centre for Quantum Technology (WACQT) funded by the Knut and Alice Wallenberg Foundation (2017.0449); Smarter Electronic Systems, a collaborative venture by Vinnova, Formas, and Swedish Energy Agency (2019-02110).

Acknowledgments

The authors would like to thank Professor Philip R. Hemmer for the loan of the slow light filter crystal. We would also like extend our gratitude to Yujia Huang for publishing her simulation source code on GitHub.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References

1. L. V. Wang, “Ultrasonic modulation of scattered light in turbid media and a potential novel tomography in biomedicine,” Photochem. Photobiol. 67(1), 41–49 (1998). [CrossRef]  

2. J. Gunther and S. Andersson-Engels, “Review of current methods of acousto-optical tomography for biomedical applications,” Front. Optoelectron. 10(3), 211–238 (2017). [CrossRef]  

3. J. Li and L. V. Wang, “Methods for parallel-detection-based ultrasound-modulated optical tomography,” Appl. Opt. 41(10), 2079–2084 (2002). [CrossRef]  

4. F. Ramaz, B. C. Forget, M. Atlan, A. C. Boccara, M. Gross, P. Delaye, and G. Roosen, “Photorefractive detection of tagged photons in ultrasound modulated optical tomography of thick biological tissues,” Opt. Express 12(22), 5469–5474 (2004). [CrossRef]  

5. K. Zhu, B. Zhou, Y. Lu, P. Lai, S. Zhang, and Y. Tan, “Ultrasound-modulated laser feedback tomography in the reflective mode,” Opt. Lett. 44(22), 5414–5417 (2019). [CrossRef]  

6. Y. Li, P. Hemmer, C. Kim, H. Zhang, and L. V. Wang, “Detection of ultrasound-modulated diffuse photons using spectral-hole burning,” Opt. Express 16(19), 14862–14874 (2008). [CrossRef]  

7. A. Kinos, Q. Li, L. Rippe, and S. Kröll, “Development and characterization of high suppression and high etendue narrowband spectral filters,” Appl. Opt. 55(36), 10442–10448 (2016). [CrossRef]  

8. A. Bengtsson, D. Hill, M. Li, M. Di, M. Cinthio, T. Erlöv, S. Andersson-Engels, N. Reistad, A. Walther, L. Rippe, and S. Kröll, “Characterization and modeling of acousto-optic signal strengths in highly scattering media,” Biomed. Opt. Express 10(11), 5565–5584 (2019). [CrossRef]  

9. H. Zhang, M. Sabooni, L. Rippe, C. Kim, S. Kröll, L. V. Wang, and P. R. Hemmer, “Slow light for deep tissue imaging with ultrasound modulation,” Appl. Phys. Lett. 100(13), 131102 (2012). [CrossRef]  

10. C. Venet, M. Bocoum, J.-B. Laudereau, T. Chaneliere, F. Ramaz, and A. Louchet-Chauvet, “Ultrasound-modulated optical tomography in scattering media: flux filtering based on persistent spectral hole burning in the optical diagnosis window,” Opt. Lett. 43(16), 3993–3996 (2018). [CrossRef]  

11. W. Leutz and G. Maret, “Ultrasonic modulation of multiply scattered light,” Phys. B 204(1-4), 14–19 (1995). [CrossRef]  

12. M. Kempe, M. Larionov, D. Zaslavsky, and A. Z. Genack, “Acousto-optic tomography with multiply scattered light,” J. Opt. Soc. Am. A 14(5), 1151–1158 (1997). [CrossRef]  

13. L. Wang, “Mechanisms of ultrasonic modulation of multiply scattered coherent light: a Monte Carlo model,” Opt. Lett. 26(15), 1191–1193 (2001). [CrossRef]  

14. L. Wang, “Mechanisms of ultrasonic modulation of multiply scattered coherent light: an analytic model,” Phys. Rev. Lett. 87(4), 043903 (2001). [CrossRef]  

15. S. Sakadžić and L. V. Wang, “High-resolution ultrasound-modulated optical tomography in biological tissues,” Opt. Lett. 29(23), 2770–2772 (2004). [CrossRef]  

16. Y. Huang, M. Cua, J. Brake, Y. Liu, and C. Yang, “Investigating ultrasound–light interaction in scattering media,” J. Biomed. Opt. 25(02), 1–12 (2020). [CrossRef]  

17. S. Sakadzic and L. Wang, “Modulation of multiply scattered coherent light by ultrasonic pulses: an analytical model,” Phys. Rev. E 72(3), 036620 (2005). [CrossRef]  

18. G. D. Mahan, W. E. Engler, J. J. Tiemann, and E. Uzgiris, “Ultrasonic tagging of light: theory,” Proc. Natl. Acad. Sci. 95(24), 14015–14019 (1998). [CrossRef]  

19. J.-B. Laudereau, A. A. Grabar, M. Tanter, J.-L. Gennisson, and F. Ramaz, “Ultrafast acousto-optic imaging with ultrasonic plane waves,” Opt. Express 24(4), 3774–3789 (2016). [CrossRef]  

20. P. Lai, R. A. Roy, and T. W. Murray, “Quantitative characterization of turbid media using pressure contrast acousto-optic imaging,” Opt. Lett. 34(18), 2850–2852 (2009). [CrossRef]  

21. M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (7th Edition) (Cambridge University Press, 1999), 7th ed.

22. B. O. Enflo and C. M. Hedberg, Theory of nonlinear acoustics in fluids, vol. 67 of Fluid mechanics and its applications (Kluwer Academic Publishers, 2002).

23. A. Pierce, Acoustics: An Introduction to Its Physical Principles and Applications (Springer International Publishing, 2019).

24. G. ter Haar, “Ultrasonic imaging: safety considerations,” Interface Focus. 1(4), 686–697 (2011). [CrossRef]  

25. E. Fubini-Ghiron, “Anomalies in acoustic wave propagation of large amplitude,” Alta Frequenza 4, 530–581 (1935).

26. W. A. Riley and W. R. Klein, “Piezo-optic coefficients of liquids,” J. Acoust. Soc. Am. 42(6), 1258–1261 (1967). [CrossRef]  

27. R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli, and G. Valentini, “A solid tissue phantom for photon migration studies,” Phys. Med. Biol. 42(10), 1971–1979 (1997). [CrossRef]  

28. Q. Li, “Quantum memory development and new slow light applications in rare-earth-ion-doped crystals,” Ph.D. thesis, Lund University (2018). https://lup.lub.lu.se/search/publication/a16f47ae-8cef-45d4-95d8-db4e7c72dcd3.

29. D. Hill, “Tagged light simulator,” Bitbucket2021, https://bitbucket.org/qilund/tagged-light-simulator.

30. S. Prahl, M. Keijzer, S. Jacques, and A. Welch, “A Monte Carlo model of light propagation in tissue,” Proc. SPIE 10305, 1030509 (1989). [CrossRef]  

31. R. Michels, F. Foschum, and A. Kienle, “Optical properties of fat emulsions,” Opt. Express 16(8), 5907–5925 (2008). [CrossRef]  

32. Y. Huang, “Ultrasound tagging,” https://github.com/yjhuangcd/ultrasound-tagging.