Ultra-wideband optical coherence elastography from acoustic to ultrasonic frequencies

Visualizing viscoelastic waves in materials and tissues through noninvasive imaging is valuable for analyzing their mechanical properties and detecting internal anomalies. However, traditional elastography techniques have been limited by a maximum wave frequency below 1-10 kHz, which hampers temporal and spatial resolution. Here, we introduce an optical coherence elastography technique that overcomes the limitation by extending the frequency range to MHz. Our system can measure the stiffness of hard materials including bones and extract viscoelastic shear moduli for polymers and hydrogels in conventionally inaccessible ranges between 100 Hz and 1 MHz. The dispersion of Rayleigh surface waves across the ultrawide band allowed us to profile depth-dependent shear modulus in cartilages ex vivo and human skin in vivo with sub-mm anatomical resolution. This technique holds immense potential as a noninvasive measurement tool for material sciences, tissue engineering, and medical diagnostics.


INTRODUCTION
Measurement of the mechanical properties of materials is routinely performed in many places in sciences, engineering and industries, as well as hospitals [1][2][3][4] .Widely used tools include strain-stress testing 5 and dynamic mechanical analysis 6 for measuring bulk properties, and atomic force microscopy (AFM) 7,8 and microrheology 9,10 for local measurements.In clinical medicine, elastography has been adopted for disease diagnosis [11][12][13][14] .Elastography allows to noninvasively measure the elasticity of tissues in normal and abnormal states using medical imaging modalities, such as ultrasound and magnetic resonance imaging (MRI).In all these tools, samples under test are deformed by some force, their responses are measured, and the mechanical properties are calculated from the data 15,16 .The timescale or frequency range of the measurement spans from quasi static (as slow as 10 -5 Hz) to acoustic (as fast as 10 3 Hz) ranges.Higher speeds up to a few tens of kHz have been used in AFM 17 and in our recent work on dynamic optical coherence elastography (OCE) 18 .In the latter, elastic waves of short impulse or continuous-wave monotones are generated, their propagation within a tissue is visualized, and from the data the wave velocities are determined and related to the elastic moduli of the tissue 15,19 .
While the material analysis thus far has been focused on the quasi-static to acoustic ranges, we hypothesized that a higher frequency range beyond 10 kHz can offer a window of opportunity that was previously underappreciated especially for elastography.First, the higher frequency data can reveal the viscoelastic characteristics of materials in the shorter time scale.Dynamic mechanical analysis (DMA) is widely used to characterize viscoelasticity, but it has a limited frequency range (1 to 100 Hz).The timetemperature superposition technique can mimic the high-frequency responses of simple, homogeneous materials 20 but is not applicable to composite materials and living matters 9 .Ultrasound non-destructive testing (NDT) 21,22 is an established technique using ultrasound waves, typically at 500 kHz to 20 MHz, to measure the time of flight across a material or structure with well-defined boundaries to determine bulk elastic properties, but this technique is not directly applicable to elastography and especially for complex materials such as tissues.Second, since the spatial resolution of elastography is approximately given by the wavelength of elastic waves 23 , the higher frequency can lead to the higher resolution.Third, wave velocity dispersion over an extended frequency range is useful to extract more detailed mechanical information such as the depth-dependent variation of elasticity and internal stress.The depth information is not accessible using laser Doppler vibrometry or AFM-based rheology that only probes the surface [24][25][26] .Despite all these anticipated benefits mentioned above, it is technically challenging to extend the upper frequency limit.At the same amplitude of vibration, the energy dissipation or damping of a wave increases with ~3 ( , frequency) because the wave energy is proportional to  2 and the viscosity tends to grow with .So, with increasing frequency, the wave amplitude must be reduced by ~3/2 to avoid excessive sample heating or damage.Detecting the reduced amplitude requires improved sensitivity, for example, by a factor of a thousand when comparing 100 kHz to 1 kHz.To analyze waves with submicron amplitudes, a nanometer-scale sensitivity is required to an instrument.OCE is an emerging elastography technique 11,27 built on optical coherence tomography (OCT) 28 .OCE has been developed with a variety of different system architectures and applied to various tissue types, such as cornea 18,29 , sclera 30 , breast 31 , brain 32 , and skin 33 .As OCE uses optical interferometry, it offers superior detection sensitivity in the order of nanometers, compared to micrometers for ultrasound and millimeters for MRI.The high sensitivity of OCE make it a good candidate for elastography beyond the acoustic range.
Here, we demonstrate an OCE system capable of visualizing elastic waves from the acoustic to ultrasonic frequencies up to a few MHz.We developed a novel aliasing technique, demodulation algorithm, and jitter-correction method to handle such high frequencies far beyond the typical axial line-scan (A-line) rates of OCT.After verifying the system with homogeneous materials, we apply it to dynamic shear analysis of soft materials, and stiffness mapping of complex tissues in a knee joint ex vivo and human skin in vivo with high, 3-dimensional resolution.Our results demonstrate the benefits and broad applicability of ultrawideband OCE.

Signal demodulation of high-frequency vibration in swept-source OCE
Figure 1a depicts a basic interferometer used in swept source OCT, where the optical wavenumber is tuned repeatedly with a period of :  =  0 +  1 [], where [] = (, ) − /2,  0 = 2/ (, center wavelength), and  1 is a tuning rate 34 .Imagine a mirror-like sample that vibrates with an amplitude , angular frequency   (= 2  ), and reflection coefficient  at a mean depth of  0 .As the depth is modulated by  cos (  ) , we obtain an interferometric signal () =  cos (2 0  0 + 2 1  0  + 2 0  cos(  ) + 2 1  cos(  )), where 2 1  0  and    are assumed to be multiples of 2 for convenience.Figure 1b illustrates () for  ≪  .For the typical case of  1  ≪  0 and introducing  0 = 2 1  0 , we get () ≈  cos(2 0  0 +  0  + 2 0  cos(  )).In the frequency domain, it consists of an amplitude  0 ( 0 ) at the carrier frequency  0 and two first harmonic sidelobes at  0 ±   with an amplitude of  1 ( 0 ) ≈  0 .This is depicted in Fig. 1c.In OCT, the frequency-domain analysis is performed by discrete Fourier transform () = {()} for each A-line data.The frequency resolution is about equal to the half of the A-line rate   = 1/.So, the first harmonic sidelobes overlap with the main peak when   < 0.5  .This condition has always been met in conventional OCE.The carrier and modulated components interfere with each other in the pixel at  0 with a phase  ≈ 2 0  0 + 2 0  cos(  ).From the amplitude of the phase variation,  is readily measured.This method has been commonly used in OCE to date.
We realized that the condition   <   /2 is not necessarily required.When   >   /2 , the modulation sidelobes separate from the carrier (Fig. 1d).In this case, one cannot directly apply the conventional method.We rewrite ()  0 ≪1 → { 2 0  0 [  0  +  0  ( 0 −  ) +  0  ( 0 +  ) ]} .Now that the modulation component is separated from the carrier, the vibration amplitude  0  manifest itself in the amplitude of the Fourier component at  0 ±   .In this case,  can be measured from the amplitude, not the phase.However, this is only possible for the mirror-like sample when there is no other signal (background) present at  0 ±   .Consider a tissue-like sample that has continuous scattering points along its depth with ().The problem of retrieving  from the scatterer at  0 becomes more complicated.Now the Fourier component ( 0 −   ) contains not only the vibration signal (negative sideband) originated from the scatterer at  0 but also the time-independent reflection signal (carrier signal) from a different scatterer at  0 −   , where   =   /(2 1 ) .Depending on whether ( 0 −   ) > ( 0 ) or not, the phasor circle may or may not exclude the origin, and  is embedded differently in the phase and amplitude of ( 0 −   ) .Furthermore, ( 0 −   ) also contains a vibration signal (positive sideband) originated from another scatterer located at  0 − 2  .Figure 1e (i) shows a phasor diagram illustrating this situation.
Our solution to this problem is as follows.In one embodiment, we may add a large constant bias to ( 0 −   ) to move the rotating phase circles far away from the origin and then retrieve ( 0 ) from the time-varying amplitude of phase after proper renormalization with the bias and ( 0 ) = 〈( 0 )〉, where 〈 〉 denotes time average.In another embodiment, the background 〈( 0 −   )〉 is subtracted from ( 0 −   ) .This brings the phasor circles to the origin as depicted in Fig. 1e (ii).( 0 ) and ( 0 − 2  ) can then be identified from a Fourier transform of  − 〈〉.We used this second algorithm in experiments.Its rigorous mathematical description is provided in Methods.The situation of   >   /2 is equivalent to signal aliasing that occurs when the signal frequency is greater than the sampling rate.In fact, our algorithm works for any   values regardless of whether it is aliased or not.Thus, it is a general algorithm for OCE.

Experimental demonstration of anti-aliasing demodulation
We used a swept-source OCT system previously built using a polygon-scanner wavelength-swept laser providing a center wavelength of 1307 nm, 3-dB bandwidth of 80 nm, axial resolution of 16 μm, A-line rate   of 43.2 kHz, and maximum average optical power of 12 mW on samples (Supplementary Fig. S1).To test our demodulation scheme, we used a piezoelectric transducer (PZT) actuator block as a mirror-like sample and apply a sinusoidal waveform synchronously with a master clock in the OCT system (Supplementary Fig. S2).To generate mirror-like reflection, a flat glass plate was attached to the PZT.The laser power on the sample was attenuated to generate a signal-to-noise (SNR) of 40 dB at the air-glass interface.Fig. 2a shows a typical A-line profile in the Fourier domain when   was 679.6 kHz.The signal appeared at the 107 th pixel.The interval between pixels is 21.5 µm in depth or 52.7 kHz.Two first harmonic sidelobes emerged at about 12 pixels away from the main peak.The side lobes were lower by 20 dB in signal power (10 dB in amplitude) than the main lobe.Since the power ratio is 20  10 ( 0 ), we get  =20.4 nm.The sideband position increased linearly with the vibration frequency, with the expected slope of 0.019 pixel / kHz (Fig. 2b).

Noises and time jitter correction for ultrasonic frequencies
Various noises ultimately limits the system's sensitivity to vibration and the accuracy of wave velocity measurement.The fundamental limit comes from optical SNR, , in A-lines.The length of the modulationinduced phasor is ( 0 ) 0 ( 0 ) .The minimum detectable amplitude   is then given by ( 0 ) 0   ( 0 ) =   ( 0 −   ) , where   ( 0 −   ) =   ( 0 ) = ( 0 )( 0 ) −0.5 is the minimum detectable reflection coefficient by definition 35 .So,  0   =  −0.5 , same for non-aliased and aliased regimes.We note that although ( 0 ) is obtained from sidebands at  0 ±   , its error is governed by the SNR of the main peak at  0 .In OCE,  is typically extracted from a set of M-mode data consisting of  Alines.Then, the sensitivity is improved by √ , and  0   =  −0.5  −0.5 .For the above experiment, where  = 108 and  = 10 4 , we find   =195 pm.Another fundamental noise source is 1/f electrical noise.SNR at low frequencies may be limited by the 1/f noise rather than the reflectivity-limited phase noise.Mechanical vibrations of optical components add noises.In our system, the mechanical jitters of beam scanners and rotating polygon filters were suppressed to practically acceptable levels.
Additionally, time jitters in electronics boards and signal generators generate phase noise.This noise increases in proportion to the modulation frequency and can become dominant at the ultrasonic range.
In our system, we found that the sweep period  has a Gaussian jitter with a standard deviation   of 4.9 ns (Supplementary Fig. S3).The time error accumulates through successive sweeps (Fig. 2c), resulting in a phase error   .Let   denote the period jitter in the -th sweep, where 〈  〉 = 0 and 〈  − 〈  〉〉 =   .The timing of the -th A-line is given by   =  + ∑    =1 .The output from the -th A-line is   =  0      . is determined from the inverse discrete Fourier transform of   at   .The result is  0 (1 +   ∑ ∑    =1  =1 /) .The standard variation of the second term is   ≈     √( + 1)(2 + 1)/6 2 (Supplementary Note 1).For  ≫ 1 , we get   ≈ √/3     .The total phase noise Δ is a square sum of the SNR-and jitter-induced noises: Δ = √1/ 2 + 0.58(    ) 2 .To measure the time jitter in real time and correct for it, we recorded the modulation waveform by using another digitizer channel in the system.The waveform was generated by a trigger synchronized with a particular laser wavelength, so it is supposed to contain the same time jitters as in the sweep period.We obtained the phase of the modulation waveform (applied to the PZT) and subtracted it from the phase of the demodulated signal of each M-scan (Supplementary Note 1).This phase correction was effective (Fig. 2d).The total phase noise after the correction coincides with the SNR limited noise (Fig. 2e).
Finally, the phase noise causes uncertainty in determining the wave velocity,  , from the phase variation along the propagation distance.The error Δ in velocity is given by Δ/ = Δ/2 (/) , where  is the effective measurement length.In soft materials, elastic waves are attenuated, and  is typically a couple of wavelengths.Therefore, Δ = 0.01 yields a high accuracy of Δ/ < 1%.

Ultrasonic OCE of hard materials
were excited on the surface by using a custom-made contact probe with the piezoelectric actuator (see Methods).For a 5 mm-thick acrylate plastic block, we obtained a flat dispersion curve over 1 to 2 MHz (Fig. 3a).Since the wavelength (< 1.4 mm) is smaller than the thickness, the excited wave is the Rayleigh surface wave with a velocity   = (0.862 + 1.14) (1 + ) ⁄   , where   is the pure (bulk) shear wave speed and  denotes the Poisson's ratio.The shear elastic modulus  is given by  =   2 , where  is the density.From the measured   = 1285  16 m/s and known = 1.2 g/cm 3 and = 0.37, we obtained = 2.26  0.06 GPa.This value at MHz is slightly higher than the literature value of 1.7 GPa measured by quasistatic mechanical tools 36 .The difference is presumably due to the weak viscoelasticity of the material.Next, we measured a 1 mm-thick polystyrene plate.Figure 3b shows a representative -domain plot measured at 787.6 kHz.The dispersion relations for the A0 and S0 waves were fitted using the Lamb wave model 37 .
The result was = 1.64  0.05 GPa, in comparison to quasi-static shear modulus of ~ 1.3 GPa reported in literature 38 .We also performed OCE on glass and metal (Supplementary Fig. S4).We obtained = 24.0  0.61 GPa for borosilicate glass and 45.7  1.79 GPa for copper (Fig. 3c).These values were consistent with their quasi-static shear moduli 39,40 , indicating their near-pure elasticity.We also performed OCE on the cortical surface of a bovine tibia bone specimen (= 1.6 g/cm 3 and = 0.3) and found  to vary from 2.36  0.04 GPa at 485 kHz to 3.60  0.12 GPa at 2 MHz (Fig. 3d).This frequency dependence is in part due to the viscoelasticity of the bone but also contributed by its depth-dependent variation of modulus.Later we will show how the wave velocity dispersion can be used to measure depth-dependent elasticity.

Dynamic shear analysis of uniform soft materials
Next, we measured soft polymer samples using surface waves excitation in a frequency range of 0.1 kHz to 100 kHz. Figure 4 displays the experimental and analysis data from a rubber (Ecoflex 00-50), polydimethylsiloxane (PDMS), and a tough hydrogel 41 , respectively.Figures 4 (a, d, g) show measured surface displacement profiles at 40 kHz.The profiles contain not only the Rayleigh waves but also supershear surface waves 42,43 , which is evidenced by the two separated peaks in the wavenumber domain (Supplementary Fig. S5).We isolated the Rayleigh wave and derived the complex wave number   +   of the bulk shear wave (Supplementary Fig. S6).The measured speed ( =     ⁄ ) and the coefficient of attenuation (  ) are shown in Fig. 4 (b, e, h).The complex shear modulus,  ′ +  ′′ , is related to the complex wave numbers via: The measured storage and loss shear moduli are plotted in Fig However, the dramatic increase of ′′ beyond 10 kHz is only revealed by the ultrasonic OCE.For PDMS, our data up to 130 kHz show no crossover between ′ and ′′.The viscoelastic modulus of PDMS in the high frequency range was previously estimated only by the time-temperature superposition method using DMA data below 100 Hz 45 .Among the three materials, the tough hydrogel sample showed the lowest viscosity.The elasticity of this type of hydrogel has recently been measured using DMA below 100 Hz 41 .The low viscosity is due to the water that lubricates highly-entangled polymer chains and reduces energy dissipation 46 .Our data reveals the unusually low ′′ over the entire acoustic range, suggesting that this material is good at transmitting sounds and even ultrasounds.

Depth-resolved stiffness mapping of cartilage
We applied OCE to bovine cartilage in the knee joint ex vivo (Fig. 5a).The firm, connective tissue consists of three layers with different mechanical properties: the noncalcified cartilage (NCC), calcified cartilage (CC) and the bone (Fig. 5b).In our three-layer model, the thicknesses and shear moduli of NCC, CC and bone layer are denoted by hi and μi (i = 1, 2, 3), respectively.The NCC layer is readily distinguishable in the OCT image, from which we determined h1 = 300 μm.Unfortunately, the lower part of the CC layer and the underlying bone are not visible due to the limited optical penetration depth.We estimated the thickness of the CC layer to be h2 = 2.5 mm from cross-sectional cuts of equivalent samples (see Methods) and assumed h3 ≫ h2.We measured the elastic wave motions over 10 to 500 kHz. Figure 5d shows the dispersion of the Rayleigh waves, which shows several features.Besides a dramatic transition between 20 and 40 kHz, there appears to be a weak transition between 80 and 100 kHz.The Rayleigh surface wave has a 1/e amplitude decay depth equal to, approximately, a half wavelength.Therefore, as the frequency increases, the wave is increasingly confined near the surface, and its speed reflects the average stiffness of the region.The high speed below 20 kHz is due to the hard bone (μ3 = 3.27 GPa, see Fig. 2d).The first downward transition from 20 to 40 kHz is related to the bone-CC interface as the Rayleigh wave moves away from the bone.And the downward transition from 80 to 100 kHz is attributed to the CC-NCC interface and indicates that μ1 < μ2.We developed a three-layer finite element analysis (FEA) model (see Method) and plotted the best-fit result in Fig. 5c, from which we determined μ1 = 5.6 ± 0.2 MPa and μ2 = 13.2 ± 0.9 MPa.
Our demodulation method allowed us to generate subsurface wave maps (Supplementary Fig. S7).
Figure 5d shows cross-sectional displacement images obtained at 100, 140, and 500 kHz, respectively.It is apparent that the wavelength is longer in NCC than CC.We performed FEA using the geometrical and elastic parameters obtained from the curve fitting to generate displacement maps.The simulation reproduces the measured displacement fields quite well.

Depth-resolved stiffness mapping of in vivo human skin
We next examined the stiffness of human fingertip skin in vivo over 100 Hz to 500 kHz.We have previously measured the moduli of the epidermis, dermis (D), and hypodermis (H) in the forearm skin by OCE at acoustic frequencies 33 .With ultrasonic frequencies, we were interested in resolving the stratum corneum (SC) 47 and the viable epidermis (VE) layers within the epidermis (Fig. 6a). Figure 6b show cross-sectional displacement images at different frequencies.At 20 kHz, the Rayleigh wave is confined predominantly in the SC and VE layers.Above 300 kHz, the majority of the elastic energy resides in the SC layer.
The wave velocity dispersion curve has a strong frequency dependence in the normal skin condition (Fig. 6c, red circles).To relate the dispersion to shear modulus, let us consider a depth-profile of elastic shear modulus, ().There are no analytic solutions for wave propagation for arbitrary nonuniform modulus.For small nonuniformity, the wavenumber can be estimated from , where () is the wave function in a uniform medium with a space-averaged modulus: |()| 2 ≈ 2   −2   , where   ≈ 0.31  .We may extend the perturbation theory to large variation with a simplified form: where   denotes the penetration depth of the Rayleigh wave.We may choose   = 1/2  ≈ 0.25  or more generally,   =  = /, where  is a fitting constant.By differentiating both sides of Eq. ( 3) with respect to , we obtain where  = / ,  = 1.048   is bulk shear velocity obtained from the measured Rayleigh wave velocity   (), and  ′ = /.Figure 6d shows () obtained from the experimental data.We find () to decay rapidly from 56 MPa at the surface to 100 kPa at the asymptotic depth (using = 0.25).
This finding is consistent with the fact that water content decreases exponentially from the VE to the surface of SC 48 and that the elastic modulus decreases exponentially with depth in SC.
We measured the same skin immediately after the fingertip was immersed in warm water for 20 min.
The thickness of the SC layer increased from 0.31 mm to 0.41 mm due to swelling (Fig. 6e).The wave velocity at high frequencies decreased significantly (Fig. 6c, yellow dots), indicating a softening of the SC layer.The Rayleigh waves became overdamped at frequencies above 50 kHz due to high viscous damping.Using Eq. ( 4), we obtained = 2.6 MPa at the surface of the hydrated SC.
To verify, we performed FEA simulation by modeling the skin as a four-layer structure.We measured the shear moduli of the dermis ( 2 = 9-17 kPa) and hypodermis ( 3 = 2-3 kPa) layers from wave dispersion below 1 kHz 33 .The shear modulus in the SC layer at the pre-hydrated condition was assumed to have an exponentially decreasing function from  0 = 50 MPa at the surface to  1 = 2 MPa at the interface between the SC and VE layers (Supplementary Fig. S8).In the hydrated condition, the stiffness of the swollen SC layer was assumed to be uniform and the same as the VE layer ( 0 =  1 = 2 MPa).In Fig. 6c we plot the dispersion relations obtained by the FEA simulation, which captures the salient features of the dispersion relations.

DISCUSSION
We have demonstrated an OCE system capable of visualizing elastic waves within various materials and tissues over a wide range of shear modulus from 1 kPa to 10 GPa and across an unprecedented frequency range from static to a few MHz.The fundamental reflectivity-limited performance of the system provided sufficient SNR to determine wave velocities with high accuracy.The measured speed dispersion over the wide frequency range was critical to extract depth-dependent stiffness from multi-layered tissues commonly found in various tissues.The contact area of the PZT probe was optimized to be less than half the wavelength to maximize acoustic impedance matching from the actuator motion to the elastic waves.
Throughout the experiments, the peak wave amplitude ranged from several hundreds of nm at acoustic frequencies to as small as 10 nm at ultrasonic frequencies.
The demodulation scheme was highly effective and allowed the upper frequency limit to be no longer limited by the A-line rate.This relaxes the system requirement.The demodulation scheme should be applicable to any swept source OCT systems regardless of their A-line rates.It should be noted that the aliasing scheme does not work for spectra-domain OCT because signal modulation during the A-line acquisition period is averaged out due to the integrating nature of CCD spectrometers (equivalent detector bandwidth is only   /2).This is analogous to the phase washout motion artifact in spectral-domain OCT 34 .
Swept-source OCT is free from the phase washout as long as   is lower than the detector bandwidth, which is typically 1000 times   .
We envision ultra-wideband OCE to be a powerful tool for mechanical characterizations of soft and hard tissues with high spatio-temporal resolution.For example, ultra-wideband OCE may be useful in capturing the full spectrum of the complex rheology of tissues [49][50][51] and biomaterials 52,53 .With the capability of 3D high-resolution stiffness mapping, ultra-wideband OCE has the potential to enable novel biomechanics-based clinical diagnostics [54][55][56] .Besides, the ability to detect up to MHz excitation may be useful to study fast neuronal activities through neuro-mechanical coupling and develop novel therapeutics to treat brain diseases 57 .Furthermore, ultra-wideband OCE may be used to extract mechanical stress in thin structures without prior knowledge of material properties 58 .Finally, the emergence of deep learning offers a new approach to apply the information-rich elastography data to multi-parameter inverse problems 59 .

MATERIALS AND METHODS
Swept-source OCT signal.The photodetector current () of the OCT interferometry from a point scatter located at  0 can be expressed as () = () cos(2 0  0 + 2 1  0 ), (5)   where  denotes the reflectance amplitude, () denotes the optical power profile,  =  0 +  1  denotes the wavenumber of the swept source.The Fourier transformation of the photocurrent acquired during time interval −  2 ⁄ ≤  ≤  2 ⁄ with respect to  ̂= 2 1 , denoted by (), is given by where  is the period of wavelength sweep (A line period).We assume that the optical power profile has a Gaussian shape with a full width at half maximum (FWHM) , i.e., ( where 'c.c.' denotes the complex conjugate.When  < 1 , we can approximate the integral range to infinity and get where we introduce  0 () that is defined as where   ≡ 4ln2  1  corresponds to the FWHM axial resolution.
We define The Fourier transform of the photocurrent ( ) is where '*' denotes the convolution.
( 0 ;  0 − when  is an odd number and we get the apparent frequency ( ̅  − 0.5  ).
Note that the sum ∑   ( 0 ±   ;  0 )  =1 = 0, given that     ⁄ is an integer, where  is the number of A lines.We have =   ( 0 ;  0 ), the static OCT signal of the scatter located at  0 .Subtracting this term from ℱ  ( 0 ) we can get, The left side of Eq. ( 21) contains the left sidelobe of the scatter S' and the right sidelobe of scatter S''.Here we show the extraction of the vibration for S' as an example.Note that =   ( 0 +   ;  0 +   ), the static OCT signal of S'.
( 0   0 +  ) (  +) +   ( 0 −  ; 0 ) where   0 +  denotes the vibration amplitude at  0 +   .The first term of the left side in Eq. ( 22) gives algorithm to extract the surface displacement profiles.Next, we performed a 1-dimensional Fourier transform to move the data from time  domain to frequency  domain.The frequency domain data was filtered at the driven frequency (when   <   /2) or the aliased frequency (when   >   /2) to obtain lower noise waveforms.After we obtained the displacement profiles over the  coordinate, we performed another 1-dimensional Fourier transform to move the data from the spatial  domain to the wavenumber  domain.The wavenumber  of the surface wave was then determined from the plot by selecting the peak corresponding to the Rayleigh surface wave.This filtering in the - domain is critical to remove other higher-order modes especially at high frequencies.The phase velocity is then given by  = 2  ⁄ .To obtain the cross-sectional displacement image, the same method used for obtaining the surface wave displacement was applied throughout the entire pixels along the depth .In addition, the motion artefact caused by the surface wave motion were corrected within the sample 63 .
Measure complex wave number from surface waves.For soft materials, the near-field displacement is primarily dominated by the Rayleigh and supershear surface waves in high frequency regime.An analytical approximation solution for the surface waves excited by a cylindrical probe can be obtained.The surface displacement at lateral coordinate  is 62 where  1 is the Bessel function of the first kind of order 1,  0 (1) is the Hankel function of the first kind,  is the radius of the cylindrical probe,  0 is the pressure amplitude exerted to the sample by harmonic vibrations of the probe.  is the vibration frequency.,   and   are wavenumbers of the shear, Rayleigh surface, and supershear surface waves, respectively.  and   are roots of the secular equation ℱ() = (2 2 −  2 ) 2 − 4 2 √ 2 ( 2 −  2 ) sign{Re( 2 −  2 )} = 0 .So, we have   = 1.047 and   = (0.4696 − 0.1355).ℱ ′ () = dℱ d ⁄ .To derive  from the experiments, we fit the real and imaginary displacements simultaneously with Eq. ( 23) using the least-squares method.
Representative fitting results can be found in Supplementary Figure S7.
Finite element analysis.Finite element analysis (FEA) was performed with Abaqus 6.12 (Dassault Systèmes).We used plane strain models and performed time domain simulations.A localized harmonic surface pressure was applied to excite elastic waves, simulating the PZT actuator in the experiments.Other sides of the model were completely constrained.The geometry and total time of the model was scaled by the wavelength and stimulus frequency   , respectively.The time increment and total time were 0.1   ⁄ and 20   ⁄ , respectively.The width and height of the model were about 16 folds of the maximum wavelength to avoid reflection from the boundaries.The model was divided into uniform layers, each with different thickness and elastic parameters.To model a layer with a gradient property, we assigned a temperature-dependent elastic modulus to the layer and then prescribed the temperature field with an analytical distribution.The gradient mesh was adopted to reduce the computation costs.At the surface the mesh size was smaller than one tenth of the wavelength and at the bottom the mesh size was about half of the wavelength.The element type used in this study was 8-node biquadratic element (CPE8RH).
The convergence of the model was checked by making sure the results independent of the mesh size.
Preparation of hard materials.The acrylate plastic block (McMaster-Carr) presented in Fig. 3a had a thickness of 5 mm and an area of 50 mm×50 mm.The polystyrene petri dish (Fisher Scientific) presented in Fig. 3b had a bottom thickness of 1 mm and a diameter of 150 mm.The borosilicate glass coverslip (Fisher Scientific) measured in Supplementary Fig. S5a had a thickness of 0.15 mm and an area of 24 mm×50 mm.The copper foil (All Foils Inc.) measured in Supplementary Fig. S5b had a thickness of 0.3 mm and an area of 15 mm×40 mm.The glass coverslip and the copper foil were placed on a lens holder so that the sample was bounded by the air on two sides.
Preparation of soft materials.The bulk silicone rubber has a diameter and height of 60 mm and 12 mm, respectively.It was prepared from Ecoflex 0050 material (Smooth-On Inc) by mixing the Ecoflex 1A and 1B at 1:1 ratio by weight.The mixture was poured into a mold and cured at room temperature overnight.The material was then post-cured in an oven at 80ºC for 2 hours.The bulk PDMS sample has a diameter and height of 120 mm and 50 mm, respectively.It was prepared by using a 2:1 mixing ratio of base elastomer and curing agent (Sylgard 184, Dow Corning) and cured for 45 mins at 85 C.The hydrogel sheet had a size of 40 mm×40 mm and a thickness of 1.3 mm under fully swollen, highly entangled condition 41 .All soft materials were assumed to have a mass density of  ≃ 1000 kg/m 3 .
Cartilage tissues.Two fresh bovine tibia bones from two juvenile calves were obtained 1 h post-mortem (Research 87 Inc., Boylston, MA).The tissues were wrapped in a wet towel to keep them well-hydrated until use.Three measurements were performed on the cortical surface of the bone.Four measurements were performed on the femoral condyle region covered by articular cartilage.To measure the thickness of calcified cartilage layer, we made a cross-sectional cut after the OCE experiments.
Human subject.We measured the skin of the left index finger on a 31-year-old male subject.The study was conducted at the Massachusetts General Hospital (MGH) following approval from the Institutional Review Board (IRB) of Massachusetts General Hospital and the Mass General Brigham Human Research Office.Written informed consent was obtained from both subjects prior to the measurement and all measurements were performed in accordance with the principles of the Declaration of Helsinki.The measurement site was marked by a surgical grade skin marker.For the hydration test, the whole index finger was immersed in warm water for 20 min, wiped with a clean towel, and then measured immediately.Depending on whether ( 0 −   ) is greater or smaller than the magnitude of rotating phasor ( 0 ) 0 ( 0 ), the pixel phase undergoes an oscillatory or diverging pattern, respectively.For demodulation, the rotating phase circles (i) is brought to the origin (ii) by subtracting the offset ( 0 −   ), and ( 0 ) is determined from the magnitude of the counterclockwise rotating phasor.

Fig. 2 Fig. 3 Fig. 4 :
Fig. 2 Implementation of aliasing in OCE.a, Signal from the surface of a stationary PZT in the frequency (depth) domain without modulation (top) and with modulation at 679.6 kHz (bottom).b, Measured position of modulation peak at different modulation frequency.c, Several M-scan traces at 679.6 kHz before and after time jitter correction.d, Phase data of the M-scan traces showing noises before and after jitter correction.e, Phase noise performance measured with optical SNR of 40 dB and = 108.The SNR-limited phase noise corresponds to -94 dB/Hz.The 1/f noise is also seen.Note that this graph is for 40 dB SNR.For lower SNR, the optical phase noise level would increase.For example, for 20 dB optical SNR, the optical noise level is moved up to 0.07 rad, crossing the 1/f noise line at 200 Hz.

Fig. 5 Fig. 6
Fig. 5 Depth-resolved stiffness mapping of joint tissue ex vivo.a, Picture of a bovine joint sample and the PZT contact probe.b, Cross-sectional OCT image with three layers labeled: noncalcified cartilage (NCC), calcified cartilage (CC), and the bone.The shear moduli and thickness of the NCC, CC and bone layers are denoted by μi and hi (i = 1, 2, 3), respectively.c, Measured Rayleigh surface wave speeds (circles) and the best-fit FEA result using the three-layer model (dashed curve).Error bars represent the standard deviation of four measurements from two samples at two locations.d-f, Side by side comparison of the measured cross-sectional displacement image and the FEA simulation at 100 kHz, 140 kHz, and 500 kHz, respectively.Dashed lines indicate the approximate boundaries of the NCC layer.

Fig. S2 .
Fig. S2.Timing diagrams of the OCE system.The FBG optical clock provides the sync signal for the I/O board that controls the stimulus signal, the OCT beam scan, and data acquisition.For stimulus frequencies below 21.6 kHz, the I/O board generates the harmonic stimulus directly.For stimulus frequencies above 21.6 kHz, the I/O board generates a stimulus trigger to externally trigger a function generator to send the harmonic stimulus.

Fig. S3 .
Fig. S3.Gaussian distribution of the time jitter in sweep period .Data represents 900 M-scans, with each M-scan consisting of 108 A-lines.By evaluating the phase difference between two consecutive A-lines, we obtained the time jitter over 900*107 = 96,300 data points.The distribution is fitted with a Gaussian curve (dashed line).The standard deviation of the time jitter is 4.9 ns.

Fig. S4 .
Fig. S4.Ultrasonic OCE of the glass and metal.a, Dispersion relation of the coverslip glass (dot) and theoretical fitting with a Lamb wave model (line).b, Dispersion relation of the copper sheet (dot) and the theoretical fitting with a Lamb wave model (line).

Fig. S5 .
Fig. S5.eavenumber k domain plots obtained by the Fourier Transform.Comparison is made at 40 kHz.The Rayleigh surface wave mode (R) and the supershear surface wave (SS) mode are labeled on the plots.

Fig. S7 .
Fig. S7.Cartilage OCE.Cross-sectional displacement images at 40 -140 kHz for a, eithout applying the demodulation algorithm for inner scatters, and b, After applying the demodulation algorithm for inner scatterers.

Fig. S8 .
Fig. S8.FEA model for the skin.The depth dependent profile of shear modulus in the human skin for the normal, pre-hydrated state (margenta) and hydrated state (turquoise).
44 323 kPa, respectively, over 0.2 kHz to 50 kHz.′′growsat a faster rate above 10 kHz and becomes greater than ′ above 40 kHz.The crossover of ′ and ′′ causes strong wave attenuation (     ⁄ = 0.41).When     ⁄ > 1 , waves become overdamped, and their velocities cannot be determined.At low frequencies, our data agreed with previous measurements by MRI elastography, which showed ′ and ′′ to increase from 30 to 80 kPa and from 10 to 40 kPa over 0.2 to 7.8 kHz44.
4 (c, f, i).They show characteristic frequency dependencies due to viscoelasticity.For the rubber, ′ and ′′ increase from 50 to 160 kPa and from 9