Spatial resolution in dynamic optical coherence elastography

Abstract. Dynamic optical coherence elastography (OCE) tracks elastic wave propagation speed within tissue, enabling quantitative three-dimensional imaging of the elastic modulus. We show that propagating mechanical waves are mode converted at interfaces, creating a finite region on the order of an acoustic wavelength where there is not a simple one-to-one correspondence between wave speed and elastic modulus. Depending on the details of a boundary’s geometry and elasticity contrast, highly complex propagating fields produced near the boundary can substantially affect both the spatial resolution and contrast of the elasticity image. We demonstrate boundary effects on Rayleigh waves incident on a vertical boundary between media of different shear moduli. Lateral resolution is defined by the width of the transition zone between two media and is the limit at which a physical inclusion can be detected with full contrast. We experimentally demonstrate results using a spectral-domain OCT system on tissue-mimicking phantoms, which are replicated using numerical simulations. It is shown that the spatial resolution in dynamic OCE is determined by the temporal and spatial characteristics (i.e., bandwidth and spatial pulse width) of the propagating mechanical wave. Thus, mechanical resolution in dynamic OCE inherently differs from the optical resolution of the OCT imaging system.


Introduction
Optical coherence elastography (OCE) uses optical coherence tomography (OCT) to detect local tissue displacement following an applied force, enabling the quantification of biomechanical properties within a region of interest. OCE can potentially identify pathologies opaque to classical OCT using local maps of tissue mechanical properties (namely Young's modulus), which can vary greatly between healthy and disease states. Detecting spatial heterogeneities in a mechanical structure can help monitor and treat disease progression in numerous fields, including ophthalmology, 1-5 dermatology, 6 cardiology, 7 gastroenterology, 8 and oncology. [9][10][11] While it is a general belief that dynamic OCE is capable of high-resolution (on the order of the imaging modality, OCT) mapping of tissue mechanical properties, the spatial resolution in wave-based optical elastography has not yet been characterized.
Elastography, like all medical imaging modalities, can be characterized in terms of spatial resolution (how small a feature can be seen), contrast resolution (how small a difference in the imaging parameter, shear modulus for the case of OCE, can be detected), and artifact level (how well are spatially resolved, high contrast features related to true tissue properties). In general, spatial resolution in an image is defined as the minimum distance between two objects at which they can be considered individuals. [12][13][14][15][16] The conventional definition of spatial resolution based on the Rayleigh criterion is often adopted in OCE and characterized in terms of the system point spread function (PSF) (i.e., the apparent size of a point source in an image). This definition was recently applied in the field of compression OCE. 17 Using this approach, the propagating elastic wave used to infer modulus in dynamic OCE is typically assumed to behave as a step function at a material interface perpendicular to the propagation direction. Thus, the spatial resolution in dynamic OCE is usually defined by the lateral (perpendicular to light propagation direction) resolution of the optical measurement system. As dynamic OCE does not image isolated sources (propagating mechanical waves are tracked to determine local speed), the image depends on factors such as the size of the object being imaged and the characteristics of the mechanical load. Consequently, the image cannot be considered simply as the convolution of an input with a PSF, and the general definition of spatial resolution must be modified.
A more appropriate definition for spatial resolution in dynamic OCE is thus the minimum size of an inclusion, in a homogeneous background, for which propagation speed (and thereby, elasticity) can be differentiated from the background. Although more appropriate for OCE, this definition is still incomplete because it is independent of image contrast. That is, an inclusion may be differentiable from the background, but an elasticity image based on the local wave speed for dynamic OCE in most cases does not present the true elastic modulus contrast between inclusion and background. Indeed, in many cases where bounded geometries are imaged, even relative wave speed differences (i.e., relative contrast) are influenced not only by shear elastic modulus but also by geometric factors such as medium thickness (e.g., in the cornea). 18,19 Indeed, it is not clear whether high spatial resolution maps of the mechanical wave speed represent faithful images of the true tissue elastic modulus.
In this study, we identify and resolve propagating waves from both sides of an interface to define spatial resolution for the simplest case encountered in dynamic OCE. The resolving power of the dynamic OCE system is quantified for the case of highly accurate contrast resolution (i.e., image contrast is quantitatively close to elastic modulus contrast) with minimal artifacts (i.e., all image features are closely related to true heterogeneity in the elastic modulus) around a single interface normal to a propagating plane wave. In application to dynamic OCE, we define resolution as the smallest possible inclusion that can be accurately detected without sacrificing quantitative contrast, which considers the probing mechanism (elastic wave propagation), detection mechanism (OCT), and the mechanical model.
In OCE, the probing mechanism determines the appropriate model used to estimate mechanical properties. In dynamic OCE, the typical assumption of an infinite, homogeneous, nearly incompressible elastic material would suggest that the elastographic resolution can be described by the OCT image resolution alone, yet, spatial resolution has been reported to range from the micron scale to multiple millimeters for similar reconstruction methods based on elastic wave propagation. [20][21][22][23][24][25][26][27][28] Recently, it was suggested that the mechanical model of the medium is directly related to resolution in the closely related field of compression OCE, 17 thus indicating that spatial resolution is not always defined by the resolution of the imaging method (OCT).
Here, we raise the main question of this paper-how do elastic wave characteristics and the detection capabilities of OCT influence the overall resolution in dynamic OCE? It has been previously demonstrated that the lateral resolution of acoustic radiation force-based ultrasound (US) elastography is defined by both excitation and detection transducers, i.e., the convolution of excitation and detection beam width profiles. 29 In dynamic OCE, the PSF of a typical OCT system is much narrower (usually tens of micrometers) than the PSF of the excitation source (from hundreds of micrometers to a few millimeters). Thus, the convolution between excitation and detection beam width profiles will largely resemble that of the excitation beam. If the excitation beam width truly determines spatial resolution in OCE, OCT cannot provide an advantage in resolution over conventional US elastography if we imagine that the same pushing source is used for elastic wave excitation. As we will show, spatial resolution in dynamic OCE is determined primarily by the properties of propagating mechanical waves and not by limitations of the typical OCT imaging system, indicating that highresolution dynamic OCE can only be achieved by spatial and temporal shortening of the excitation beam width.
To explore OCE resolution under highly controlled conditions, we used both detailed simulations and measurements on a specific model system. Specifically, propagating elastic waves with various spatial pulse widths (PW) were both numerically simulated and experimentally measured. The numerical simulation was designed to mimic experimental conditions in terms of wave excitation and model geometry, highlighting wave behavior close to an interface between two regions with different shear moduli. We primarily focused on Rayleigh waves propagating in the subsurface layer along an air-medium interface to explore scattering (reflection and transmission) effects on estimated local wave speeds near a vertical boundary (i.e., boundary perpendicular to the free surface) between the two regions. Experimentally, Rayleigh wave signals were generated within a tissue-mimicking phantom using acoustic microtapping (AμT) with a home-built, air-coupled ultrasound source 30 and detected by a home-built M-B mode, phase-sensitive, spectral-domain OCT (SD-OCT) system. Group velocity maps within multipart phantoms based on simulated and measured displacement fields were reconstructed and used to quantify lateral (direction perpendicular to the vertical boundary) resolution based on wave speeds near a boundary for varying OCE operating conditions.
As will be demonstrated below, spatial resolution is fully defined by the spatial PW of the mechanical wave, the spatial sampling used to track the mechanical wave (OCT-limited), and the assumptions of the mechanical model used to infer modulus. For sampling typically used in OCT, we show that the spatial resolution in dynamic OCE is primarily determined by a characteristic wavelength of the mechanical wave pulse and further complicated by complex wave-interface interactions.

Background
When a temporally short axial load is applied within a spatially infinite, isotropic, purely elastic solid, local energy is converted into a pulsed elastic wave (i.e., wave packet) propagating within the medium. The goal of all wave-based elastography systems in the elastic limit is to determine the bulk shear wave speed based on the displacement of a propagating elastic wave.
Starting with an infinite, homogeneous, nearly incompressible elastic material, consider a plane longitudinal (l-) or bulk shear (s-) wave propagating across the boundary (at plane x ¼ 0) between two isotropic elastic media characterized with Lamé constants λ 1;2 and μ 1;2 , correspondingly [ Fig. 1(a) Now consider plane waves propagating along the x-axis (normal to the z-axis), located within the bulk material. At the plane interface (x ¼ 0), propagating waves are both reflected and transmitted according to the reflection and transmission coefficients primarily determined by the ratio C ð2Þ l;s C ð1Þ l;s , 32 and the propagation speed is discontinuous at the boundary. Only plane waves propagating in the AEx direction are allowed under these conditions, and no other wave modes will be generated at the interface. The amplitude of counterpropagating wave packets is simply given by the transmission and reflection coefficients. These waves will not change their temporal shape; however, superposition of the counterpropagating plane waves in the x < 0 region close to the interface will affect the apparent shape of the incident wave packet. Preservation of temporal wave shape is exploited in traditional time-of-flight reconstruction to locally estimate propagation speed and, thereby, infer the elastic modulus. 33 In the spatial domain, the original waveform can stretch or compress as a function of local propagation speed. This scaling has been exploited to locally estimate propagation speed based on the wavenumber, thereby inferring elastic modulus. 24,34,35 To accurately reconstruct the propagating wave speed near the boundary for this simple case, the two counterpropagating modes must be differentiated to locally estimate the propagation speed of the wave packet traveling in the transmission direction. For shear-wave imaging in dynamic elastography systems, a directional filter 36 is usually applied to wavefield data gathered by the imaging system to extract one of the wave packets. If an ideal directional filter can be designed that does not distort the wave packet, then the spatial resolution of the elasticity map can approach the resolution of the primary imaging system (i.e., OCT-scale resolution for dynamic OCE). However, practical directional filters applied to data with finite bandwidth (BW) and finite signal-to-noise ratio (SNR) cannot preserve the spatial resolution of the primary imaging system even for this highly simplified case of plane wave propagation in an infinite medium. While this degradation is only amplified for more practical conditions encountered in dynamic OCE, it is not difficult to design a directional filter that improves reconstruction accuracy without limiting resolution to the finite BW of the filter. In other words, directional filtering will inevitably limit spatial resolution in the ideal case of wave packets traveling only in the AEx direction but does not serve as the limiting factor in typical dynamic OCE.
Typical dynamic OCE systems are confined to imaging depths of about 1 to 2 mm due to the penetration limits of OCT. Consequently, mechanical waves tracked in OCE are launched and propagate in highly bounded media, typically with a free surface (air) defining one of the boundaries. Many mode types are possible in tissue with such boundary conditions, namely, Rayleigh, Love, Stoneley, Scholte, or other types of guided waves. 5,32 Here, we consider the simplest case of boundary effects typically encountered in dynamic OCE [ Fig. 1(b)], where a temporally short axial load applied to the free surface (z ¼ 0) between air and the propagation medium (i.e., soft tissue) is used to generate a bulk l-and s-wave that travels into the medium, as well as a Rayleigh wave that travels along the free surface.
First, consider Rayleigh waves propagating along the x-axis far from the excitation source where they can be easily isolated from bulk waves. For a nearly incompressible material such as soft tissue, the Rayleigh wave speed is simply related to the shear-wave speed: 32 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 4 5 9 C R ≅ 0.955C s ; (2) and, hence, the shear elastic modulus can be easily computed from the Rayleigh-wave speed. 5 Assuming negligible viscosity, the temporal shape and amplitude of the propagating surface wave do not change with distance when viewed a sufficient distance from the excitation zone. This means simple time-of-flight measurements can accurately track propagation speed and directly convert the measurement to modulus. If this wave impinges on a vertical boundary similar to the bulk case shown in Fig. 1(a), then at first it would appear to represent the same problem. However, this is not the case. When a surface wave hits a vertical interface [ Fig. 1(b)], its energy is not fully distributed between reflected and transmitted Rayleigh waves; new modes are possible, 37 which greatly complicate the mechanical wave packet in this region.
To understand the effects of mode conversion in the neighborhood of the interface, the mechanical boundary conditions on strain (σ) must first be specified at (x ¼ 0, z > 0) and (x, z ¼ 0) interfaces E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 2 3 0 σ xz ð1;2Þ j z¼0 ¼ 0 σ zz ð1Þ j x¼0;z>0 ¼ σ zz ð2Þ j x¼0;z>0 : The boundary conditions must be applied to solutions of the general mechanical equations of motion in the medium E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 1 2 3 ρ ∂ 2Ũ ∂t 2 ¼ μ ð1;2Þ ΔŨ þ ½λ ð1;2Þ þ μ ð1;2Þ grad divŨ; (4) whereŨ is the displacement vector of the propagating wave. For the given boundary conditions, not only are transmitted and reflected Rayleigh waves produced, but also additional modes such as a Stoneley wave propagating along the z-axis and bulk shear waves traveling at an angle to the z-axis are launched. 37 The specific energies converted from the pure incident Rayleigh wave to secondary modes (longitudinal, shear, Stoneley, etc.) can be difficult to quantify and will vary greatly for different geometries and stiffnesses. 38 Because the values of each converted mode differ depending on factors such as the ratio μ 2 ∕μ 1 , it is difficult to generalize how each mode will specifically affect the reconstruction with depth for the wide range of conditions typical of OCE. Scattered body waves emerging in a particular direction create a finite region where, even for different mixed-mode amplitudes, the signal will be distorted. This region of undefined wave behavior, i.e., a mixture of multiple propagating modes, will corrupt Rayleigh wave tracking and, hence, its speed reconstruction. We also note that the Rayleigh wave waveform itself changes when passing the interface by transferring energy to other modes. As we show below, it is hard to separate the Rayleigh wave from other modes near the interface.
Finding an analytical solution for all modes near the x ¼ 0 interface is not a trivial problem, especially for broadband Rayleigh-wave pulses. Thus, to explore OCE spatial resolution under these conditions, we used both detailed simulations and experimental measurements on model systems. The specific geometry is designed to highlight wave-behavior close to an interface between two lateral regions with different shear moduli. Specifically, Rayleigh waves with varying pulse width were used to determine how close single-mode surface waves could approach a vertical interface without corruption due to scattering and mode-conversion. The ability to uniquely identify Rayleigh waves from both sides of the interface enabled quantitatively accurate reconstruction of the elastic modulus.
We present a detailed analysis based on the time-of-flight method to reconstruct changes in wave velocity around a vertical boundary and demonstrate the effects of mixed-mode propagation. Here we demonstrate that in a noiseless scenario, the Rayleigh wave pulse width primarily determines the spatial resolution of elasticity maps preserving the true elastic modulus contrast between media. Based on this analysis, we quantify lateral (x-direction for this model) resolution as the width of the transition zone around a vertical interface (x ¼ 0, z > 0) where the Rayleigh wave velocity cannot be accurately measured.
This lateral resolution is related to the minimal size of an inclusion that can be detected with accurate elastic modulus contrast.

Numerical Simulation
A finite element method (FEM) simulation was designed to provide an ideal, noise-free comparison between theoretical and experimental data. The elastic wave equation [Eq. (4)] was solved using OnScale (OnScale 2018, onscale.com, Redwood City, California). The domain consisted of two linearly elastic regions separated by a vertical plane (Fig. 2). Simulations were run with input bulk s-wave speeds (C s ) on either side of the vertical plane as C s ð1Þ ¼ 4.4 m∕s (19.36 kPa) and C s ð2Þ ¼ 2.5 m∕s (6.25 kPa), typical values for biological tissue. Density (ρ ¼ 1000 kg∕m 3 ) and l-wave speed (C l ¼ 158.2 m∕s) were uniform throughout the domain. These properties correspond roughly to a Poisson's ratio (v) of 0.4995. Numerical tests showed that above this longitudinal wave speed, solutions do not change appreciably (See Sec. 8, Appendix A).
Along the top boundary of the domain (z ¼ 0), we applied a spatially and temporally varying pressure load to simulate the acoustic push applied to the surface by the AμT transducer in our experiments. To define the spatial shape of the pressure load, we acquired a raster scan of the acoustic field [ Fig where R and C air are the reflection coefficient and l-wave speed of air, respectively (R ≈ 1, C air ¼ 340 m∕s), P 0 ¼ 5 kPa is the amplitude of the applied pressure, and d ¼ 600 μm is the push width, defined here as the FWHM of the acoustic intensity. The temporal push profile is given by a super-Gaussian function where T is the push duration, defined as the full-width-at-halfmaximum of the envelope function sðtÞ. A time delay (t 0 ) was introduced to avoid impulsive loading at time t ¼ 0. The push profile was centered, and full symmetry was assumed at the left boundary.
For a transient excitation that strongly satisfies the stress confinement condition (i.e., infinitesimally short push), the mechanical wave BW is determined by relaxation of the stress across the acoustic excitation beam, which occurs at the speed, C s : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 6 3 ; 3 8 9 BW ∝ C s d ; where d is the width of the push and defines the spatial width of the generated mechanical wave. If the push intensity is not confined during the time of stress relaxation [push duration (T) exceeds stress-relaxation time], BW is constrained by T.
The spatial pulse width (PW) of the Rayleigh wave will obey the inverse relationship in Eq. (7), expanding in space as a function of the temporal pulse duration. Limiting the BW of the pulsed wave will increase the spatial pulse width of the generated mechanical wave, regardless of what pushing method is chosen. For our experiment, the spatial profile IðxÞ was set to have an equivalent d ¼ 600 μm, and the push duration (T) was varied (250, 300, 450, and 545 μs) to reduce spectral content in the simulated mechanical wave.
The bottom and right boundaries were set to be absorbing layers to minimize reflections. The domain was discretized with linear FEMs on a regular rectangular grid with a minimum of 40 elements per wavelength (equivalent to a grid spacing of 2.75 μm). Because underintegrated linear elements are prone to spurious solutions, we applied Belytchko-Bindeman strain hourglass suppression. 39 The equations were integrated in time using an explicit time-stepping method to generate the full vibration velocity vector (time derivative of the dynamic displacement vector) at each space and time point. Only the vertical component of the vibration velocity was analyzed for direct comparison to experimental results (Fig. 4, Video 1).
It is worth noting that the numerical solution to Eq. (4) for the imposed boundary conditions includes wave energy propagating ahead of the Rayleigh-mode. This disturbance is evidence that the Rayleigh wave equation has additional solutions besides the primary Rayleigh root. One possible solution is a supershear evanescent wave coupled to a plane shear wave in the medium 40,41 and has been observed experimentally in multiple studies. [42][43][44] For this analysis, we treat the primary Rayleigh root solution as the incident wave-mode and ignore incident wave energy leading the primary disturbance. We also note that this "fast-mode" is generated due to mode-conversion at the vertical interface and serves as a further disruptor when tracking the Rayleigh wave speed in the near-field region close to a boundary.

Tissue Mimicking Phantom
Boundary effects were explored experimentally using an elastic tissue-mimicking material [polyvinyl alcohol (PVA)-based] with controllable mechanical properties. PVA-based phantoms were created using protocols adapted from Ref. 45. They were made by first mixing a 4∶1 ratio of dimethyl sulfoxide (DMSO, CAS: 67-68-5, EMD Millipore Corp. Billerica, Massachusetts) to water using a stir plate. The optical properties were tuned by introducing titanium dioxide microparticles. A concentrated stock solution of 0.3 wt. % titanium fine powder suspended in water was prepared and a tip sonicator (Digital Sonifier 450, Branson, Danbury, Connecticut) was used for 3 min (30% duty cycle) at an amplitude of 30% to help disperse the microparticles in solution. The concentrated microparticle solution was then added to the DMSO and water solution to achieve a concentration of 0.025 wt. %. 4 or 6 wt. % PVA (146 to 186 kDa, >99% hydrolyzed, CAS: 9002-89-5, Sigma-Aldrich Corp., St. Louis, Missouri) was then added to the solution. The solution was covered and stirred on a hot plate maintaining a temperature of 120°C for ∼1 h to dissolve the PVA. Once fully dissolved, it was degassed using a vacuum chamber to remove any bubbles before casting in molds and stored at −20°C until the phantom hardened. The entire mold was 10 cm × 10 cm × 10 cm.
Phantoms with dissimilar mechanical properties were prepared using this protocol but with varying concentrations of PVA (4% and 6%). The 4% PVA solution was allowed to cool Journal of Biomedical Optics 096006-5 September 2019 • Vol. 24 (9) down to 35°C in a mold with a central barrier prior to casting the 6% PVA solution. The microparticle concentration was slightly increased in the 4% PVA solution to visualize separate regions. The mold barrier was removed and the 6% PVA solution was cast into the pre-existing phantom, resulting in a 10 cm × 20 cm × 10 cm welded-surface, two-part material (Fig. 5). Once cast, the mold was stored at −20°C for up to 12 h to solidify. Finally, the casted, hardened phantom was placed in a water bath. During this process, the DMSO slowly diffused out of the phantom in exchange for water. The water bath was regularly changed for a minimum of 48 h to remove all DMSO. The phantoms were stored in deionized water to prevent phantom dehydration.

Optical Coherence Tomography Wave Detection
To generate a planar Rayleigh wave, a cylindrically focused, aircoupled acoustic radiation force (AμT) was directed at the phantom material to induce a localized displacement. 30,46 The AμT transducer effectively applied a successive number of localized line "pushes" perpendicular to the sample surface over a wide  Journal of Biomedical Optics 096006-6 September 2019 • Vol. 24 (9) region relative to the propagation distance of interest, resulting in an elastic wave that can be approximated as a plane wave (reducing geometric diffraction effects at distances shorter than the AμT excitation line length). The full-width-at-halfmaximum lateral focus of the ultrasound transducer centered at the axial focus was measured as ∼0.35 mm. 30 In practice, the ultrasound transducer was fixed at an angle (∼45 deg to vertical) above the phantom to avoid blocking the probe beam from the OCT system, which resulted in ∼0.6 mm push width, equal to that of the numerical simulation. The ultrasound excitation was focused on the sample surface 5 mm in the −x direction from the initial OCT observation point to allow Rayleigh waves to fully separate from simultaneously generated bulk l-and s-waves. Four pulse durations (100, 200, 400, and 600 μs) were used in the experiment.
Note that for AμT (same as for laser excitation of ultrasound), 47 when a surface wave signal is excited by an infinitesimally short (in time) push, the spectral characteristics of the propagating signal are defined by the spatial width of the push. When the push duration is not infinitesimally short in time, signal spectral characteristics will be determined by the convolution of the spatial and temporal push functions. In the experiment, we define the spatial function as ∼0.6 mm wide, as determined by the transducer field and its tilt angle to the sample. To make experimental propagating signal shapes and spectra match those of the simulations, the push duration can be adjusted. Thus, simulation parameters were tuned slightly so that experimental and simulated propagating signals had close spectral and temporal shapes.
Elastic waves were tracked in the phantom material using a phase-sensitive OCT (PhS-OCT) system designed to precisely capture elastic wave displacements in both space and time. For this experiment, M-B mode OCT was used to track wave propagation. The M-B mode PhS-OCT system employed an SD-OCT setup, which has been described in a previous study. 48 Briefly, it consisted of a broadband super-luminescent diode (1310 nm center wavelength, ∼46 nm spectral BW), a beam splitter, a stationary reference arm, a sample arm integrated with a set of galvo-scanners, and a high-speed spectrometer. The A-line rate of the system (temporal resolution), determined by the sampling rate of the 1024-pixel line-scan InGaAs array, was 46.5 kHz, sufficient to acquire induced mechanical waves inside the phantom without aliasing. The optical resolution was ∼15 μm axially and ∼24 μm laterally.
To track Rayleigh wave propagation on the sample surface, an external TTL trigger was used to synchronize the PhS- where λ was the center wavelength of the broadband light source,n was the refractive index of the medium, and f s was the sampling frequency.

Data Processing
Multiple raw vibration velocity field measurements were acquired and averaged to suppress Gaussian-distributed system phase noise. The number of repeats for each push duration (100, 200, 400, 600 μs) were 20, 16, 12, and 8, respectively so that the average SNR in each data set was approximately the same (50 dB). SNR was quantified using the noise floor [Amp noise , defined as the maximum amplitude of the phase variation along signal-averaged M-mode sequences with no wave propagation (static)], relative to the maximum amplitude of the vibration speed from the incident Rayleigh wave, Amp signal , in each signal-averaged data set, and quantified using SNR ¼ 20 log 10 Amp signal Amp noise . Once signal-averaged displacement sets were acquired, both simulation and experimental data underwent identical preprocessing that included a directional filter to suppress wave reflection at the boundary by masking the second and fourth quadrants of the two-dimensional (2-D) fast Fourier transform (FFT) spectrum. 50 Wave amplitude was normalized at each lateral location and a temporal bandpass filter was applied to remove system noise and side-lobe excitation components using a Gaussian window centered at the peak frequency and bandlimited at the −20 dB cutoff.
Typically, in time-of-flight reconstruction, the maximum of the normalized cross-correlation (NCC) function is used to determine the time delay (Δt) between a reference and delayed waveform spaced by Δx to estimate wave velocity. For experimental data sets corrupted by noise, jitter error will reduce the accuracy of the time lag associated with the maximum magnitude of the NCC. 51 By assuming sample homogeneity over spatial regions of interest, reconstruction kernels can be expanded over small spatial domains to generate more accurate speed estimates. However, expanding the reconstruction kernel size will reduce lateral resolution around any inhomegeneity. 52 To optimize the tradeoff between resolution and accuracy, a multiplekernel-weighted-averaging approach (MKWA) was applied to improve reconstruction accuracy and minimize the loss in spatial resolution.
The MKWA method can reduce jitter error by calculating time-delays for all available kernel sizes [bound by x AE k∕2, where the detector spacing (Δx) in each case was denoted as the kernel size, k] centered at all locations. The local velocity for each kernel spacing was determined using the NCC function for waveforms assuming elastic wave propagation parallel to the surface. Quadratic interpolation of the lag term about the max was employed to refine the lag-time measurement. The wave speed at all locations for an overdetermined list of velocities is expressed as the sum of the product of wave speed measurements using different kernel sizes and their corresponding weights E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 3 2 6 ; 1 5 8 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 1 0 1 ; i∈ ½1; 2; : : : n; where C R i ðx; zÞ represents the wave speed at location ðx; zÞ reconstructed using the i'th kernel. Here, w i denotes the weight Journal of Biomedical Optics 096006-7 September 2019 • Vol. 24 (9) for each kernel spacing, and n is the number of kernels used for wave speed reconstruction. The linear weighting function assigned higher weights to smaller kernels [Eq. (10)], which provided an ∼8 dB improvement in time-delay reconstruction accuracy while maintaining nearly full resolution for an equivalent kernel size of one-tenth ( 1 10 ) the spatial pulse width of the Rayleigh wave in the soft region (detailed in Sec. 9, Appendix B). For consistency, noisefree simulation data sets were reconstructed using the MKWA method with the same equivalent kernel size.
A one-dimensional edge profile was generated using a depthaveraged set of MKWA wave speeds for both simulation and experiment. To minimize the influence of the nonvertical interface, the boundary between the two different media for each depth was aligned before performing depth averaging. The two-dimensional velocity field was centered at each z-location, and the mean taken over a depth equal to one-third of the minimum measured pulse width.
Antisymmetric sigmoid fitting was applied to each depthaveraged edge profile to quantify the transition zone between media. A different fitting function was employed to define the transition region on either side of the x ¼ 0 boundary, where each transition profile was combined with its centrosymmetric copy to form a single-part edge profile for sigmoid fitting. The sigmoid function on either side of the boundary is expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 4 7 7 C R fit ¼ Half of the −6-dB cut-off of the first-order derivative of C R fit with respect to x was used to quantify the transition size on both sides of the vertical boundary. The total transition zone was equal to the sum of the measured transition on either side of the boundary. Similar sigmoid-based fits have been used to describe elastic resolution in other studies. 17,22,29,[52][53][54] 5 Results

Medium with Vertical Interface
The spatial pulse width of the Rayleigh wave in experimental and simulated data is fully defined in the hard (PW 1 ) and soft (PW 2 ) regions as the full-width-at-half-maximum of the signal envelope [Figs. 6(a) and 6(b)]. The spectral content is defined by the center frequency (normalized first moment of the spectral amplitude bound by −20-dB), 55 F c , and the maximum frequency (upper frequency at −6-dB), F max , of the power spectrum following an FFT of the wave disturbance far from any source or boundary [ Fig. 6(c)]. As wave speed is constrained by the material properties of the medium, wavelengths are long for low-frequency excitation. Conversely, the spatial pulse width is shorter for broadband (high-frequency) excitation.
In a linear elastic material, the frequency of a harmonic wave must be conserved as it passes through the interface between two media. As a pulsed wave is simply the superposition of signals over a wide range of frequencies, the pulsed signal spectra will also be conserved, and the pulse width will scale in space as a function of the wave speed. The measured pulse width in the soft portion approximately (due to mode conversions at the interface) scales linearly with the change in wave speed C s ð1Þ ≅ ð 2.5 m∕s 4.4 m∕s Þ · PW 1 ) for all data sets. To simplify the analysis of resolution at the interface, here we introduce a mean pulse width ðPW mean Þ as PW mean ¼ 1 2 ðPW 1 þ PW 2 Þ, defining a characteristic wavelength of the Rayleigh wave in the region close to the interface. As we will show below, this characteristic wavelength better captures the mechanical resolution of OCE than the optical resolution of the OCT imaging system. Two-dimensional group velocity maps reconstructed from simulated and experimental data with different pulse widths are shown in Fig. 7  As demonstrated, mode-converted waves on the order of the pulse width affect the estimated wave speed close to the boundary, resulting in a finite transition size before stabilizing to the true Rayleigh wave-speed. The mode-converted bulk shear and Stoneley waves traveling into the medium degrade reconstruction accuracy, manifesting as features in the group velocity map that are not reflective of the actual geometry (i.e., artifact). The directional filter can mitigate this effect in the negative-x direction (−x) but is still insufficient to remove scattered wave-packets that travel at an angle normal to the z ¼ 0 surface, particularly in the positive x direction (þx). In addition, the temporal wave-shape change due to boundary traction forces will limit the NCC from accurately reconstructing a propagation speed directly on either side of the boundary. The structural OCT image confirmed that the vertical boundary in the phantom was normal to the surface (<1 deg) to ensure that no unexpected wave modes were generated. The depth-averaged group velocity measurements and the corresponding sigmoid-fit profiles are normalized and shown in Fig. 8. The transition size was defined as the full-width-athalf-maximum of the first derivative of the sigmoid fit, with respect to the lateral direction of the fitting function. Note that the small "bump" appearing on the þx side of the boundary right after the transition zone [see Fig. 8(a)] is an artifact resulting from averaging the 2-D velocity reconstruction at a depth relative to the wavelength. Because the reconstructed group velocities at each depth were centered before averaging, and depth averaging was performed relative to the elastic pulse width, any artifact in the reconstruction will affect each data set in the same way. Aligning the boundary prior to depth averaging also ensures that transition zone "blurring" was not a result of depth averaging.
Results show that the total size of the transition zone region scaled linearly as a function of incident and transmitted spatial pulse widths (Fig. 9). The spatial width of the pulsed wave changed as a function of the local velocity, thus the quantification based on the first derivative of the asymmetric fitting function determined the transition size relative to the pulse width uniquely in the hard (−x) and soft (þx) portion of the model. When a single Rayleigh wave-mode was incident on a vertical boundary, experimental and simulated results demonstrated a linear relationship between the pulse width and the size of the blurred region, where the transition out of, or into, the other material occurs in space at approximately 1 4 the size of the spatial pulse width in each material. This result indicates that to accurately differentiate between two regions with different moduli, a finite region equal to at least 1 2 the mean pulse width, ½ðPW mean ¼ 1 2 ðPW 1 þ PW 2 Þ must be considered. Thus, the ultimate resolution in dynamic OCE may be approximated by 1 2 PW mean .

More Complex Media Boundaries
To examine the relationship between resolution and contrast for more complex internal boundaries, the wave field generated by a Rayleigh wave passing through a three-part, horizontally layered material was simulated. The two outer layers were assigned shear wave speeds of C s ¼ 2.5 m∕s, whereas the inner layer, representing a stiff inclusion, was assigned a shear wave speed of C s ¼ 5 m∕s. The simulated push excitation was applied in the first outer layer, generating a Rayleigh wave with a spatial pulse width of 0.5 mm. This corresponds to a maximum potential spatial pulse width of 1 mm in the hard inclusion. For the wave speeds in this arrangement (2.5 and 5 m∕s), considering that frequency is conserved, the mean pulse width is PW mean ¼ 0.75 mm. We varied the width of this middle layer (1, 0.5, 0.25, Journal of Biomedical Optics 096006-9 September 2019 • Vol. 24 (9) and 0.1 mm) to examine whether group velocity reconstruction can accurately resolve the inclusion. The reconstruction presented the inclusion at nearly full contrast with minimal artifact [ Fig. 10(a) and 10(b)] when the inclusion width is greater than the mean pulse width. The contrast begins to deviate proportionally when the inclusion size reduced relative to PW mean {as seen for the cases of 0.5, 0.25, and 0.1 mm inclusion widths [Figs. 10(c)-10(h)]}. At an inclusion size equal to one-tenth the PW mean , [Figs. 10(g) and 10(h)], the reconstruction detected the vertical inclusion but with very poor contrast. We conclude that the minimum inclusion size that can be reconstructed with full contrast is approximately equal to the mean pulse width PW mean . Note here that the minimum inclusion size is twice the transition zone between two media because the inclusion consists of two vertical interfaces, with a transition region of 1 2 PW mean at each. A corollary of the previous analysis is that full resolution at maximum contrast may be achieved for smaller inclusion sizes using shorter pulse width (i.e., higher BW) waves.
As more complex geometries are introduced, boundary conditions become even more disruptive to reconstruction accuracy. Due to the elliptical nature of particle motion in a Rayleigh wave, local stresses transfer energy along both lateral and axial directions. These multidimensional traction forces can generate complicated wave fields at all interfaces. Once a single-mode wave interacts with an inclusion boundary, wave components are guided along the boundary, mode-convert, diffract, and interfere with the original wave shape, reducing correlation and negatively affecting spatial resolution in complex boundary conditions.
To illustrate this point, a 5-m∕s semicircular inclusion (6 mm wide and 3 mm deep) was simulated at the center of a 2.5-m∕s media [ Fig. 11(a)] and probed using a Rayleigh wave with mean pulse width, PW mean , of 0.75 and 4.5 mm, respectively (see Video 2). The boundaries of the semicircular inclusion are better defined when reconstructed using PW mean ¼ 0.75 mm [ Fig. 11(b)] compared with PW mean ¼ 4.5 mm [Fig. 11(c)]. Note that boundary regions are corrupted in the reconstruction even for PW mean ¼ 0.75 mm due to mode conversions. For longer pulse widths (lower frequency) waves, the reconstructed image identified the inclusion but could not accurately quantify modulus, highlighting the importance in defining both image contrast and resolution.  ) is also shown for both experiment and simulation.

Discussion
The typical dynamic OCE procedure is to: (1) launch a broadband mechanical wave within a sample using a temporally short, spatially sharp excitation force (akin to shear wave elasticity imaging, 56 (3) use local displacement (vibration velocity) fields to estimate mechanical properties. 5,21,62 Under highly controlled conditions, the estimated wave velocity can define the shear modulus (stiffness) within a small region. 5,20,56,61,63,64 Even with typical limitations imposed by a practical measurement system (i.e., measurement noise will introduce additional imperfections), we have shown using numerical simulations and experimental studies in tissue-mimicking phantoms that OCE spatial resolution is fundamentally limited by the properties of propagating mechanical waves. As shown above, the fundamental resolution limit in dynamic OCE is determined by the pulse width of the elastic wave used to probe a material. As complicated boundaries can change wave shape and generate interfering wave-modes, a quantitative, artifact-free time-of-flight reconstruction cannot be done until the parent mode has separated from the complex displacement field. It follows that the local material properties will play a role in the achievable resolution. Thus, resolution cannot be defined simply by the detection mechanism of OCT. It is interesting that, although Qian et al. 29 reached similar conclusions and correctly defined resolution in conventional US elastography, they incorrectly reported that the effective lateral resolution in dynamic OCE is defined by the capabilities of the OCT imaging system. As shown in this work (as well as by Hepburn et al. 17 ), the resolution in both wave-based and compression OCE is not defined solely by the detection mechanism of OCT.
If a Rayleigh wave encounters an inclusion with a modulus different from that of the original propagation medium, the pulse width of the incident wave, and any generated wave-modes, will scale as a function of the relative differences in C s , and the resolution will be defined by the local pulse width. As we have shown, the minimum resolvable feature (e.g., any interface with elastic modulus differential) in a system will be approximately equal to ½ the mean spatial pulse width at all locations within the imaging field. The minimum inclusion size reconstructed at  full contrast is twice this limit and corresponds approximately to the mean pulse width at the inclusion interface. Note that when the inclusion modulus is unknown, the resolution in OCE can be defined only approximately because the mean pulse width depends on both the inclusion and surrounding medium moduli. Defining the resolution by the probe pulse width is not accurate and can either underestimate or overestimate it depending on inclusion properties.
Even with noise-free conditions near a boundary, mode-conversion will reduce correlation between reference and delayed wave forms and lead to errors in wave speed estimates. For an elastic plane wave traveling in a known direction within a nearly incompressible homogeneous material bound by a single surface, the maximum of the NCC between any reference and delayed signal will be very close to the theoretical limit of one at all locations, enabling direct quantification of the shear modulus. 65 When this wave approaches a heterogeneity, the discontinuity imposes boundary conditions [Eq. (3)] producing multiple mode conversions near the interface (Fig. 4, Video 1) changing the shape of the transmitted Rayleigh mode. When multiple modes contribute to the displacement field, errors due to decorrelation will shift the maximum, producing an inaccurate measurement of local time-delay.
As bulk shear and Stoneley waves launched at an interface travel in different directions from the parent Rayleigh wave packet, mode mixing will be spatially dependent, producing an artifact close to the interface in the wave-speed reconstruction (Figs. 7, 10, and 11). This effect has been recognized as an "edge-effect" in MRE, 66 and modulus quantification is often only considered accurate when viewed approximately ½ of a pulse width from any boundary or source. 67,68 Thus, image details with spatial variation less than about a half of the elastic pulse width cannot be used to directly determine local elastic modulus assuming a single-propagating wave-mode.
System noise, relative to wave amplitude (SNR), will further reduce reconstruction resolution and accuracy. In practice, phase-stability of the OCT system determines the minimum detectable displacement (vibration velocity) based on the time-delay between reference and delayed signals. For a homogeneous material, noise will shift the maximum of the cross-correlation function away from the actual peak, referred to as jitter error. The predicted jitter magnitude is the minimum error achievable by any unbiased time delay estimation algorithm. As SNR decreases, jitter will increase. It has been shown that time-delay estimation error is worse for very short window sizes and can be improved with higher SNR. 51,52 Assuming single-mode propagation, detector spacing must be wide relative to the mechanical wavelength to minimize jitter using a single kernel window size. For an inclusion, the effective detector spacing to accurately reconstruct time delays determines OCE resolution. The lateral resolution of OCT (determined by the numerical aperture of the objective lens in the imaging system, typically tens of micrometers) can provide spatial sampling much finer than typical acoustic wavelengths, 69 allowing an overdetermined set of time-delay terms to greatly improve the accuracy of local wave-speed measurements, 70 a feature used here to improve reconstruction accuracy. While current systems can detect nm-scale displacements, [71][72][73] both OCE resolution and reconstruction accuracy may be optimized by sub-nm-scale displacement sensitivity, allowing accurate timedelay reconstruction of low amplitude waves with very small detector spacing.
Another way to improve resolution is to utilize more complicated reconstruction methods to quantify material stiffness. To accurately recover spatial resolution to a fraction of the pulse width for an arbitrary imaging environment, full wave field reconstruction should be performed at a fine spatial scale similar to current coarse spatial scale methods in magnetic resonance elastography (MRE). 66,[74][75][76][77] The general case of wave field reconstruction often assumes local homogeneity, a valid assumption when probing regions spaced by approximately ½ the wavelength of a harmonic elastic wave. Consequently, extrinsic variables can significantly influence the spatial resolution as defined in this study.
Directional filtering can reduce the influence of reflected and scattered wave-packets to improve reconstruction accuracy by attempting to isolate a single wave packet traveling in an assumed direction. 36 However, wave scattering can make it difficult to isolate a Rayleigh wave, particularly when there are mixed-mode projections traveling in the same direction as the parent wave. Practical directional filters applied to data with finite BW and finite SNR cannot preserve the spatial resolution of the primary imaging system. However, practical directional filters can be designed such that the BW and SNR preserve the parent wave pulse width, reducing artifact without limiting spatial resolution. As a thought experiment, if a filtering approach could be designed to isolate and remove independent wavemodes, it may be possible for OCE to approach OCT resolution. Of course, this requires that all possible modes are known.

Conclusions
In conclusion, the spatial resolution in dynamic OCE when tracking mechanical wave propagation near an interface is ultimately limited by the characteristic wavelength of the mechanical wave or by the pulse width for broadband signals. The ability to accurately reconstruct an inclusion with contrast nearing that of the physical modulus breaks down when the mean pulse width is greater than the size of the physical inclusion. While the ultimate resolution in dynamic OCE struggles to match that of its parent imaging modality, OCT, an advantage that cannot be understated is that high-resolution structural information can be acquired simultaneously. Thus, OCE can provide unique information on both mechanical modulus and tissue light scattering at a scale and speed difficult to achieve with other imaging modalities. It follows that OCE and OCT can be combined with optical microangiography 78 to potentially provide quantitative image contrast based on tissue elasticity, optical properties, and microvasculature in a single acquisition that can be done within seconds. The combination of relatively coarse quantitative elastic modulus mapping with the high-resolution capabilities of OCT and OCT angiography may provide valuable clinical information that will likely continue to drive the future development of this technique. Soft biological tissues typically exhibit a longitudinal wave speed (≈1540 m∕s) orders of magnitude higher than the shear wave speed (1 to 10 m∕s). This property presents a number of computational challenges including long simulation times, increased floating point error, and the formation of spurious solutions due to hourglassing 39 (i.e., artificial, zero-energy strains). In practice, it is not necessary to directly model the true longitudinal wave speed, as numerical solutions will tend to converge in the limit of C s ≪ C l when the longitudinal wave speed is sufficiently high. To verify this, we compared numerical solutions for different longitudinal wave speeds from 12.2 to 353.6 m∕s with the theoretical solution predicted by the elastodynamic Green's function convolved with the spatial and temporal push profiles. 79,80 In all cases, the shear wave speed was set as 2.5 m∕s. Far-field waveforms show close agreement to the theoretical solution, particularly when the longitudinal wave speed exceeds 158.2 m∕s (Poisson's ratio >0.4995) (Fig. 12). The difference between numerical and exact solutions is quantified using the l 2 error. Results indicate that it is not necessary to set the true value of C l ¼ 1540 m∕s (typical for soft tissues) in the simulations. Any C l corresponding to Poisson's ratio larger than 0.4995 will produce the same wave propagation over the time scales used in these simulations. Using this fact, simulation times can be dramatically reduced, spurious solutions due to hour glassing can be virtually eliminated, and floating point error in the numerical solution can be minimized.

Appendix B. Multiple-Kernel-Weighted-Averaging
Due to jitter error in time-of-flight reconstruction, it is practical to increase detector spacing to improve the accuracy of the measured wave speed. However, increasing the kernel size will reduce spatial resolution. To accurately reconstruct wave speed in experimental data, while minimizing spatial resolution loss, the MKWA method described in Sec. 4.3 was employed. The trade-offs between resolution and accuracy in MKWA and single kernel methods are explored below. For direct comparison between MKWA and the single kernel method, an equivalent kernel size based on the weighting function [Eq. (10)] was defined as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 2 ; 3 2 6 ; 5 3 3 where w i represents the weight for the i'th kernel, k i is the spacing of the i'th kernel, and n is the number of kernels used. Twodimensional group velocity maps were calculated in experimental and simulated data sets using single and equivalent kernel sizes ranging from the minimum OCT scan interval to the maximum spatial pulse width of the elastic wave. Reconstruction accuracy as a function of increasing weighted averages [increasing n in Eqs. (9) and (10)] was compared to the single kernel time-of-flight reconstruction using the experimental data set corresponding to PW mean ≅ 1.54 mm [ Fig. 6(a), Video 1]. Assuming a quasi-semi-infinite homogeneous phantom (i.e., forward propagating Rayleigh wave only), with wave speed values generated using an unbiased estimator, variation in the local wave velocity, C R , can be used to define reconstruction accuracy as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 3 ; 3 2 6 ; 3 3 2 where C R is the average group velocity in a defined region, and σðC R Þ represents the corresponding variance of the group velocity. For this analysis, SNR C R was calculated in the experimentally equivalent C s ð1Þ and C s ð2Þ regions for all regions in z, excluding velocity values within one pulse width of the x ¼ 0 boundary. The measured SNR C R from equivalent kernel values is summarized in Table 1.
Increasing the kernel size to compute spatial gradients along the transverse dimension improves accuracy at the expense of OCE spatial resolution. The loss (%) in OCE resolution as a function of increasing kernel size for both single kernel and MKWA methods was quantified using the antisymmetric sigmoid fitting function described in Sec. 4.3 for the equivalent noise-free simulated data set. The minimum kernel size, equal to the grid spacing (2.75 μm), is representative of super-resolution OCT. 81,82 The minimum kernel size will not reduce resolution and was used as a baseline for determining any relative reduction in resolution with increasing k. Spatial resolution losses for all equivalent kernel sizes are summarized in Table 1. The tradeoff between SNR C R gain and spatial resolution loss using MKWA versus single kernel processing is summarized in Table 1.
For a kernel size 1 10 th the minimum pulse width, MKWA provides an 8.1 dB increase in reconstruction accuracy and a negligible increase in the measured transition zone. Conversely, the single kernel method reduced spatial resolution at the equivalent kernel size and is more affected by noise-induced jitter error. Based on this analysis, all wave speed images constructed from experimental data used MKWA processing with an equivalent kernel size of 1 10 th the minimum pulse width.

Disclosures
The authors declare they have no competing financial interests.  (9) optical coherence tomography with surface acoustic wave elastography and shear wave elastography for clinical application.
John J. Pitre Jr. received his BSE in biomedical engineering from Tulane University and MS and PhD degrees in biomedical engineering from the University of Michigan. He is currently a postdoctoral researcher at the University of Washington in the Department of Bioengineering. His research interests include ultrasound and optical coherence elastography, tissue biomechanics, and elastic wave propagation.
Liang Gao graduated from the College of Optical Sciences, University of Arizona. His PhD thesis was about ultrasound elasticity imaging of the human posterior tibial tendon. As a senior fellow in Department of Bioengineering, University of Washington, his research focuses on medical ultrasound imaging, shearwave elastography, and optical coherence elastography.
David Li received his BSE, MSE, and PhD degrees in biomedical engineering from the University of Michigan. He is currently a postdoctoral researcher at the University of Washington in the departments of bioengineering and chemical engineering. His current research areas are in photoacoustic imaging, ultrasound-based contrast agent design, and nonlinear imaging modalities blending the use of optics and acoustics.
Ivan Pelivanov is a researcher and associate professor of the Department of Bioengineering at the University of Washington. His projects include areas of ultrasonics, laser physics, laser-ultrasonics, and photoacoustics, both fundamental and applied, in NDT and biomedical fields, where light and sound could be used for diagnostics and imaging. He has spent years in designing ultra-broadband contact and noncontact transducers, from single element probes of different sizes and geometries to multi-element arrays.
Shaozhen Song received his doctoral degree from the University of Dundee, Scotland, United Kingdom. He is currently a postdoctoral fellow in the Biophotonics and Imaging Laboratory (BAIL), University of Washington, Seattle, USA, a biophotonics research group directed by Professor Ruikang Wang. His current research interests include biophotonics and imaging, especially in swept source optical coherence tomography (OCT), optical microangiography, optical elastography, as well as their applications in ophthalmology, dermatology, and cancer.
Chunhui Li is currently the programme director of biomedical engineering and a lecturer of mechanical engineering at the School of Science and Engineering, University of Dundee. Her main research interests are in photonics and ultrasound systems in biomedical applications. Her team is currently focusing on functional optical imaging using coherence gating (OCT) development for tissue stiffness and microvascular mapping.
Zhihong Huang is a professor of biomedical engineering at the University of Dundee. She received her PhD in mechanical engineering from the University of Glasgow. Her research interests include ultrasound and photonics instrumentation for ultrasonic surgical tool design and characterization, sono-and optical-elastography for medical diagnosis, and image guided and robotic assisted interventions. She is a fellow of IET.