## 1.

## Introduction

Because many illnesses are accompanied by a change in tissue elastic properties, medical practitioners often use palpation to “feel” the presence of disease. This underlies the motivation behind the development of biomedical imaging modalities that contrast soft-tissue elasticity,^{1} which has particular relevance in diseases of the liver,^{2} abdomen,^{3} blood vessels,^{4} and breast.^{5} Current methods for elastography include acoustic radiation force impulse imaging,^{6, 7} sonoelastography,^{8} and magnetic resonance (MR) elastography,^{9} the latter two of which work by monitoring low-frequency (typically, <1 kHz) elastic waves inside the tissue. In comparison to ultrasonic and magnetic resonance techniques, optical detection techniques have the potential to be lower cost and provide real-time imaging over a larger tissue surface area. In particular, coherent optical techniques, such as optical coherence elastography (OCE),^{10, 11, 12} detect nanometer tissue displacements and thus have the potential to provide improved sensitivity and resolution. Recent methods in OCE have employed sinusoidal actuation for higher specificity^{13} and for fitting the tissue response to elastic wave models^{14, 15} in analogy to the ultrasonic and MR methods. However, although OCE maps internal tissue deformations, its depth penetration is limited to only a few millimeters.^{16} In order to circumvent this limitation to direct optical measurement of internal tissue deformation, we propose to probe the depth-dependent tissue stiffness by measuring the elastic waves only on the tissue surface. This method is analogous to nondestructive testing techniques in materials analysis^{17, 18, 19, 20} and seismology,^{21} where internally mapped stiffness properties are inferred from surface acoustic waves (SAWs).

In fact, the use of SAWs for elastography of soft tissues is a growing area of interest,^{22, 23} but many of these experiments have been accomplished using a laser Doppler vibrometer, which cannot rapidly assess the entire tissue surface. Early studies using stroboscopic photography, while lacking nanoscale motion sensitivity clearly reveal the radial propagation of SAWs in human tissue, under point-source excitation.^{24} In comparison, holographic interferometry is capable of imaging nanoscale deformation over a large sample surface and has been widely used for vibrometry in materials analysis.^{25, 26} This technique has also been employed on biological tissues, where it was suggested that the method may be suitable for elasticity measurements.^{27} Related work using speckle pattern imaging demonstrates the detection of skin lesions based on their distinct mechanical properties.^{28} These techniques may also be suitable for breast cancer detection, because most malignant breast phenotypes are associated with an increased elastic modulus.^{29} The potential for holographic detection of breast abnormalities was demonstrated in human patients under quasi-static deformation.^{30} Under sinusoidal excitation, a more recent study of breast tissue phantoms revealed variations in the surface displacement as a function of the position of a tumorlike inhomogeneity.^{31} Although these results are promising, currently only empirical results have been demonstrated. Furthermore, because tumor location and breast size are highly variable, a versatile method with the ability to image at different depths is desirable.

In this paper, we describe a digital holographic interferometry method for elastography based on the quantitative description of SAWs as Rayleigh waves on an elastic half-space. Although this is only a first step toward 3-D reconstruction of elastograms, it provides needed intuition as to the nature of SAWs. In contrast to previous work,^{25, 26, 27, 30, 31} we employ a continuous, rather than a pulsed (i.e., stroboscopic), laser source, which is enabled by the increasing availability of high-speed digital cameras to capture fast tissue motion. We also developed a phase-shifting algorithm in the time domain for use with our on-axis holography setup, instead of the traditional Fourier transform or spatial phase-shifting techniques, which require careful control of the fringe spacing via off-axis holography.^{25, 32} A distinguishing feature of this algorithm is that it directly provides the phase of the mechanical wave (modulo π) using *a priori* knowledge of the SAW frequency, without the need to shift the reference beam phase. However, the key difference between this and previous work is that we measure the frequency-dependent velocity (i.e., the velocity dispersion of SAWs). This is important because the depth-dependent attenuation of the out-of-plane motion of SAWs is proportional to the SAW wavelength.^{33, 34} By sweeping the wavelength of SAWs, we therefore probe the tissue at varying depths, which enables reconstruction of the depth-dependent tissue elasticity. In this paper, we study the feasibility of this system for quantitative elastography using homogeneous and two-layer soft-tissue phantoms. We also demonstrate empirical results in phantoms with a local inhomogeneity to investigate the potential of this holographic interferometry technique for breast cancer detection.

## 2.

## Methods

## 2.1.

### Tissue-Phantom Preparation

Tissue phantoms with optical and mechanical properties comparable to that of soft tissues were used in these experiments, similar to those used previously.^{35} For the homogeneous phantoms, they were fabricated from a mixture of pure non–cross-linking polydimethylsiloxane (PDMS) fluid (50-cSt viscosity, ClearCo, Inc., Bensalem, Pennsylvania), and a cross-linking room-temperature PDMS base and its associated curing agent (Sylgard 184 A and B, respectively, Dow Corning, Inc. Midland, Wisconsin). The Young's moduli of the phantoms were adjusted over the range from 2 to 40 kPa by changing the ratio of the pure PDMS and the PDMS base over the range from 16:1 to 5:1, whereas the ratio between the PDMS base and its curing agent was held constant at 10:1. The relation between the mixture ratio and Young's modulus was estimated based on previous work.^{36} Titanium dioxide (TiO_{2}) powder [Sigma-Aldrich Co., LLC (St. Louis, Missouri), No. 224227, mean size 1 μm, <5 μm] was embedded at a concentration of 4 mg/g in the tissue phantoms to function as tissuelike optical scatterers for holographic imaging. The phantom solutions were mixed thoroughly on an ultrasonicator for several tens of minutes at room temperature and then poured into a 126-mm-diam cylindrical container, with a thickness of ∼30 mm. The phantom dimensions were chosen to be large compared to the attenuation length of Rayleigh waves to eliminate significant boundary reflections. To further eliminate the boundary reflections, a few millimeters of the phantom adhering to the sidewall were removed. All phantoms were cured in an oven at 80°C for one day and subsequently at room temperature for another day. To prepare the two-layer phantoms, the same process was repeated for the top layer. The solution mixture was poured onto the cross-linked lower layer, and the whole phantom was cured as for the homogeneous phantoms. To construct a phantom with a local inhomogeneity, we prepared two homogeneous phantoms with the same mixing ratio and stacked them together, between which a thin metal ring with a 6.6 mm o.d., 3 mm i.d., and thickness of 0.5 mm was sandwiched.

## 2.2.

### Digital Holographic Interferometry

Figure 1 shows the in-line, image-plane holographic interferometry system for monitoring SAWs. A 633-nm helium-neon laser of 5 mW power (Thorlabs) was divided into sample and reference arms by a 92:8 beam splitter, where the majority of the optical power was directed to the sample arm. Two lenses in the sample arm acted as a 2.5× beam expander with the lens spacing adjusted so that the output of the last lens was slightly divergent, illuminating the phantom with a beam nominally 30 mm in diameter. A 2–in.-diam mirror (downward reflecting mirror in Fig. 1) directed the beam downward onto the horizontal surface of the phantom and subsequently directed the diffusely reflected light from the phantom through a 200-mm focal length, 2–in.-diam imaging lens, which imaged the light onto the sensor of a fast CMOS camera (Photron SA3). The phantom was mounted on a vertical stage with which the height was adjusted to focus the image. The camera frame rate was adjusted between 1 and 7.5 kHz to sample images at a rate ≥10× the frequency of the SAW. At these frame rates, the corresponding pixel array size varied between 1024 × 1024 and 256 × 256, respectively. The pixel size is 17 μm, and the imaging area for the full 1024 × 1024 pixel array was 16 × 16 mm^{2}. Because there are highly spatially modulated speckles in the image, an adjustable aperture was placed close to the imaging lens. The aperture size was adjusted until the speckles could be distinguished by the camera sensor. The reference arm beam was passed through a spatial filter (i.e., the combination of a short focal length lens and a pinhole) and a 150-mm focal-length positive lens (f_{4} in Fig. 1), before being directed into the camera. The distance from the pinhole to the lens f_{4} was adjusted so that interference fringes between the reference and sample could be identified over the whole measurement area. Because the optical path of an edge ray from the sample to the camera is greater than that of a central ray, the output light from the lens f_{4} was chosen to be divergent. Ideally, the optical path difference (OPD) between edge rays and central rays should be balanced for the two arms. A neutral density filter was inserted in the reference arm to optimize the interference signal within the dynamic range of the camera. A retroreflector was used in the reference arm to adjust the OPD between the two arms so that it was within the coherence length of the laser, ∼200 mm.

The image of the phantom interacted with the reference beam, forming a local holographic image. As SAWs propagated across the phantom surface, surface displacement along the axis that bisects the propagation directions of the incident and detected scattered fields (approximately the axis normal to the phantom surface) induced a phase shift in the interferogram due to the changing OPD. As such, this system was predominantly sensitive to out-of-plane sample motions. For generation of SAWs, a piezoelectric transducer (Kinetic Ceramics, Inc., Hayward, California) with a travel range of 20 μm was placed into contact with the phantom immediately to one side of the illumination area. The transducer contact area could be regarded as a point because its diameter is much less than the SAW wavelengths in our experiments. The transducer was driven by a high-voltage power supply controlled by a function generator. Typically, the transducer was initiated at a fixed frequency and the system was allowed to reach a steady state (∼1−2 s), after which the camera recorded the holographic images in real time. On-board image data were transferred to the computer for further processing, as described later.

## 2.3.

### SAW Phase Extraction

The SAWs propagate outward from the point of actuation, which cause temporally modulated phase shifts in the holographic image. When the phantom's thickness is much greater than the SAW wavelength, we expect the SAW to be characterized as a Rayleigh wave.^{33} In that case, the displacement in the *z* direction (out of plane) can be expressed as
[TeX:]
$\Delta z\left({\vec r,t} \right) = A\cos ({\omega t + \vec k \cdot \vec r})$
$\Delta z\left(\stackrel{\u20d7}{r},t\right)=A\mathrm{cos}\left(\omega t+\stackrel{\u20d7}{k}\xb7\stackrel{\u20d7}{r}\right)$
, where *A* is the SAW amplitude (assumed to be slowly varying in
[TeX:]
$\vec r$
$\stackrel{\u20d7}{r}$
, roughly as the inverse square root of
[TeX:]
$\left| {\vec r} \right|$
$\left|\stackrel{\u20d7}{r}\right|$
),
[TeX:]
$\vec r$
$\stackrel{\u20d7}{r}$
is the position in the (*x*,*y*) surface plane, ω is the angular frequency of the SAW, and
[TeX:]
$\vec k$
$\stackrel{\u20d7}{k}$
is the spatial frequency of the SAW. The intensity recorded on the camera can then be written as follows:

## Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} I\left({\vec r,t} \right) = I_{{\rm DC}} + I_{{\rm AC}} \cos [ {\varphi _{{\rm SAW}} \cos ({\omega t + \vec k \cdot \vec r}) + \phi ({\vec r})}], \end{equation}\end{document} $$I\left(\stackrel{\u20d7}{r},t\right)={I}_{\mathrm{DC}}+{I}_{\mathrm{AC}}\mathrm{cos}\left[{\varphi}_{\mathrm{SAW}}\mathrm{cos}\left(\omega t+\stackrel{\u20d7}{k}\xb7\stackrel{\u20d7}{r}\right)+\phi \left(\stackrel{\u20d7}{r}\right)\right],$$_{SAW}= 4π

*A*/λ is the amplitude of the optical phase shift induced by the SAW, [TeX:] $\phi \left({\vec r} \right)$ $\phi \left(\stackrel{\u20d7}{r}\right)$ is the optical phase difference between the reference and sample fields in the absence of the SAW,

*I*

_{DC}is the sum of the reference and sample beam intensities,

*I*

_{AC}is the cross-spectral density of the reference and sample beams, and λ is the laser wavelength.

To extract the phase of the SAW, that is,
[TeX:]
$\beta _{{\rm SAW}} \left({\vec r,t} \right)\break = \omega t + \vec k \cdot \vec r$
${\beta}_{\mathrm{SAW}}\left(\stackrel{\u20d7}{r},t\right)=\omega t+\stackrel{\u20d7}{k}\xb7\stackrel{\u20d7}{r}$
, in the limit of small SAW amplitudes (Δ*z* ≪ λ) and sufficiently fast temporal sampling (ωΔ*t* ≪ 2π, where Δ*t* is the time between successive frames), it is possible to use a temporal phase-shifting algorithm similar to Refs. 37 and 38, such that

## Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \beta _{{\rm SAW}} ({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t}) &=& \tan ^{ - 1} \left[\! {\frac{{I({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t {-} \Delta t}) {-} I({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t {+} \Delta t})}}{{I({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t {-} \Delta t}) {+} I({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t {+} \Delta t}) {-} 2I({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t})}}}\right.\hspace*{-5pt}\nonumber\\ && \times\left. {\tan \left({\frac{{\omega \Delta t}}{2}} \right)} \vphantom{\left[\! {\frac{{I\left({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t {-} \Delta t} \right) {-} I\left({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t {+} \Delta t} \right)}}{{I\left({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t {-} \Delta t} \right) {+} I\left({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t {+} \Delta t} \right) {-} 2I\left({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t} \right)}}}\right.}\right], \end{eqnarray}\end{document} $\begin{array}{ccc}\hfill {\beta}_{\mathrm{SAW}}\left(\stackrel{\rightharpoonup}{r},t\right)& =& {\mathrm{tan}}^{-1}\left[\frac{I\left(\stackrel{\rightharpoonup}{r},t-\Delta t\right)-I\left(\stackrel{\rightharpoonup}{r},t+\Delta t\right)}{I\left(\stackrel{\rightharpoonup}{r},t-\Delta t\right)+I\left(\stackrel{\rightharpoonup}{r},t+\Delta t\right)-2I\left(\stackrel{\rightharpoonup}{r},t\right)}\right.\hfill \\ & & \times \left.\mathrm{tan}\left(\frac{\omega \Delta t}{2}\right)\right],\hfill \end{array}$where a group of three successive frames at times *t* – Δ*t*, *t*, and *t* + Δ*t* are processed to produce the SAW phase image at time *t*, and
[TeX:]
$\mathord{\buildrel{\lower3pt \hbox{\scriptscriptstyle\rightharpoonup}}\over r}$
$\stackrel{\rightharpoonup}{r}$
is the (*x*,*y*) pixel position in the frame. Briefly, we can understand why this works by rewriting the phase-shifting algorithm in terms of time derivatives of *I* and using the small angle approximation,

## Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} \beta _{{\rm SAW}} ({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t}) &=& \tan ^{ - 1} \left[ {\frac{{I'({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t})}}{{I''({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t})}} \cdot \frac{2}{{\Delta t}}\tan \left({\frac{{\omega \Delta t}}{2}} \right)} \right]\nonumber\\ & \approx& \tan ^{ - 1} \left({\frac{{I'({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t})}}{{I''({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t})}}\omega } \right). \end{eqnarray}\end{document} $\begin{array}{ccc}\hfill {\beta}_{\mathrm{SAW}}\left(\stackrel{\rightharpoonup}{r},t\right)& =& {\mathrm{tan}}^{-1}\left[\frac{{I}^{\prime}\left(\stackrel{\rightharpoonup}{r},t\right)}{{I}^{\prime \prime}\left(\stackrel{\rightharpoonup}{r},t\right)}\xb7\frac{2}{\Delta t}\mathrm{tan}\left(\frac{\omega \Delta t}{2}\right)\right]\hfill \\ & \approx & {\mathrm{tan}}^{-1}\left(\frac{{I}^{\prime}\left(\stackrel{\rightharpoonup}{r},t\right)}{{I}^{\prime \prime}\left(\stackrel{\rightharpoonup}{r},t\right)}\omega \right).\hfill \end{array}$*I*in Eq. 1 about φ

_{SAW}(in the limit of small SAW amplitudes), then we obtain

## Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} I\left({\vec r,t} \right) &\approx& I_{{\rm DC}} + I_{{\rm AC}} \cos [ {\phi \left({\vec r} \right)} ]\nonumber\\ && -\, I_{{\rm AC}} \varphi _{{\rm SAW}} \sin [ {\phi \left({\vec r} \right)}]\cos ({\omega t + \vec k \cdot \vec r}), \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill I\left(\stackrel{\u20d7}{r},t\right)& \approx & {I}_{\mathrm{DC}}+{I}_{\mathrm{AC}}\mathrm{cos}\left[\phi \left(\stackrel{\u20d7}{r}\right)\right]\hfill \\ & & -\phantom{\rule{0.16em}{0ex}}{I}_{\mathrm{AC}}{\varphi}_{\mathrm{SAW}}\mathrm{sin}\left[\phi \left(\stackrel{\u20d7}{r}\right)\right]\mathrm{cos}\left(\omega t+\stackrel{\u20d7}{k}\xb7\stackrel{\u20d7}{r}\right),\hfill \end{array}$$## Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \begin{array}{l} I'({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t}) \approx \omega I_{{\rm AC}} \varphi _{{\rm SAW}} \sin [ {\phi ({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r} })}]\sin ({\omega t + \mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over k} \cdot \mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r} }), \\[8pt] I''({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r},t}) \approx \omega ^2 I_{{\rm AC}} \varphi _{{\rm SAW}} \sin [ {\phi ({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r} })} ]\cos ({\omega t + \mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over k} \cdot \mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}} \over r} }). \\ \end{array} \end{equation}\end{document} $$\begin{array}{c}{I}^{\prime}\left(\stackrel{\rightharpoonup}{r},t\right)\approx \omega {I}_{\mathrm{AC}}{\varphi}_{\mathrm{SAW}}\mathrm{sin}\left[\phi \left(\stackrel{\rightharpoonup}{r}\right)\right]\mathrm{sin}\left(\omega t+\stackrel{\rightharpoonup}{k}\xb7\stackrel{\rightharpoonup}{r}\right),\hfill \\ {I}^{\prime \prime}\left(\stackrel{\rightharpoonup}{r},t\right)\approx {\omega}^{2}{I}_{\mathrm{AC}}{\varphi}_{\mathrm{SAW}}\mathrm{sin}\left[\phi \left(\stackrel{\rightharpoonup}{r}\right)\right]\mathrm{cos}\left(\omega t+\stackrel{\rightharpoonup}{k}\xb7\stackrel{\rightharpoonup}{r}\right).\hfill \end{array}$$Plugging these back into Eq. 3 reveals how the extracted phase β_{SAW} using the phase-shifting algorithm in Eq. 2 is equivalent to the phase of the mechanical wave,
[TeX:]
$\omega t + \vec k \cdot \vec r$
$\omega t+\stackrel{\u20d7}{k}\xb7\stackrel{\u20d7}{r}$
. It should be noted that our method is nontraditional in that the reference beam is not phase shifted and the optical phase is not obtained directly; instead, we employ *a priori* knowledge of the mechanical wave frequency to obtain the mechanical phase. This has the advantage of skipping the unnecessary evaluation of the optical phase, in favor of directly assessing the mechanical phase. One consequence is that we do not know at any instant in time whether the SAW is in its upswing or downswing and, therefore, we can only measure β_{SAW} modulo π. In our phase images, therefore, we represent β_{SAW} = –π/2 to π/2 as a gray-scale value from black to white, respectively.

Figure 2a (Video 1) shows an example interference pattern measured by the camera,
[TeX:]
$I({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}}\over r},t})$
$I\left(\stackrel{\rightharpoonup}{r},t\right)$
. In the video, the evenly distributed fringes aligned from lower right to upper left are background fringes caused by the protection glass in front of the camera sensor due to multiple reflections (Fabry–Perot). These are stationary and do not affect the measurement result using the temporal phase-shifting algorithm. The SAW phase image,
[TeX:]
$\beta _{{\rm SAW}} ({\mathord{\buildrel{\lower3pt\hbox{\scriptscriptstyle\rightharpoonup}}\over r},t})$
${\beta}_{\mathrm{SAW}}\left(\stackrel{\rightharpoonup}{r},t\right)$
, which was extracted from Fig. 2a according to Eq. 3, is shown in Fig. 2b (Video 2). It reveals a SAW propagating away from the actuator position located outside of the image to the right. Note that the phase-wrap position from β_{SAW} = π/2 to –π/2 provides a highly contrasted way to visualize the SAW wavefront, which we also verified with simulations. We then obtain the SAW phase velocity by tracking a specific wavefront in the image series. Although we found human inspection to be sufficient for SAW velocity measurement within the scope of this work, in future work we plan to develop an automated method for SAW velocity extraction.

It is important to note that modulation of the speckle is the dominant effect giving rise to the SAW phase signal. This is because the phantoms were designed so that the diffuse scattering from the TiO_{2} micropowder dominates over the specular reflection from the surface in order to more accurately mimic the optical signal from real tissues. This experimental configuration can actually be thought of as a specific case of digital holography known as electronic speckle pattern interferometry.^{37, 39, 40}

## 3.

## Results and Discussion

## 3.1.

### Surface Acoustic Waves in Homogeneous Phantoms with Varying Young's Moduli

Our hypothesis is that SAWs on a soft-tissue–like, homogeneous, elastic half-space will behave as Rayleigh waves.^{33} consistent with other observations.^{23} Here, we study this phenomenon using digital holographic interferometry. We generated tissue phantoms with Young's moduli ranging from 2 to 40 kPa to match the range of elastic moduli encountered in soft tissue such as the human breast.^{29} According to elastic wave theory, the phase velocity of Rayleigh waves in a homogeneous half-space depends on the Young's modulus (*E*) and density of the material (ρ) according to^{23}

## Eq. 6

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} c_{\rm R} \approx \frac{1}{{1.05}}\sqrt {\frac{E}{{3\rho }}}, \end{equation}\end{document} $${c}_{\mathrm{R}}\approx \frac{1}{1.05}\sqrt{\frac{E}{3\rho}},$$To produce SAWs, we excited the phantoms with a 150-Hz sinusoidal signal on the piezoelectric transducer. This choice of frequency was such that the corresponding SAW wavelength (and, thus, penetration depth) was smaller than the phantom depth (30 mm). This produces a surface wave with constant phase velocity due to the elastic homogeneity of the material. Using the SAW phase-extraction method previously described, we quantified the wave speed as a function of Young's modulus and plotted the results in Fig. 3. Each data point was obtained by averaging eight measurements for a single sample. The solid line is the theoretical curve calculated with Eq. 6 using ρ = 960 kg/m^{3}. The dashed line is a least-squares best-fit curve with one fitting parameter, *k*
_{0}, such that *c*
_{R} = *k*
_{0}√E. The fitting value of *k*
_{0} = 0.525 is within 6% of the theoretical value of 0.560. The *R*
^{2} value of the fit was 0.99. For additional validation, we measured the SAW speeds for one homogeneous phantom while sinusoidally driving the transducer at different frequencies (75−300 Hz) and found that the measured velocity remained constant with a standard deviation of ∼7.5% of the mean. This demonstrates that the induced out-of-plane SAW motions measured by holographic interferometry are well described as Rayleigh waves.

## 3.2.

### Surface Acoustic Wave Velocity Dispersion in Two-Layer Phantoms

Because the penetration depth of Rayleigh waves is proportional to the wavelength, it is known that a variation in the material's elasticity versus depth gives rise to dispersion in the SAW velocity.^{17, 18, 19, 20, 21} We hypothesize that SAWs in soft tissues will exhibit a velocity according to Eq. 6 given by the average material elasticity *E* within the surface depth probed by the wave. In other words, we expect the long-wavelength (low-frequency) SAWs to probe the tissue response over a large depth range, while the short-wavelength (high-frequency) SAWs should probe only the upper tissue layers.

In order to test this hypothesis, we generated two two-layer phantoms, one with a soft upper layer and stiff lower layer, and a second with a stiff upper layer and soft lower layer. For the first phantom, the Young's modulus *E* of the upper layer was 7 kPa [corresponding to a theoretical velocity of 1.4 m/s according to Eq. 6], and the lower layer was 22 kPa (corresponding to a velocity of 2.6 m/s). The top layer was 4 mm thick, and the bottom layer was 10 mm thick. The measured SAW speed at different driving frequencies is shown in Fig. 4a, where the dashed lines indicate the theoretical velocities in a homogenous medium composed of the upper or lower layer. For the second phantom, the Young's modulus *E* of the upper layer was 38 kPa (corresponding to the velocity of 3.4 m/s) and the lower layer was 11 kPa (corresponding to the velocity of 1.8 m/s). The top layer was 6 mm thick, and the bottom layer was 30 mm thick. The measured SAW speed at different driving frequencies is shown in Fig. 4b.

From the two plots, it can be seen that the SAW velocity gradually transitions from close to the velocity of the lower layer to that of the upper layer with the increase of the driving frequency. This expected behavior is observed in both cases, with the stiffer layer either being on top or underneath, such that a positive dispersion is indicative of decreasing stiffness versus depth, and a negative dispersion that of increasing stiffness versus depth. In future work, this measured dispersion profile can be used to quantitatively predict the depth-dependent elasticity either in the perturbative regime^{18, 19} or by treating the exact solution as an inverse problem.^{17, 20}

Qualitatively, we can compare the transition wavelength between lower and upper wave velocities to the layer thickness. This is obtained at the point at which the slope d*v*/d*f* is the maximum and is observed to occur at approximately 100 and 450 Hz in each of the phantoms, respectively. The associated SAW wavelengths, using the velocities measured at each of these frequencies, are 18 and 7 mm, respectively. In the first phantom with a 4-mm-thick soft upper layer, the transition wavelength of 18 mm is much larger than the layer thickness, suggesting significant decay of the wave has reduced its effective probing depth. In the second phantom with a 6-mm-thick stiff upper layer, the transition wavelength was ∼7 mm. Depending on Poisson's ratio, out-of-plane Rayleigh wave motion is reduced to ∼20–33% of its surface value at one wavelength in depth.^{34} We might therefore expect the transition to occur when the wavelength is on the order of the upper layer thickness. However, comparing the transition wavelength with upper layer thickness in the two phantoms suggests that the softer phantom material exhibits greater decay than the stiffer material. This suggestion is also supported by the observation that the SAW velocity in the phantom in Fig. 4b does not quite reach that of the soft lower layer at low frequency. If the attenuation in the soft layer is too rapid, then the wave will tend to bias toward a velocity given by the stiffer upper layer. Further investigation is needed to determine whether the softer phantom material exhibits more dissipative loss, which might be expected because the softer phantom contains a lower concentration of cross-linking PDMS (which is dominantly elastic) relative to the non–cross-linking PDMS fluid (which is dominantly viscous).

In experiments, we also note that the observable range of SAW frequencies (frequencies where SAWs had sufficient amplitude to provide a detectable phase variation) were limited in practice. No SAWs were observed at <20 Hz, while at >80 Hz, SAWs were observed in every experiment. We found the high-frequency limit to lie between 300 Hz for phantoms with a soft upper layer and 800 Hz for phantoms with a stiff upper layer. Improvements in phase sensitivity or actuator sweep may increase the detectable frequency range, in order to increase the depth penetration. Thus, we expect the spatial resolution for elastographic reconstruction to be fundamentally limited by the SAW bandwidth and signal-to-noise ratio. To speed the data acquisition, a chirped (frequency-swept) SAW waveform can be applied to rapidly measure the dispersion curve; the rate of chirp can be slowed to provide the desired frequency resolution, depending on the needed spatial resolution.

## 3.3.

### Surface Acoustic Wave Propagation in Phantoms with a Local Inhomogeneity

To investigate the potential of this technique for detecting local inhomogeneities, we prepared an otherwise homogeneous soft phantom (7 kPa) with an embedded thin metal ring as described in Sec.
2.1. The embedding depth was 8 mm below the surface. Figure 5 shows the phase reconstruction of the SAW propagation at frequencies of 65 Hz (Video 3) and 89 Hz (Video 4) using a frame rate of 1 kHz, and 200 Hz (Video 5) using a frame rate of 2 kHz. These frequencies nominally correspond to wavelengths of 22, 16, and 7 mm, respectively. The measurement area was 16 × 16 mm^{2} at 1 kHz and 10 × 10 mm^{2} at 2 kHz. The effect of the ring is clearly observed as a deformation of the SAW propagation, and the shape of this deformation appears to change at different frequencies due to the varying penetration depth. In particular, at the frequency of 200 Hz, the impact of the ring is minimized because its wavelength of 7 mm is shorter than the ring depth of 8 mm, while the impact is larger at the shorter wavelengths. The impact appears to be largest at 89 Hz, although significant wavefront deformation is observed in the upper left corner of the images at 65 Hz as well. This demonstrates that this technique is sensitive to local inhomogeneities at least as deep as 8 mm. The capability to identify local inhomogeneities is useful because tumors and other disease states are often characterized by a palpable, local elastic abnormality.

## 4.

## Conclusions

This paper demonstrates the potential for a new elastography technique for soft tissue by monitoring the out-of-plane motion of SAWs with digital holography. In homogeneous and two-layer soft-tissue phantoms, SAWs display the characteristics of Rayleigh waves. SAW propagation is also clearly affected in response to a stiff inclusion placed 8 mm beneath the tissue phantom surface. As such, the frequency-dependent velocity of SAWs contains information about the depth-dependent tissue elasticity (Fig. 4), whereas the lateral propagation direction of SAWs contains information about the lateral profile of the tissue elasticity (Fig. 5). By deconstructing this system into simplified phantoms, we have provided needed intuition as to how SAWs behave in response to elastic inhomogeneities.

Before this system can be used in living subjects, motion artifacts (e.g., from cardiopulmonary function) will need to be addressed. Previous experiments using modulated contrast methods *in vivo* with phase-resolved optical coherence tomography^{42} suggest that bandpass filtering around the excitation frequency may be sufficient to reject these types of motion artifacts.

We note that current depth-profiling methods using SAWs^{17, 18, 19, 20} have been performed in elastic media, and future work is needed to extend these for use in viscoelastic media (where the mechanical properties of tissue may be a function of frequency).^{43} It may be possible to reduce the viscoelastic parameter space by characterizing the scaling of the tissue frequency dependence with appropriate models such as that in Ref. 44 to enable depth profiling in viscoelastic tissues. In future work, we hope that this will lead to reconstruction methods based on elastic scattering theory to convert the observed SAW behavior into 3-D tissue elastograms.^{45} This purely optical technique is noninvasive, takes advantage of the high displacement sensitivity afforded by interferometry, and provides depth-dependent elasticity information potentially several centimeters into soft tissue. In particular, holographic interferometry allows us to assess a large surface area simultaneously. Our hope is that it may eventually lead to a new clinical tool for quantitative palpation.

## Acknowledgments

We thank Vasilica Crecea, at the University of Illinois at Urbana-Champaign, for providing calibration data for the Young's modulus of the silicone tissue phantoms. This work was supported by startup funds from the University of North Carolina at Chapel Hill.

## References

*In vivo*visualization of abdominal malignancies with acoustic radiation force elastography,” Phys. Med. Biol., 53 279 –293 (2008). https://doi.org/10.1088/0031-9155/53/1/020 Google Scholar

*In vivo*three-dimensional optical coherence elastography,” Opt. Express, 19 (7), 6623 –6634 (2011). https://doi.org/10.1364/OE.19.006623 Google Scholar

*in vivo*acousto-optical elastography,” Opt. Express, 14 (21), 9770 –9779 (2006). https://doi.org/10.1364/OE.14.009770 Google Scholar

*in vivo*optical coherence tomography,” Opt. Express, 13 (17), 6597 –6614 (2005). https://doi.org/10.1364/OPEX.13.006597 Google Scholar

*In vivo*magnetomotive optical molecular imaging using targeted magnetic nanoprobes,” Proc. Natl. Acad. Sci., 107 (18), 8085 –8090 (2010). https://doi.org/10.1073/pnas.0913679107 Google Scholar