A Conformable Ultrasound Patch for Cavitation‐Enhanced Transdermal Cosmeceutical Delivery

Increased consumer interest in healthy‐looking skin demands a safe and effective method to increase transdermal absorption of innovative therapeutic cosmeceuticals. However, permeation of small‐molecule drugs is limited by the innate barrier function of the stratum corneum. Here, a conformable ultrasound patch (cUSP) that enhances transdermal transport of niacinamide by inducing intermediate‐frequency sonophoresis in the fluid coupling medium between the patch and the skin is reported. The cUSP consists of piezoelectric transducers embedded in a soft elastomer to create localized cavitation pockets (0.8 cm2, 1 mm deep) over larger areas of conformal contact (20 cm2). Multiphysics simulation models, acoustic spectrum analysis, and high‐speed videography are used to characterize transducer deflection, acoustic pressure fields, and resulting cavitation bubble dynamics in the coupling medium. The final system demonstrates a 26.2‐fold enhancement in niacinamide transport in a porcine model in vitro with a 10 min ultrasound application, demonstrating the suitability of the device for short‐exposure, large‐area application of sonophoresis for patients and consumers suffering from skin conditions and premature skin aging.


Introduction
Transdermal drug delivery offers an attractive alternative to conventional drug delivery methods of intravenous [1,2] and oral frequency, and each frequency regime enhances skin permeability through different mechanisms. [14,32] Traditionally, HFS has been the modality of choice and has been shown to improve the penetration of small molecule drugs and their efficacy in vitro and in vivo. [33,34] HFS can be used to nucleate cavitation within the stratum corneum, where the bubble size (radius < 20 μm) is comparable to the intercellular distance of the lipid bilayers. However, it was later found that LFS, which generates larger bubbles  μm in radius), could be more effective than HFS [35,36] due to stronger transient phenomena such as shockwave generation, bubble collapse, and microjetting. [24,32] While the cavitation threshold is smaller with LFS, HFS generates higher bubble count densities and velocities at the expense of electrical power. Further, there is a diminishing return on bubble count after a certain frequency. [37] Thus, IFS offers an interesting compromise between the convective effects and intensity of cavitation.
HFS treatments can be applied with the probe in direct contact with the skin, but the enhancement effects with lower frequencies cannot be obtained without a fluid coupling medium in between. Currently, this involves completely submerging the target in the fluid or attaching flanged cylindrical fixtures to create "cavitation chambers" on the skin surface, [38,39] which is not a feasible or scalable solution for non-medical applications. Further, manual positioning of the probe introduces high variability in treatment efficacy. This is not surprising given that acoustic cavitation is a highly sensitive stochastic phenomenon, where multiple systems as well as environmental parameters such as transducer properties, spatial geometry, [40] and choice of coupling medium, [41] dissolved gas content, [42] and heterogeneity and separation of target membrane profoundly impact each instance of cavitation. There is a need for a platform that can ad-dress the above issues and support the high-enhancement mechanisms offered by low-frequency sonophoresis while maintaining conformal and flush contact with the skin. A summary of the existing solutions and our device is presented below in Table 1.
In this study, we report a hands-free IFS system, the conformable ultrasound patch (cUSP), with small (0.8 cm 2 ) piezoelectric discs (PZT-Ds) optimally positioned within a conformable poly(dimethylsiloxane) (PDMS) substrate. [49][50][51] Leveraging a conformal platform allows tight control of the relative spacing between the ultrasound source and the curvy target membrane. The resulting cavity provides a 1 mm-deep reservoir space for the coupling fluid medium wherein inertial cavitation, convective mixing, and microjet formation can be induced on the surface of the skin. The cUSP system removes the need for large, heavy hand-held ultrasonic probes and cavitation chambers and achieves delivery of IFS through tuned PZT-Ds and self-contained small-volume (400 μL) drug reservoirs integrated within the patch (Table 1).
First, the device design and system parameters such as the working frequency, applied voltage range, and acoustic pressure distribution are determined using a fully coupled electromechanical finite element (FEM) simulation software (COMSOL). The modeling is performed at multiple levels with results experimentally verified at each stage, starting with bare actuator dynamics in air and evolving to include fluid-structure interactions of the piezoelectric transducers within the substrate. Together, the simulation model, acoustic characterization, and high-speed imaging of the cavitation bubble dynamics provide a holistic characterization of the cavitation regime that are used to guide parameters for device design such as the separation distance of the piezoelectric transducer from skin boundary. Niacinamide (NIA) is chosen as the candidate cosmeceutical drug in conjunction with the Figure 1. Device concept. a) Photograph of the 2D-array of conformable ultrasound patch (cUSP) (diameter 5 cm, thickness 2 mm) attached to a curvilinear glass cylinder;scale bar: 1 cm. b) Exploded view of the 2D-array of cUSP (5 cm × 5 cm × 2 mm) showing each constituent layer; scale bar: 5 mm (the inset shows the fluid cavity formed between the PZT-D element and skin; scale bar: 2 mm). c) Schematic illustration of the cUSP on skin, showing the cavitation mechanism within the cavity between the device and the skin, and resulting drug penetration through stratum corneum. Proportions of stratum corneum thickness to device size are distorted for the sake of clear illustration of the working mechanism. d) Acoustic pressure distribution in the device operating in radial mode at a frequency of 212 kHz and applied voltage of 50 V, showing an undamped pressure zone (≈100 kPa) within the medium in the cavity, which is enough to nucleate cavitation. cUSP, and the permeation enhancement provided is evaluated in vitro using a porcine model in a Franz diffusion cell. With further optimization and miniaturization of the cUSP and extensive in vivo human testing, we envision this novel approach to bring forth advancements in the pharmaceutical and cosmeceutical industry for better treatment and repair of the skin for patients and consumers suffering from skin conditions, damaged skin barriers, or premature skin aging.

Device Design
The conformable ultrasound patch (cUSP) consists of 4 piezoelectric elements (lead zirconate titanate (PZT)) embedded in a poly(dimethylsiloxane) (PDMS) substrate. First, each PZT disc (PZT-D), with a diameter of 10 mm and thickness of 2 mm is separately coated with an ultra-thin PDMS layer using a precisely controlled dip-coating method. This creates a thin-film, waterproof encapsulation which provides minimal damping of the PZT-D motion in the fluid medium to allow maximum energy transfer. No backing or matching layers are added to the PZT-D as the objective of the design is to create large undamped peak pressures at a single driving frequency very close to the PZT-D surface, rather than maximize the energy transfer over a large frequency bandwidth and decrease the Q-factor, as is desirable for imaging transducers. [52] Four PZT-Ds are then arranged in a 2×2 array configuration with a center-to-center separation of 20 mm along with the serpentine metal (copper, Cu) electrodes that establish their electrical connections. The system is placed in a custom 3D printed poly(lactic acid) mold and fully encapsulated in a 2 mm-thick PDMS layer. The PDMS encapsulation has square dimensions of 5 cm × 5 cm when it is initially demolded. For attachment onto the skin, the patch is trimmed to a circle with a diameter of 5 cm. The CAD file is provided in the Supporting Information Data. Figure 1a,b shows the optical photograph of the cUSP attached to a glass cylinder and an exploded view of the multilayer construction of the patch, respectively.
The distance of the ultrasound transducer from the target membrane shows an inverse relation with the count of bubbles. [53] Thus, the piezoelectric ultrasound transducer discs (PZT-D) are positioned in the substrate such that the disc surface is maintained at a small distance of 1 mm from the skin, while forming a cavity wherein the fluid medium can be contained (see inset in Figure 1b). An annular gap of 1 mm width is maintained between the sides of the PZT-D and the surrounding PDMS substrate to allow undamped oscillation of the PZT-D in the radial direction, which is the selected operating mode for the transducer as discussed in later sections. Ultrasound energy is applied to induce cavitation within this region as shown in the schematic illustration in Figure 1c. An FEM simulation map of the acoustic pressure field with this geometry reveals that the pressure amplitude in the cavity is sufficient to nucleate cavitation [54,55] and trigger violent bubble collapse [31] (Figure 1d). It must be noted that while a theoretically predicted value of ≈100 kPa is required to counter the atmospheric pressure (101.3 kPa) and trigger the degassing of dissolved air from an aqueous solution, in practice, the same effect may be achieved with lower pressures if impurities and gaseous nuclei are already present. [56,57] This has been demonstrated in the subsequent sections.

Electromechanical Characterization and Modeling of the Device
To understand the relation between the design parameters and the acoustic performance, we have developed finite element models (FEM) and experimentally verified them. The device was first modeled in air, where the mechanical boundary condition considered in the FEM is free motion on the edges of the device. The impedance spectrum from the FEM simulations is calculated based on the ratio between voltage and current for the electrodes over the desired range of frequencies. For our specific PZT material with N p of 2200 m s −1 , N t of 2070 m s −1 (see Experimental Section), diameter of 10 mm, and thickness of 2 mm, the theoretical resonance frequencies are calculated to be 220 kHz and 1.04 MHz. The impedance graph in the frequency domain measured by the FEM simulation is shown in Figure 2a. An impedance analyzer (E4990A, Keysight Technologies) was used to experimentally determine the impedance characteristics of the device in air before fluid-structure interactions. The measured impedance spectrum in the frequency domain is shown in Figure 2a along with that obtained by the FEM simulation. From the impedance magnitude spectra, clear resonant modes emerge for 5 frequencies. The lowest impedance values are obtained at 217 kHz and 1 MHz experimentally and at 225 kHz and 1.01 MHz from simulation. The slight shift (<4%) in the resonant frequencies between simulation and experiment likely results from the damping effect of the copper electrodes on the piezoelectric transducer surface.
For the entire range of frequencies examined, the impedance values predicted by FEM simulations agree remarkably well with the experimental data as shown in Figure 2a. The resonant frequencies measured experimentally are slightly smaller than the FEM predicted values, which can be attributed to higher damping in the real-world scenario due to differences in air density and the copper-tape soldered electrical contacts made on the electrodes. In water, the resonant frequency of the device is further lowered owing to the viscous damping effect of the fluid medium, as expected [58] (Figure S1a, Supporting Information). It should be noted that there is an inherent device-to-device variation (<1%) in the radial resonance frequency in the part as received from the supplier. For each time the disc is used, the exact resonance frequency is scoped in the neighborhood of 220 kHz to find the frequency corresponding to the lowest impedance value for the given disc, coating, and environment, and this frequency is then selected for further operation of the device to obtain optimal energy transfer for that environment.
The minimum impedance values are found to correspond to the radial (f r , 217 kHz) and thickness resonance modes (f t , 1 MHz). Accordingly, the displacement profiles are characterized and compared for these two modes going forward. The device is suspended in air and water under the 3D laser doppler vibrometer (LDV, MSA-500, Polytec) measurement system to measure the deflection across the top surface of the PZT-D as shown in Figure 2d,e. The frequency is swept from 0.2-1 MHz and the shape modes are reconstructed for the 2 resonant frequencies f r , f t as given in Figure S2, Supporting Information. The displacement plots from the simulated and experimental model in air are overlaid in Figure 2b for a frequency sweep of 0.2-1 MHz. The displacement under the effect of a continuous sinusoidal signal at the two resonance modes (f r , f t ) is considered as the final system will be driven in a pulsed or continuous mode. Figure 2c Figure S2c,d, Supporting Information. This can be attributed to the energy being focused in a narrow frequency band for multiple cycles for a continuous sine as compared to the application for a single cycle at each frequency step in the case of a sweep. The shape mode for the maximum deformation of the top PZT-D surface presents the same number of nodes (3 anti-nodes at f t , 1 anti-node at f r ( Figure S2c,d, Supporting Information)) as predicted by the FEM model ( Figure S2a,b, Supporting Information). A similar measurement is performed with the PZT-D submerged in water, and a negligible shift in resonance frequency is observed ( Figure S1a, Supporting Information).
The model is now extended to include the effect of the PDMS substrate on the vibrational characteristics of the device. Figure 2d,e shows the displacement maps at 220 kHz for a PZT-D with a 1 mm-thick PDMS backing. This configuration is chosen as the PZT-D element is embedded on a 1 mm-thick PDMS substrate forming the bottom encapsulation of the patch, while the top surface is free within the patterned cavity. Figure 2d depicts the displacement on the top surface of the PZT-D. The distribution is similar to that observed in Figure S2c, Supporting Information. The displacement scale is capped at 20 nm in the plot to allow an appreciation of the displacement distribution in the PZT-D, while discounting the disproportionately high displacements experienced in the PDMS. These are subsequently shown in Figure 2e, which depicts a log plot of the displacement on the bottom of the PDMS substrate. These values are 4 orders of magnitude larger than the PZT-D displacements, which follows from the ratio of the elastic moduli of the two materials (PDMS: 0.6 MPa, PZT: 60 GPa). However, there is no significant shift in the resonant frequencies of the PZT-Ds in the two cases (with and without PDMS backing substrate), and the displacements attained on the top surface of the PZT-D (the surface of interest coupled with the fluid medium and skin) remains undamped with the addition of PDMS (comparing with Figure S2a, Supporting Information), reaching a peak value at f r of 20 nm in both cases. It is to be noted that in the final device, there is an additional thin encapsulation layer for each PZT-D element before it is assembled onto the cUSP substrate. (see Experimental Section). However, the thickness of this layer is in the range of 70-85 μm ( Figure S13, Supporting Information) which is much smaller than the 1 mm-thick PDMS backing layer considered in Figure 2d,e and thus has a negligible impact on the PZT-D performance. Next, Figure 2f shows an overlay of the impedance spectrum of a single bare PZT-D, PZT-D with 1 mm-thick PDMS coating, and the cUSP (4 PZT-D elements in a 2D array, parallel configuration). An interesting observation can be made here. Since the 4 PZT-D elements are attached in parallel, we expect the impedance of the equivalent circuit to be 1/4 th that of a single element. Far away from the resonance peaks, it can be seen that the patch impedance is indeed shifted down to about 0.25× the single-element impedance values (for example, in the frequency range of 100-150 kHz in Figure 2f). However, near the resonant peaks, for example, at 220 kHz, the patch impedance (67 Ω) is greater than that of a single element (30 Ω). This behavior is caused due to a change in the relative contributions of the electrical and mechanical components of the system in the different operating modes. Off-resonance, the curve is dominated by the high impedance values of the PZT-D discs. At resonance, the impedance of the PZT-D discs is dramatically reduced, and the impedance reading is dominated by the mechanical damping offered by the PDMS substrate, which is larger in the case of the cUSP as compared to PDMS-coated single element. It is critical to take into account this electro-mechanical coupling effect between the PZT-D and the substrate in order to correctly characterize the impedance value of the system and optimize the power transfer characteristics. The FEM results match well with the experiment and show that the numerical method can accurately characterize the vibrational behavior of the device, and can thus be extended to simulate fluid interactions. We investigate the pressure regime for both modes (f r , f t ) in the following section to further consolidate the choice for the final operating mode of the device.

Acoustic Pressure Field Modeling and Characterization
The section above establishes the fidelity of the COMSOL model for accurately modeling the mechanical behavior of the PZT-D. Next, we proceed to understand the mechanical coupling and energy transfer between the PZT-D and the fluid medium. It is desirable to characterize the pressure gradient very close to the surface of the PZT-D to predict the location of onset of the cavitation regime. However, this is challenging to achieve experimentally with a hydrophone as the measurement becomes confounded with large electromagnetic coupling and strong reflections that interfere with the original pressure field when the hydrophone tip is brought very close to the PZT-D surface. Thus, a combined approach is employed wherein the experimentally measured pressure in the far-field is validated with the simulated values and the FEM model is then used to extrapolate the pressure field very close to the PZT-D surface.
Unlike 1D techniques of vibration analysis of piezoelectric discs which assume motion only in the thickness direction and . Acoustic pressure characterization of the device inside a water tank. a) Experimental measurement and numerical simulation results of acoustic pressure at f r at different distances away from the device along the z-axis at 10 V peak-to-peak. b) Experimental measurement and numerical simulation results of acoustic pressure under different applied voltages at f r 1 cm away from the surface of the device in the water. c) 2D contours of acoustic pressure distribution in a water column directly above the PZT-D surface for radial extension at f r (left) and for thickness extension at f t (right). Due to the symmetry, only half of the water column is considered. The applied voltage is 50 V peak-to-peak in each case. d) Comparison of the center-point pressure and average acoustic pressure in the water medium in the frequency range of 0.2-1.1 MHz. The centerpoint pressure represents the pressure directly above the center of the PZT-D (r = 0, z = 0) which is expected to be the location of maximum pressure for both resonance modes, and the average acoustic pressure is calculated across a 1 mm-thick slice above the entire surface of the PZT-D. e) Pressure distribution in the water cavity and skin shown in cross-section of the cUSP device taken diametrically through the center of 2 PZT-Ds.
ignore other modes of vibrations, the FEM analysis accounts for the 3D deformation of the PZT-D while calculating the acoustic pressure in the surrounding medium. As the geometry is rotationally symmetric, it is possible to set up the model as a 2D axially symmetric problem in the r-z plane where r is the radial and z is the thickness direction for vibrational patterns of radial stretching (f r ) and thickness compression (f t f t ) (refer to Experimental Section for further details on the FEM model).
To experimentally evaluate the acoustic pressure profile, the PZT-D is immersed into a tank filled with deionized (DI) water. A hydrophone (HGL-0200, Onda Corp.) pointed directly along the central axis of the PZT-D is used to measure the generated ultrasound pressure at different distances away from the device. The ultrasound pressure field is generated by applying a continuous sinusoidal signal to the PZT-D. The schematic of the experimental setup is shown in Figure S4, Supporting Information.

Figure 3a
shows the acoustic pressure measured at 1, 2, and 3 cm away when the PZT-D was driven at a constant frequency at the radial mode f r and sinusoidal voltage of 10 V peak-to-peak. A low voltage is selected in order to minimize the electromagnetic coupling noise in the system and allow accurate measurement of the pressure field. The pressure amplitude decreases from 7.6 to 2.0 kPa with the increase in distance from 1 to 3 cm. Later, we vary the applied frequency to the PZT-D while keeping the distance between device and hydrophone constant. The pressure profile at various distances for the frequency range 100-500 kHz is shown in Figure S3, Supporting Information. There is a sharp drop-off in pressure from 1 to 2 cm, suggesting that the pressure field is confined within a small radius around the PZT-D. The pressure amplitudes calculated from the FEM simulation show a good agreement with the measured values.
Next, we investigate the influence of applied voltage on acoustic pressure. We vary the applied voltage to the transducer and measure the acoustic pressure using a hydrophone at a constant frequency of f r and a constant distance of 1 cm between the device and the hydrophone (See Experimental Section and Figure S4, Supporting Information, for detailed experimental setup). As the applied voltage is changed from 45 to 100 V, the output acoustic pressure value increases from 51.3 to 116.4 kPa, maintaining a linear relationship (Figure 3b). Overall, there is excellent agreement between the experimental and simulated absolute pressure values for the far-field measurements at different driving voltages. This allows a confident extrapolation of pressure values at 50 V using the FEM model for the pressure field within 1 cm of the PZT-D surface. Figure 3c shows the 2D contours of the acoustic pressure distribution in water in the vicinity of the device at the two resonance modes. The results indicate the specific areas of the water that exhibit the highest acoustic pressure. The pattern and magnitude of the acoustic pressure distribution are distinctly different for the two main modes of vibration. For radial extension, the wave propagation in the water has a wide central lobe with a maximum pressure of ≈100 kPa occurring near the middle of the PZT-D. In thickness deflection, the maximum pressure of 450 kPa occurs within a narrow axial band directly above the center of the PZT-D and fluctuates significantly with an eventual drop-off towards the edges. Figure S5c,d, Supporting Information, shows the damping of pressure while moving outwards from the center of the PZT-D surface radially (r = 0) to the edge (r = 5 mm), and in the z-direction (z = 0) to the water surface (z = 10 mm).
Further, we evaluate the pressure field for the two resonant modes to determine the operating frequency for the system. Despite the occurrence of a higher peak pressure at 1 MHz, the spatial average of the pressure over a volumetric slice of 1 mm thickness over the top-surface of the PZT-D is more comparable for the radial mode and thickness mode as shown in Figure 3d. The cavitation threshold depends not only on the acoustic pressure but also on the irradiation time, and is thus typically lower for larger cycle times or lower frequencies. [31,32,56] It has been experimentally verified that the cavitation threshold at f t is almost double that at f r , [55] which is also observed by performing high-speed videography of the cavitation regime. Bubbles are nucleated at a threshold voltage of 50 V at the radial mode frequency (f r ). However, no bubbles are formed at this voltage when the device is operated at the thickness mode frequency (f t ). If the frequency is immediately switched from f r to f t , the streaming bubbles nucleated at the lower frequency come to a sudden standstill and oscillate in place (Video S5, Supporting Information). Therefore, the first radial mode of resonance (f r ) is selected for the operation of the device going forward to maximize both convective effects due to acoustic streaming as well as inertial cavitation events at a low threshold voltage.
The acoustic pressure studies presented above consider the absolute pressures generated by the device in a large water tank. However, the final use case for the device involves positioning an array of PZT-D elements in a restricted PDMS cavity, where the dimensions and material properties of the cavity walls and skin, as well as the spacing between the active PZT-D elements, can impact the acoustic pressure field. Figure 3e depicts an acoustic pressure map of the cross-section of the cUSP taken across the center of 2 PZT-D elements placed against the skin and driven in parallel at 50 V. The distribution represents the rms values of the pressure at a steady state for a 220 kHz driving signal. The pressure values are commensurate with the axis-symmetric model of a single element considered in Figure 1d, which suggests that there is little acoustic interaction between the PZT-D elements at this transducer spacing of 20 mm. However, interactions of the pressure fields may arise with smaller transducer spacings, and this effect is explored in Figure S13, Supporting Information. There is minimal distortion in the pressure field at the waterskin interface, which is expected as the acoustic impedances of the two materials are similar at 1.48 MRayl and ≈1.6 MRayl respectively. However, there is some degree of reflection of ultrasound energy both from the skin as well as the PDMS walls in the cavity, which creates a zone of constructive interference close to the top surface and sides of the PZT-D, and causes a spike in the peak pressure to 140-160 kPa, as compared to the peak pressure of 100 kPa observed for the case of a bare PZT-D in a water column as shown in Figure 3c. These simulations demonstrate that for the given 2D-array patch configuration, the pressure fields are well-contained above each PZT-D element, and the geometry of the PDMS cavity and application of the device against skin in fact serves to compound the peak acoustic pressures beyond the cavitation threshold.

Cavitation Detection and Bubble Dynamics
We proceed to study the bubble dynamics induced by the PZT-D system. In this experimental study, the onset and intensity of cavitation are first characterized by measuring the subharmonic response of the acoustic pressure field in the fluid medium. Through simultaneous collection of spectroscopic and acoustic data, Johnston demonstrated a clear correlation between bubble radius and shockwave period with the measured acoustic signal as a subharmonic of the driving frequency. [59] An increase in broadband noise and subharmonic signal strength are representative of inertial cavitation and the emission of pressure shockwaves. [31,60] A needle hydrophone (HNC-0400, Onda Corp.) placed close to the PZT-D is used to measure the subharmonic component of the acoustic field due to cavitation. The hydrophone is positioned perpendicular and out of the direct field of oscillation of the PZT-D to reduce the amplitude of direct vibrations from the PZT-D and allow pick-up of small signals of the cavitation field, which can be assumed to be spherically symmetric. Figure 4a shows a schematic representation of the cavitation measurement setup driven at the radial mode. The threshold voltage for cavitation is determined by gradually increasing the driving voltage until a subharmonic (f/2) or ultra-subharmonic (3f/2) signal occurs in the fast Fourier transform (FFT) spectrum of the hydrophone output signal that is significantly higher than the noise floor ( Figure S6a, Supporting Information). The ultrasubharmonic signal appears first, which is attributed to the fact that the sensitivity of the hydrophone is higher in this frequency range. The ultra-subharmonic signal and broadband noise are obtained by averaging 32 frames of FFT spectrum for 3 different devices. As seen in Figure 4b and Figure S6b, Supporting Information, the magnitude of the ultra-subharmonic peak increases almost linearly with increasing voltage. The standard deviation in the ultra-subharmonic signal magnitude starts increasing significantly when the applied voltage is 50 V and above, indicating sporadic appearance of ultra-subharmonic peaks. Given the stochastic nature of inertial cavitation, its dependence on existing bubble nuclei populations, [61] and short lifetime of inertial bubbles, [62] we see large variations in the measured ultrasubharmonic signal even over a long sampling duration (averaged over 32 frames, corresponding to a total sampling period of 1.6 ms (see Experimental Section)). However, since the total "ontime" for each ultrasound pulse for the final system (≈50 ms) is far greater than this sampling period, there is a much greater likelihood of obtaining inertial cavitation action for a significant portion of the application time. Thus, we refer to the acoustic pressure generated when the applied voltage is 50 V as the cavitation threshold, which corresponds to around 95 kPa. At an applied voltage of 90 V and above, the data becomes significantly different (p < 0.05) from that of 10 V, implying more frequent cavitation events occurring every frame. Similarly, the broadband noise magnitude increases almost linearly with increasing voltage (Figure 4c), and there is a sixfold increase in slope after crossing the cavitation threshold of 50 V (from 1.06 × 10 −8 to 6.42 × 10 −8 ). This suggests a higher count of bubbles and/or more violent bubble collapses at higher voltages. High-speed videography (Phantom v2512, Vision Research Inc.) is performed to provide further insight into the cavitation regime in the vicinity of the PZT-D. This reveals that indeed a larger number of bubbles at higher streaming speeds are generated as the voltage is increased, as seen in Figure 4c-e and Videos S3 and S4, Supporting Information. ImageJ is used to threshold the images and count the average number of bubbles across a stack of frames (see Experimental Section). The bubble count and streaming speed increase by 149% and 144% respectively as the voltage is ramped from 50 to 100 V. Correspondingly, a 15 dB increase in ultra-subharmonic signal strength is obtained for the same voltage gain. We expect bubble nucleation at the locations of maximum pressure above the PZT-D as predicted by the simulated pressure maps in Figure 3c. In radial mode, this corresponds to the center of the top surface and the side-walls of the PZT-D, and indeed a stream of bubbles can be seen originating from both these locations in Figure 4d and Video S4, Supporting Information. The relation between the applied voltage as www.advancedsciencenews.com www.advmat.de ultrasound pulse and the measured acoustic pressure maintains linearity up to 120 V, and no step change in the pressure field is observed. This suggests that the generated bubbles are continuously streaming and do not accumulate in a given region to create a pressure-dampening effect as seen with larger bubble clouds. [63] After identifying the bubble nucleation zone, 2 teleconverters are used to image the bubbles at a higher magnification (4×). The bubble radius is characterized from these images and is found to be in the range of ≈95-115 μm. A histogram of the distribution of the bubble radii is presented in Figure  S11, Supporting Information. The anticipated bubble size at f r is ≈15 μm from the well-known relation between frequency and bubble radius (f (kHz) × R (μm) = 3000). [32,64] However, after nucleation, the bubbles coalesce and grow within a few energy cycles and often the measured bubble size is significantly larger than the theoretically predicted bubble size. [65,66] This is a feasible assumption given that the maximum frame rate of the imaging setup is 10 000 fps, and the cycle length of the applied ultrasound is around 20 times smaller than this value.
The large standard deviation in the ultra-subharmonic data can be attributed to several factors in the experimental setup. The experiment is carried out with 3 different devices which may have minor differences in coating thickness and distribution across the body of the PZT-D, impacting the pressure profile in the water surrounding the PZT-D. Further, the total dissolved gas in the water changes with time, which influences the initial nucleation seeds available for the onset of cavitation. The FEM results presented in Figure 2d,e, indicate that it is less likely that the large standard deviations stem from varying coating thicknesses, as even a 1 mm-thick polymeric backing has negligible impact on the PZT-D displacements and the coating thickness used in the studies is significantly smaller than this value (<100 μm, Figure S13, Supporting Information). The dependency and variation with total dissolved gas should thus be carefully reduced by incorporating pre-existing nucleation seeds in the coupling medium for future embodiments of the system. Despite these discrepancies, the overall observed trend with the data across the 3 devices uniformly increases with increasing voltage. The phenomena observed above confirm a repeatable onset of inertial cavitation within the medium above a threshold voltage of 50 V. While the above characterizations have been performed on a single PZT-D element, the results are largely representative of the cavitation action expected in the 2D array patch. Figure 3e clearly shows that the high acoustic pressure fields are concentrated directly above each element, with minimal interaction between the PZT-D elements in the current configuration of the patch. For future embodiments with closer packing of PZT-Ds, it would be beneficial to study the cavitation dynamics directly within the patch to account for any constructive/destructive interferences in the pressure fields generated by each element (Figure S13, Supporting Information).

In Vitro Permeation Study
In order to evaluate the efficacy of the PZT-D as a sonophoresis device, a permeation study is performed on a porcine model in vitro to evaluate NIA transport in a Franz diffusion cell. A num-ber of studies have prescribed a well-developed protocol to assess the extent of diffusion using high-performance liquid chromatography (HPLC), tape-stripping (TS), and confocal Raman spectroscopy (CRS) and demonstrated correlations between the various methods. [41,[67][68][69] The in vitro permeation setup is shown in Figure 5a, and details of the process are provided in the Experimental Section. A 3% weight/weight (w/w) solution of NIA in phosphate-buffered saline (PBS) is used as the drug formulation. The choice of PBS as the solvent is guided by an experimental study comparing the cavitation threshold for different solvents. Among the formulations of PBS, glycerol, and propane-1-2-diol, PBS displayed the least viscous damping and minimum cavitation threshold ( Figure S7, Supporting Information). The porcine skin is actively treated with the PZT-D for 10 min within the cell at 32°C with the drug formulation in the donor chamber and PBS in the receptor chamber. The PZT-D is maintained at a distance of <2 mm from the surface of the skin using positioning rods in the Franz cell. To monitor the diffusion process, the concentration of the receptor solution is sampled at fixed time intervals throughout a 1 h period while the skin resistance is simultaneously measured using in-cell integrated silver/silver chloride (Ag/AgCl) electrodes. As a control, the same study is set up to allow passive penetration of the solution across the skin without the application of ultrasound.
We anticipate heat generation in the cell due to the dissipation of a portion of the applied ultrasound energy in the fluid medium and the tissue. An initial assessment with an infrared thermal camera revealed that the temperature of the solution in the Franz cell rises significantly above room temperature. (Figure S8, Supporting Information). An in-depth study was subsequently conducted using a calibrated commercial thermocouple logger to measure the temperature profile of the solution in the donor chamber using a PBS solution and PDMS membrane to simulate the skin. In this case, a temperature rise of 1.1°C min −1 was recorded when a 100% duty cycle was applied. Based on previous literature, applying a smaller duty cycle is often preferred in order to provide a cool-down time for the ultrasound transducer and control the heat produced within the system. [70][71][72] Figure  S9, Supporting Information, shows the details of the thermal investigation and Table S1, Supporting Information, lists the temperature rise associated with various duty cycles. While a small duty cycle is desirable from a heating perspective, the ON period should be sufficiently long to allow the nucleation, growth, and collapse of bubbles to take place. A 50% duty cycle with 5000 ON cycles was deemed a reasonable compromise, and was further validated in the Franz cell with the porcine skin. A total local temperature rise of 15°C over 10 min is measured. The average skin temperature is 32°C (91°F), which maps to a maximum temperature of 47°C after 10 min of ultrasound. While this rise in temperature is significant, there are multiple approaches to mitigate this effect, which are presented subsequently in the Discussion.
However, from the perspective of experimental control, it is important to characterize the permeation enhancement provided by the ultrasound as distinct from that provided by the application of heat. A thermal control study is designed and performed using a microcontroller and a resistive heating element to mimic the thermal ramp observed with the ultrasound transducer. The setup of the thermal control Franz cell study is provided in Figure  S5c, Supporting Information.  3). c) Amount of niacinamide (NIA) permeated with the application of ultrasound and heat (thermal control) versus control (passive permeation) for a single PZT-D element. The error bars represent the standard deviation of the measurements (n = 3). d) Custom Franz cell created for in vitro studies with the cUSP device. The porcine skin is clamped between the top (purple) and bottom (pink) chambers, and the cUSP is laid face down (currently shown face up) on the surface of the porcine skin; scale bar: 10 mm. e) Amount of NIA permeated with application of ultrasound versus control (passive permeation) for the cUSP array. *Indicates that the data set is significantly different from the control (p < 0.05). † Indicates that the skin area considered for the calculation has been adjusted to equally discount the neighborhood of skin left uncovered for each study (see Experimental Section). For all cases, ultrasound is applied for 10 min (blue-shaded region). f) Amount of NIA permeated with application of ultrasound versus microneedles. The cUSP array achieves similar NIA permeation within 30 min as compared to the 4 h diffusion curve with microneedles. g,h) Multiphoton confocal microscopy images showing RhB penetration into a vertical section (7 μm thick) of a porcine skin sample for passive diffusion (g) and 10 min ultrasound treatment (h); scale bar: 200 μm. i) Percentage of tissue exhibiting fluorescence across 6 samples of porcine skin (n = 3, passive diffusion; n = 3, 10 min ultrasound). The confocal images are thresholded for fluorescence using ImageJ, details and raw images are given in Experimental Section and Figure S15, Supporting Information.
The trend of skin resistance and amount of NIA permeation from the in vitro permeation study with a single PZT-D element are shown in Figure 5b,c respectively for 10 min ultrasound, thermal control (10 min heat), and passive control. From t = 5 min onwards, the permeation amount of samples treated with ultra-sound is significantly different from that of the passive control samples, while the permeation amount of thermal control samples does not show a significant difference compared to that of passive control for the whole duration of permeation monitoring. The samples treated with ultrasound show a 19.2-fold increase in permeability as compared to the passive control at the end of 60 min (Figure 5c). Application of heat provides a 2.2-fold increase in permeability as compared to the passive control, which indicates that the acoustic streaming and cavitation account for a significant part of the enhancement mechanism provided by the sonophoresis device. The total dose delivered in each case is calculated using a calibration curve shown in Figure S10c, Supporting Information. The absolute dose delivered in a 1 h period is 0.68 ± 0.20 mg cm −2 . This dosage is significantly higher than that obtained with a previous study under pseudo-infinite dose conditions with full-thickness skin wherein negligible permeation was observed within the first 4 h. [73] High NIA transport (0.20-0.80 mg cm −2 ) in a 1 h period in observed in passive diffusion studies with synthetic membranes, [69,73] but this can be attributed to the smaller thickness (7-160 μm) and tailored physicochemical and pharmacokinetic properties of these membranes as compared to the highly heterogeneous 2-3 mm-thick porcine skin. It must also be noted that a smaller enhancement may be observed with human skin. [73] The skin resistance decreases with application of both ultrasound as well as heat, with a larger decrease observed in the case of ultrasound application. While the results are normalized with the initial skin resistance value, large standard deviations are still observed in the skin resistance data (Figure 5b). This can be attributed to the bio-variability of the specimens. The overall thickness of the samples, the relative composition of fat to dermal tissue, and the presence of small hairs and lesions on the skin samples all impact the resistance value measured across the skin sample. However, the 3 test cases (10 min ultrasound, 10 min heat, and passive control) each occupy distinct ranges with almost no overlap. On average, the skin resistance after 60 min decreases by 6.3% in the case of passive control, 22.7% in the case of heat application, and 62.3% in the case of ultrasound.
As a final step to test the device efficacy, we perform a skin permeation study with the 2D array patch in a custom Franz cell. The patch was laid flat on the skin and the donor NIA solution was added to completely fill the cavities around the PZT-Ds until the patch was just submerged to ensure no trapping of air bubbles in the fluid reservoir. The treatment protocol was identical to the single-element study, and the permeation profile was monitored for 1 h. Considering that heat application accounts for a relatively small amount of the overall permeation enhancement resulting from ultrasound, thermal control was not conducted for this series of 2D array permeation experiments.
The cumulative amount of NIA permeation from the in vitro permeation study with the cUSP array is shown in Figure 5e for 10 min ultrasound and passive control. From t = 10 min onwards, the permeation amount of samples treated with ultrasound is significantly different from that of the passive control samples. The cumulative dose delivered in a 1 h period is 0.35±0.09 mg cm −2 for control (passive permeation) and 9.14 ± 4.00 mg cm −2 for ultrasound application with the cUSP. The samples treated with ultrasound show a 26.2-fold increase in permeability as compared to the passive control at the end of 60 min. Since the primary mechanism responsible for the permeation enhancement is inertial cavitation, we expect the effect to be highly localized around the PZT-D. The skin area is adjusted to discount the area of the untreated skin in both setups (see Experimental Section). In the single-element study, the PZT-D covers 25% of the skin directly exposed to the ultrasonic field, while only 12.3% of the total skin area is covered by the 4 PZT-D elements in the larger Franz cell ( Figure S10a, Supporting Information). Once the results have been compensated for treated versus untreated skin regions, the enhancement ratios obtained at the end of 1 h with the patch and with the single element are comparable at 26.2-fold and 19.2-fold respectively (Figure 5d). The higher enhancement ratio provided by the cUSP can be attributed to the larger acoustic pressure generated within the PDMS-bounded cavity walls. We therefore demonstrate the feasibility of the patch design and the PDMS encapsulation in allowing undamped application of ultrasound energy to the skin over larger areas in a conformal, unobtrusive form factor.
To explore how our IFS device performs compared to other permeation enhancement techniques, we conducted in vitro drug permeation studies with microneedling treatment, which is the most common and widely accessible transdermal permeation enhancement method. The porcine skin was pre-treated with a commercial derma roller (200 μm length) and the NIA permeation profile was monitored for up to 6 h (see Experimental Section). The cumulative amount of NIA permeation from the in vitro permeation study with the cUSP and microneedle is shown in Figure 5f. From t = 10 min onwards, the NIA permeation of samples treated with ultrasound is significantly different from that of the passive control samples, while that of the samples treated with microneedle remains similar to that of the control throughout the 1 h permeation monitoring period. The cumulative dose delivered in a 1 h period is 4.00 mg cm −2 for microneedle treatment and 9.14 ± 4.00 mg cm −2 for ultrasound application with the cUSP. The results show that a 10 min ultrasound application with the cUSP is able to deliver 20.8-fold of NIA into the skin compared to microneedling treatment in a 1 h period, and the drug permeation enhancement may continue to increase over time. From Figure 5f, we observe that the cumulative dosage obtained with ultrasound treatment at t = 30 min (2.94 mg cm −2 ) can only be achieved with microneedle treatment at a much later time point, t = 240 min (2.88 mg cm −2 ). The permeation rate resulting from ultrasound application is approximately 13.2 times larger compared to that with microneedle treatment (slope of 0.15 and 0.01 mg cm −2 min −1 , respectively). We therefore demonstrate the superior drug permeation enhancement provided by ultrasound treatment with the cUSP as compared to other permeation enhancement techniques.
As the fate of the drug under consideration is in the skin itself (within the viable epidermis and dermis layer), we further performed a short study to ensure the successful retention of the drug within these layers. 1 cm 2 pieces of pig skin were fixed in a diffusion cell with 0.75 mL of standard stock solution (1 mg mL −1 in PBS) of Rhodamine B in the top chamber, and PBS in the bottom chamber to maintain skin hydration. The same ultrasound treatment was applied at 50 V, 50% duty cycle for a 10 min duration. The cells were allowed to stand for 24 h at 32°C. A crosssection of the skin was analyzed using a multi-photon confocal microscope. Details of the study, histology, and imaging protocols are given in Experimental Section. Figure 5g-i shows the fluorescence response in the tissue for the treated and untreated skin samples. At the end of 24 h, the dye penetrates the epidermis both in the case of passive diffusion as well as ultrasound as can be seen in Figure 5g,h. However, there is a significant increase in dye permeation to the dermis in the case of ultrasound as can be visually appreciated in Figure 5h and Figure  S15b,c, Supporting Information. The percentage of tissue exhibiting fluorescence is plotted in Figure 5i. Skin treated with ultrasound exhibits a 3.5-fold increase in dye concentration as compared to passive permeation. The results are significant with a p-value of 0.002.

Discussion
Bulky and power-consumptive equipment, [74] long exposure times, and high variability between the permeability enhancement observed in vivo and in vitro [32] have been cited as the main challenges to the successful deployment and adoption of transdermal drug delivery systems. Recently, wearable system concepts have been proposed that demonstrate good performance in vivo [48] using HFS but few studies have characterized the working physics of the enhancement mechanisms required to guarantee repeatability. In this study, we use a combined theoretical and experimental approach to characterize and precisely control the ultrasound transducer deflection, device geometry, resulting pressure field, and heat generation within the medium to induce cavitation in a repeatable manner using IFS. We characterize the threshold and strength of inertial cavitation using its acoustic subharmonic and broadband noise signature. We determine the threshold as the driving voltage for the cUSP elements which can be directly set and modulated by the user to linearly ramp the strength of cavitation. We use high-speed imaging to validate the nucleation density and size of bubbles induced, and calculate the bubble streaming velocity close to the PZT-D surface. This acts as a quantitative indicator of the convective effect provided by the ultrasound system, where previously these effects are only qualitatively or theoretically discussed. [23,75,76] Finally, we demonstrate the efficacy of the IFS device in vitro by showing a 26.2-fold improvement in niacinamide transport and decreased skin resistance in a porcine Franz cell model. We present an innovative patch design that allows consistent positioning of the ultrasound transducer and provides a controlled-volume cavity for sustained cavitation activity. As discussed in Table S1, Supporting Information, we demonstrate high enhancement ratios comparable to those obtained with large commercial, hand-held ultrasound probes with a remarkably compact and simplified ultrasound system. To the best of our knowledge, this is the first conceptual demonstration of a wearable IFS system for effective and operator-independent transdermal delivery of cosmeceuticals outside of a clinical setting.
In this system, the cUSP delivers treatment to the areas of the skin directly below the 2D array of the PZT-Ds (each 0.8 cm 2 ). This is an intermediate-zone in terms of spatial targeting as compared to highly localized treatments such as hypodermic needles or microjetting [77,78] (10-100 μm jet/needle diameter) and uniform treatments over large areas of the skin such as with microneedle patches or Dermarollers. [79] For transdermal drugs destined to take action within skin (epidermis and dermis layers), the cUSP achieves superior performance (13.2-fold larger) compared to 200 μm microneedles. While longer microneedles or pulsed microjets may enable higher drug penetration depths and volumes, there is a trade-off with pain, invasiveness, and bruising that makes them less desirable than the cUSP. Further, needle-free microjet injections require sophisticated microfluidic devices and high-power pulsed laser setups, [80] making them unsuitable for wearable applications. In this sense, the cUSP is a realization of a compact, wearable interface for effective and pain-free drug delivery to targeted zones of the skin. In future iterations, this unique capability of the cUSP can be leveraged with computer vision to extract texture and coloration information from a 3D scan of the users' face and provide auto-placement and packing of the PZT-Ds to provide targeted treatment to spots or scars. Such customized embodiments of the cUSP can then be used to conduct in vivo studies to study disease progression in human subjects over a span of months.
An immediate iteration possible on this design involves replacing the 2 mm-thick PZT-Ds with thinner discs (<1 mm thickness). As the radial mode of operation of the disc has been selected, changing the thickness of the disc will cause a negligible shift in the radial resonance mode. This reduces the overall thickness while allowing better flush contact with the skin. The conformability and fit of the patch to curvilinear tissue such as the cheek can be characterized using FEM simulations for initial predictions on stretchability and radius of curvature along various axes, and can be experimentally evaluated using a pressure-sensitive tape between the patch and skin. The design and mechanical properties of the piezoelectric components and polymer encapsulation should be tailored for better biomechanical integration with soft facial skin. [81,82] Further, dual-frequency ultrasound [83] in conjunction with heavy matching layers [52] can be applied alternatively to nucleate, grow and collapse the bubbles using the distinct effects achieved with each mode as seen in Video S5, Supporting Information, and improve the abrasive and convective action of the bubbles. Finally, the issue of heat generation within the patch can be addressed through various routes: 1) shortening the duty cycle or total ultrasound application time, and 2) adding heat-sinks in the form of a metallic or conductive polymeric mass on the backing of the PZT-Ds to dissipate the heat to the surrounding air.
This work makes an important contribution through the rigor of the framework presented for end-to-end design and performance characterization of a rigid-flexible system of "wearable ultrasound" devices. There is a need for a prescriptive framework for performance characterization within which to streamline efforts for material discovery and device engineering of ultrasoundbased drug delivery devices. In this study, we present the complete gamut of bi-modal characterization methods (simulation and experiment), fabrication techniques, and in vitro testing outcomes of a first-of-its-kind prototype of a wearable intermediatefrequency sonophoresis device. The learnings bear the potential to greatly accelerate material exploration and optimized design of all components of the system; the rigid ultrasound transducer, polymeric substrate, and fluid coupling medium. With further miniaturization of the transducer and the electronics module, we envision this device providing a game-changing alternative to oral and needle-based delivery of intramuscular small-molecule (<500 Da) drugs. [84]

Experimental Section
Piezoelectric Transducer Coating (Single-Element): The surface of the piezoelectric disc (each weighing 1.2 g, Steiners and Martins, Inc.) was cleaned using isopropyl alcohol and DI water, and a small piece of copper tape (thickness = 75 μm) was used to make flexible connections to a shielded coax cable. The device was then encapsulated using a dip-coating method. First, the PZT-D was dipped in Omnicoat (Kayaku Advanced Materials, Inc.) and baked at 150°C for 1 min. It was then dipped into SU-8 2050 photoresist (Kayaku Advanced Materials, Inc.), followed by a ramped bake from 65-95°C for 1.5 h to promote photoresist reflow and prevent the build-up of thermal stress and wrinkles in the coating. A flood UV (405 nm) dose of 350 mJ was then applied to the coating on all four sides to cure the photoresist, followed by a hard bake at 95°C for 15 min. Please note that this coating process was applicable to the single-element PZT-Ds used in the measurements for electromechanical, pressure, and cavitation characterization of individual transducers. The PZT-Ds in the 2D array patch were coated directly in PDMS to facilitate better bonding of the coated PZT-D with the PDMS base substrate. [85] The impedance of the piezoelectric device was monitored in air before coating, and in air and water after coating to estimate the damping and ensure water-tightness.
Fabrication of 2D Array Patch: The bulk piezoelectric discs were dipcoated in poly(dimethylsiloxane) (PDMS, Sylgard 184, Dow Corning) and cured at 100°C for 1 h. The electrodes were covered with a small piece of Kapton tape that was peeled-off after the dip-coating to expose the pads for connection to the circuit. For electrical connections to the device, standard photolithography was used to fabricate the copper interconnects on a 36 μm thick copper foil (EMI Copper Foil Shielding Tape 1181, 3 M). Several drops of AZP4620 photoresist (Microchemicals GmbH) were dispensed onto the copper foil attached on a silicon wafer and the whole substrate was spun at 1500 rpm for 45 s using a spin coater (PWM50, Headway Research). The substrate was then baked at 95°C for 2 min followed by 110°C for 2 min on a hotplate (VWR International). To define the serpentine pattern, the substrate was exposed to UV at 403 nm under a transparency mask for 15 s, developed using the AZ400K developer (AZ Electronic Materials), and baked at 120°C for 5 min. Next, the substrate was placed in a bath of copper etchant (Type CE-100, Transene Company Inc.) at 120°C for 1 h to remove the exposed copper. Last, to dissolve the underlying adhesive, the substrate was soaked in a bath of AZ400T stripper (AZ Electronic Materials) at 120°C until the serpentine pattern was lifted off. The electrode was then dried off and ready for subsequent assembly. To assemble the array, the coated PZT-D elements were positioned in a 3D printed (FDM Prusa i3 MK3S+) poly(lactic acid) filament mold with a 20 mm distance between adjacent elements. The serpentine electrode was then positioned in the mold and attached to the electrodes on the back side of the PZT-D elements with solder paste. Copper wires were then soldered to the top side of the copper foil for electrical connection. Part A and Part B of PDMS were mixed at 10:1 weight ratio thoroughly, degassed in a vacuum chamber, and poured into the mold to a thickness of 4 mm to completely encapsulate the PZT-D elements and the serpentine interconnects.
Electromechanical Characterization: An impedance analyzer (E4990A, Keysight Technologies) was used to characterize the electrical impedance spectrum. A 3D laser doppler vibrometer measurement system (MSA-500, Polytec) was used to measure the deflection of the device. A periodic chirp from 200 kHz to 1 MHz was applied to the device, and the PZT-D surface was scanned in 5 separate regions of 1 mm 2 sections radially outwards from the center of the device. The images were then stitched in a postprocessing software to provide a complete map of the displacement profile corresponding to the axisymmetric element from the FEM simulations (Figure 2d,e and Videos S1 and S2, Supporting Information).
Acoustic Pressure Characterization: To evaluate the acoustic pressure profile, the PZT-D was immersed into a tank (49.5 × 25.0 × 29.0 cm 3 ) filled with DI water with a hydrophone (HGL-0200, Onda Corp.) pointed directly along the central axis of the device. A burst mode sinusoid was used to determine the distance between the hydrophone and device, and a continuous sinusoidal signal was used as the excitation signal. A power amplifier (Model 2100HF, Trek) was used in conjunction with a function generator (Model 3390 Arbitrary Waveform Generator, Keithley) to achieve the desired voltage output. In order to avoid standing waves, the water tank was confined at both sidewalls by acoustic absorbing panels that provide a high echo-reduction for frequencies down to 120 kHz (AptFlex F28, Precision Acoustics). The generated ultrasound pressure was then characterized with the hydrophone at different distances away from the PZT-D (1, 2, and 3 cm) using a spectrum analyzer (NanoVNA-H, Seesii). The experimental results were filtered using a low pass filter in MATLAB.
Ultra-Subharmonic Detection: The PZT-D was suspended in a glass water tank (dimensions 26.0 × 26.0 × 26.0 cm) with an acoustic absorber (AptFlex F28, Precision Acoustics) placed diagonally through the tank to reduce ultrasound echos. A needle hydrophone (HNC-0400, Onda Corp.) was placed at a distance of 4.1 mm from the PZT-D along the axis perpendicular to the surface. A constant spacing between the PZT-D and the hydrophone was maintained by staggering the hydrophone response from the primary driving signal by a constant time of 2.8 μs, which corresponds to a distance of ≈4.1 mm (assuming a speed of sound of 1480 m s −1 in water). (Refer to Figure S12, Supporting Information) The tip of the hydrophone was aligned just outside the radial periphery of the PZT-D. This hydrophone positioning allowed the pick-up of the small signals associated with spherically propagating pressure fields produced by cavitation events without being overwhelmed by the primary vibrations of the PZT-D manifest in the direct field of the PZT-D surface. The integrity of the hydrophone signal was maintained by reducing the electromagnetic coupling in the system by introducing a 15 coil differential transformer from the voltage-driving source. An oscilloscope (Picoscope 6) was used to capture an FFT of the hydrophone signal (5 MHz bandwidth, flat-top filter). The magnitude of the ultra-subharmonic of the primary frequency (3f/2) was preferred in the measurement for its better signal strength, as the sensitivity of the hydrophone drops off significantly below 1 MHz. The data was imported as .csv and analyzed using a custom MATLAB script to characterize ultra-subharmonic signal amplitude and broadband noise with increasing voltage. For each device, the average FFT spectrum was obtained by averaging over the frequency response of 32 frames of timeseries data. The broadband noise was calculated across a wide frequency band (208 to 1144 kHz) by excluding the 40 data points around the peaks of the main driving signal, its subharmonics, and harmonics. The data was analyzed using OriginLab software. Results were presented as the mean ± standard deviation (SD) for n = 3 devices. The one-way analysis of variance (ANOVA) with Tukey's post hoc test was used to analyze data between groups. The Brown-Forsythe test was used to confirm homogeneity of variance. Statistical significance was assumed when the p value was less than 0.05.
High-Speed Imaging: A coated PZT-D was suspended in a large water tank (dimensions 26.0 × 26.0 × 26.0 cm) lined with acoustic absorber panels. Phantom v2512 high-speed camera (Vision Research Inc.) was used with two 2× teleconverter lenses (Nikon, Tamron) to achieve up to a 4× magnification of the imaging zone close to the PZT-D surface. A 150 W halogen fiber optic illumination source (AmScope) was used to provide sufficient illumination for high-frame-rate (10 000 fps) and low exposure image capture. The light source was sealed in a plastic ziploc bag and submerged under the water to provide sufficient illumination close to the PZT-D. A millimeter scale was placed at the focal plane to calibrate the pixel size in the image for bubble size and velocity calculations. The footage was recorded using the custom PCC Software and exported to a .tif stack for post-processing and analysis on ImageJ. For bubble size characterization, the frames were visually inspected to select in-focus bubbles, and a circle and bounding box were used to measure the bubble diameter. A total of 60 bubbles were selected across 3 frames and a histogram of their radii was provided in Figure S11, Supporting Information. For the bubble count analysis, a single stream of bubbles was manually selected across a substack of 250 frames. The total pixel area under analysis was kept constant at 12 000 ± 500 pixels. Contrast enhancement was used where appropriate and the figure was then thresholded to select the bubbles and apply false color. A watershed technique was applied to separate any lumped bubble outlines. A particle analyzer was then used to count the total number of bubbles within a certain size/circularity range. This was then divided by the total number of frames to give the average bubble count per frame. For the velocity measurements, a single bubble was tracked from its origin at the PZT-D surface until it reaches the first pressure node at a distance of around 1 cm where it oscillates more or less in place. A straight line was drawn and measured using the calibration scale between these points and the distance was divided by the number of frames and the duration of each frame (20 ms) to yield the velocity in m s −1 .
Skin Preparation and Processing for In Vitro Permeation Study: Porcine ear skin (Animal Technologies) (shaved, cut in pieces of approximately 3 × 4 mm, frozen at −20°C) was placed in a ziploc bag and thawed in water at room temperature for 15 min. The skin was then hydrated in phosphatebuffered saline (PBS, Corning) for 15 min. The thickness of the skin was measured with a caliper just before it was assembled in the Franz cell, and the edges of the skin outside the cell were trimmed to prevent evaporation during the experimental process.
Franz Cell Setup for In Vitro Permeation Study: A custom fixture was designed and 3D printed using poly(lactic acid) for the Franz Cell study. The top and bottom halves of the Franz Cell incorporate holders for 4 mm Ag/AgCl electrodes (BMD-4, BioMed Products Inc.) which were suspended 1 mm above the surface on either side of the skin once the cell was assembled. Once filled, a conductive PBS pathway was created across the skin and the resistance can be calculated using a voltage divider circuit. A function generator (400 3A, BK Precision) was used to apply a 10 Hz sinusoidal signal (100 mV peak-to-peak) to the skin and the following equation was used to estimate the resistivity.
where V out is the voltage at the input of the measurement circuit, and V r is the voltage drop across a known resistor R k . The Franz cell was assembled with the porcine skin tightly clamped between the donor and receptor chambers and placed in a hot water bath maintained at 32°C. A magnetic stirrer was used to agitate the solution at 100 rpm. The bottom cell was filled with PBS as the receptor solution. The choice for the solvent in the donor solution was selected based on experiments to quantify the cavitation thresholds for different solvents ( Figure  S7b, Supporting Information). The resonant frequency shift (as measured from the impedance curve) as well as the cavitation threshold was the lowest for PBS. Thus, a 3% w/w NIA/PBS was used as the donor solution. The single-element coated PZT-D was positioned at 1 mm from the skin surface using a fixture inside the Franz cell and a pulsed sinusoidal signal (f r , 50 V pp , 50% duty cycle, 5000 cycles ON/OFF) was applied. 200 μL samples were collected from the receptor at designated time intervals (t = 5 min, 10 min, 15 min, 30 min, and 1 h) and replenished with an equal amount of PBS.
For the thermal control study, an infrared thermal camera (FLIR C2, Teledyne FLIR LLC) was used to monitor the temperature rise in the Franz cell due to heat generation by the ultrasound device. An in-depth study was conducted using a calibrated commercial thermocouple logger (EL-USB-TC-LCD, EasyLog) to measure the temperature profile of the solution in the donor chamber using a PBS solution and PDMS membrane to simulate the skin. The same temperature profile was then applied using a resistive heating element (10 W, 5 Ω) and microcontroller (SAMD21E, Atmel) to mimic the thermal ramp observed with the ultrasound transducer.
To demonstrate the applicability and efficacy of the 2D array patch, a larger custom Franz cell fixture was designed and 3D printed using poly(lactic acid) ( Figure S10b, Supporting Information). The experiment protocol was identical to the Franz cell permeation study using the singleelement PZT-D. It was to be noted however that the power used to drive the 2D array patch (4 PZT-D elements) ranged from 3.6-4.8 W, while the power used to drive a single PZT-D element ranged from 1.05-1.3 W.
Dye Permeation Study: 1.5 cm × 1.5 cm pieces of skin were prepared as per the protocol given above and fixed in a diffusion cell ( Figure S13a, Supporting Information). The bottom chamber was filled with PBS to keep the skin hydrated on both ends throughout the study. A standard solution (1 mg mL −1 ) of Rhodamine B (RhB, 479 gmol −1 Sigma Aldrich) was prepared in PBS (Corning). 750 μL of the dye solution was pipetted into the top chamber. N = 3 cells were maintained for passive diffusion as controls, and were allowed to stand on the hot-plate for 24 h. N = 3 cells were treated with 10 min of ultrasound at 50 V pp , with a 50% duty cycle. A PZT-D coated in PDMS was suspended 1 mm from the skin surface and submerged into the dye solution in the top chamber. After 24 h, the cells were disassem-bled and the skin surface was cleaned with a cotton swab to remove excess dye from the surface. The skin samples were then flash-frozen in an embedding matrix (Tissue-Tek O.C.T. Compound, Sakura Finetek) at −80°C in a dry ice bath. A cryotome was used to process 7 μm-thick slices of the skin which were fixed on glass slides and imaged under a multiphoton confocal microscope (Olympus FV1000) using an 840 nm laser. ImageJ software was used to characterize the fluorescence response in the images. A binary thresholding operation was performed on both the transmitted optical image and the fluorescence image to select the regions of interest for the tissue structure and the fluorescing tissue respectively. An intersection of these pixel sets (using an AND operation in the ROI Manager) gives the fluorescence response for each sample. This set was divided by the number of pixels in the tissue structure set to yield a normalized measurement of the percentage of fluorescing tissue for each sample. The images from 6 slides (n = 3 passive diffusion, n = 3, 10 min ultrasound) are given in Figure S13b, c, Supporting Information.
Data Analysis for the In Vitro Permeation Study: To quantify the amount of NIA permeation, high-performance liquid chromatography (HPLC) was performed on a 65 μL sample using a Series 1100 HPLC machine (Hewlett Packard). The samples were run in a 20:80 methanol:water isocratic column for 10 min, with the NIA peak emerging at 2.1 min as established by a pure NIA solution as a baseline ( Figure S10a, Supporting Information). The HPLC data was recorded in amount (unit: mAU) versus time (unit: s), and the height of the peak (unit: mAU) represents the amount of NIA present in the solution. To quantify the amount permeated in standard SI units, a standard curve was established by performing a serial dilution of NIA/PBS samples with known concentrations ( Figure S10c, Supporting Information). The amount of NIA permeated was then converted to units of μg cm −2 using the standard curve. The result was scaled with the volume ratio of the receptor chamber to the extracted sample, the initial skin resistance (with the largest value normalized to 1), and the actively treated skin area in each case. In the single-element study, the PZT-D covered 25% of the skin directly exposed to the ultrasonic field, while only 12.3% of the total skin area was covered by the 4 PZT-D elements in the larger Franz cell ( Figure S10b, Supporting Information). Thus, the amount of NIA permeation with ultrasound treatment was calculated by subtracting that of control/thermal control and dividing by the skin area directly exposed to the ultrasonic field. The data were analyzed using OriginLab software. The results are presented as the mean ± standard deviation (SD). One-way analysis of variance (ANOVA) with a Tukey's post hoc test was used to analyze data between groups (ultrasound vs thermal control vs passive control for PZT-D and ultrasound vs passive control for array). The Brown-Forsythe test was used to confirm homogeneity of variance. Statistical significance was assumed when the p value was less than 0.05.
Finite Element Method (FEM): The commercial finite-element analysis (FEM) software COMSOL Multiphysics was utilized to simulate and optimize the mechanical performance and pressure characterization of the device. The tetrahedron elements were used in the solution with adaptive meshing convergence. The element size was tested to ensure the convergence and the accuracy of the simulation results. The default material properties included in the material library were used in the simulation. The model was built based on three fundamental physics interfaces of Pressure Acoustics, Solid Mechanics, and Electrostatics. The Pressure Acoustics interface solves the wave equation in the fluid media surrounding the device. The physics interfaces of Solid Mechanics and Electrostatics model the piezoelectric and solid effects to form the new interface of Piezoelectric Devices that correlate stress and strain to electric field and electric displacement via linear constitutive equations. The bottom surface of the transducer was grounded, and the AC electric potential was applied to the upper surface. The loss factor set in piezoelectric featured mechanical loss as an isotropic loss factor. The electric field was parallel to the polarization direction which was the thickness direction, and was defined as E 3 = V∕h where the piezoelectric transducer was subjected to the actuating voltage V along the thickness direction, h. Details of the theoretical basis for the multiphysics approach used in each section are given section-wise below.
Electromechanical Modeling and Characterization of the Device: The total electric current flowing through the transducer was obtained by differentiating the charge density (Q) with respect to time (t) and integrating it over the electrode area. Considering harmonic excitations, electric charge has the form of Q = Q 0 e i t , which yields the current magnitude of |I| = | Q t | = ⋅ Q 0 . The impedance magnitude was further calculated by the general Ohm's law as |Z| = |V|∕|I|. The resonance frequency of a piezoelectric transducer can be calculated by f = N∕t where N is the frequency constant of the piezoelectric material (radial mode: N p and thickness mode: N t ) and t is the dimension of the transducer (radial mode: diameter and thickness mode: thickness).
The relation for balance of body forces can be expressed as − 2 u = ∇·T where is density, is the angular frequency, u is the displacement vector, and T is the stress tensor. Voltage and displacement were related as T = c E S − e T E where S, c E , and e are the strain, elasticity, and coupling matrix, while E is the electric field vector. [86] Figure S2a,b, Supporting Information, shows the displacement profile of the device from the FEM simulation in the 2 dominant resonance modes. The two main resonant frequencies from the displacement spectra in the FEM simulations were 220 kHz and 1 MHz (Figure 2b). The in-plane z-displacement was symmetric with respect to the center radial line (z = −1).
Acoustic Pressure Field Modeling and Characterization: In the numerical model, the boundary condition of plane wave radiation was assigned to the outer surface of the water domain in order for the wavefronts to travel outward from the geometric boundaries with minimal reflection. In order to solve the wave equation in the fluid media surrounding the piezoelectric transducer, the physics interface of "Pressure Acoustics" was included in the numerical model. The wave equation of acoustic pressure in the form of Helmholtz equation was solved by the Pressure Acoustics, Frequency Domain interface as follows: where p, , f, c, and 0 represent the acoustic pressure, angular frequency and frequency of the wave propagation, speed of the wave, and density of water (assumed to be constant), respectively.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.