Bifurcation of ensemble oscillations and acoustic emissions from early stage cavitation clouds in focused ultrasound

The acoustic emissions from single cavitation clouds at an early stage of development in 0.521 MHz focused ultrasound of varying intensity, are detected and directly correlated to high-speed microscopic observations, recorded at 1 × 106 frames per second. At lower intensities, a stable regime of cloud response is identified whereby bubble-ensembles exhibit oscillations at half the driving frequency, which is also detected in the acoustic emission spectra. Higher intensities generate clouds that develop more rapidly, with increased nonlinearity evidenced by a bifurcation in the frequency of ensemble response, and in the acoustic emissions. A single bubble oscillation model is subject to equivalent ultrasound conditions and fitted to features in the hydrophone and high-speed spectral data, allowing an effective quiescent radius to be inferred for the clouds that evolve at each intensity. The approach indicates that the acoustic emissions originate from the ensemble dynamics and that the cloud acts as a single bubble of equivalent radius in terms of the scattered field. Jetting from component cavities on the periphery of clouds is regularly observed at higher intensities. The results may be of relevance for monitoring and controlling cavitation in therapeutic applications of focused ultrasound, where the phenomenon has the potential to mediate drug delivery from vasculature.

observed at higher intensities. The results may be of relevance for monitoring and controlling cavitation in therapeutic applications of focused ultrasound, where the phenomenon has the potential to mediate drug delivery from vasculature.

Introduction
Acoustic cavitation refers to the formation of bubbles and bubble-ensembles in a medium subjected to the pressure fluctuations of an intense sound field. The phenomenon is known to have a pivotal role in applications such as sonochemistry [1], material erosion [2] and acoustic cleaning [3]. Cavitation is often detected acoustically through the use of a hydrophone device, data from which can be analysed for information regarding the nature of the cavitation activity, such as stable versus inertial behaviour ( [4, chapter 4]). A seminal experimental study on this topic by Lauterborn and Cramer [5], reported a 'sub-harmonic route to acoustic chaos', whereby emission spectra with sub-harmonic features to the driving frequency, f 0 of 22.56 kHz, bifurcate to odd sub-harmonics given by n f 0 /3 and n f 0 /8 lines, then n f 0 /4 lines, with increasing drive amplitude. At a sufficiently high amplitude, the spectrum was dominated by chaotic broadband noise associated with inertial cavitation, although distinct spectral features remained apparent. The nonlinear response of single bubbles exposed to periodic pressure fluctuations, has received comprehensive theoretical treatment (see e.g. [6,7]).
Cavitation activity can also be directly observed and analysed with high-speed photography. Due to framing-rate limitations, the majority of the literature available on highspeed observation of cavitation cloud and bubble-ensemble dynamics, has been undertaken within ultrasound cleaning baths or with vibrating horns, typically driven at several tens of kHz and often in standing wave configurations [8][9][10]. Moreover, the stochastic nature of cavitation nucleation limits observations to steady state conditions, when clouds have become established and adopted quasi-stable, periodic behaviour in response to applied insonation.
High intensity focused ultrasound (HIFU) refers to the minimally-invasive application of ultrasound in the MHz range, for the administration of therapy to targeted diseased tissue. Cavitation is a common occurrence at the intensities typically employed and is the subject of a substantial research effort, seeking to exploit its potential to enhance therapeutic effects [11]. This includes inertial activity, which is characterized by broadband acoustic emissions and stable cavitation, identified by harmonic frequency content, for the delivery of therapeutic agents to tissue via extravasation and tissue permeabilization, [12,13], respectively. The latter often employs shelled contrast agent microbubbles, delivered via intra-vasculature injection, to provide cavitation nuclei in vivo. Realization of the potential for HIFU-mediated, cavitationenhanced therapy is somewhat hindered by a deficit in the fundamental understanding of cavitation cloud evolution in MHz propagating focused ultrasound, particularly over the first few hundred cycles of exposure. The introduction of controlled cavitation to tissue for enhanced therapy would further demand a mechanism through which activity can be monitored, to ensure localization to the target region. As acoustic cavitation in HIFU typically evolves extremely rapidly, detection from inception and through the early stages of cloud development is critical, to avoid unwanted collateral damage to surrounding healthy tissue.
We have recently developed a laser-nucleation technique to predetermine the onset of acoustic cavitation activity in HIFU, permitting interrogation at unprecedented spatial and temporal resolution [14]. Here, we report on the observation of bubble-ensemble dynamics during the formative phases of development in water, systematically studied under HIFU exposure across a range of typical therapeutic intensities, for the first time. Furthermore, a passive cavitation detector (PCD) was constructed, based on a series of preliminary highspeed imaging experiments, sensitive at frequencies matching the bandwidth of the observed ensemble oscillations, for parallel acoustical monitoring temporally resolved to the duration of observation. Analysis of cloud response to HIFU insonation is conducted via three approaches: the sub-harmonic frequency content of the acoustic emissions detected, which is compared to model oscillations of a single bubble of selected quiescent radius comparable to that of the cloud, and the response behaviour of the bubble-ensemble dynamics, extracted from the highspeed microscopic imaging sequences.
The results correlate acoustic emissions directly to resolved, early stage cloud dynamics in HIFU. The analysis suggests that the physical size of a cloud can be inferred from the spectra of the scattered field, for a given set of HIFU parameters, via a simple single bubble model.

The sonoptic chamber
The 'sonoptic' chamber is custom-built to accommodate the focused ultrasound field of a 100 mm diameter, spherically focused HIFU transducer, without reflection or scatter [14], impedance-matched to a drive frequency, f 0 = 0.521 MHz. Pressure maps measured with a fibre-optic hydrophone (Precision Acoustics) through the focal, near-and far-fields, in both the sonoptic chamber and a free-field scanning tank (1 × 1 × 1 m 3 ) configuration, confirmed this to be the case. The device is driven with a sinusoidal signal from a waveform generator (Agilent 33250A) passed through a power amplifier (ENI 3200L). HIFU intensities of peak negative pressure (PNP) amplitude between ∼0.7 MPa (the lowest value at which cavitation  clouds can be reliably nucleated) and 1.3 MPa are investigated. The focal region of the HIFU field is contained within a glass cavitation chamber, figure 1, which facilitates optical access for microscopic observation and pulsed laser irradiation, through long working distance objective lenses (Mitutoyo 20 × 0.28 NA and 50 × 0.42 NA, respectively). High-speed imaging is undertaken with a Shimadzu HPV-1 camera, recording 100 frames (312 × 260 pixels) at 1 × 10 6 frames per second (Mfps), with an individual exposure time of 0.25 µs. Illumination is achieved with a flash lamp coupled to a fibre-optic bundle via a condenser lens.
Acoustic cavitation is nucleated via a 6-8 ns 532 nm laser-pulse (Nano S 130-10 frequency-doubled Q-switched Nd:YAG, Litron Lasers) of energy 0.9 ± 0.1 mJ, incident to the focal region of a pre-established HIFU field. Crucially, this is below the breakdown threshold for the host medium of de-ionized water, degassed to an O 2 content below 4 mg l −1 , which avoids the comparatively large plasma-mediated vapour bubble, generally associated with pulsed laserinduced cavitation research [15,16]. A total of 160 HIFU acoustic cycles are generated, with the laser-pulse incident after 80 µs, to allow for transducer 'ring-up' to the required pressure amplitude. High-speed camera operation is triggered to capture a few frames prior to nucleation of cavitation activity, such that cloud development is observed from inception through ∼50 HIFU cycles.

Passive cavitation detector
The PCD is fabricated from piezo-ceramic composite, measuring 9 × 9 × 3 mm 3 , cut with kerfs to improve response efficiency and reduce lateral modes. The large active area of the device provides a high sensitivity to a frequency bandwidth related to its thickness. Silver epoxy provided the acoustic matching and acted as an electrode to the element. Copper tape provided the ground and shielded the device, particularly from the Q-switch of the pulsed laser. A second electrode is soldered to micro-coaxial cable and isolated from the shielding. A photograph of the device is provided, figures 2(a) and (b). The sensitivity of the device was determined by positioning the PCD in a water tank at ∼1 m from the field of a HIFU field, in degassed water. A transducer (ExAblate 2000, InSightec Ltd) was driven at a very high power to generate cavitation, which could be observed visually and produced a characteristic 'fizzing' sound in the audible range. A 1 ms signal was recorded several times for analysis in MATLAB. Three recordings in the frequency domain were averaged and the 1.17 MHz driving frequency filtered out. In this manner a sensitivity characteristic, as depicted in figure 2(c), was obtained following the boundary of the acoustic spectra.
During an experiment the PCD is positioned within the sonoptic chamber, a few millimetres from the acoustic focus, connected to an oscilloscope (Agilent MS07104A). Data is recorded at a minimum of 1 GS s −1 , saved in .bin format on an USB stick and transferred to a PC for subsequent analysis.

High-speed observations
Figures 3(a) and (b) are sample high-speed imaging data illustrating cloud behaviour in response to HIFU of PNP = 0.72 ± 0.13 (instrument error, according to manufacturer) and 1.04 ± 0.18 MPa, respectively. These are sequential images acquired over a duration of 12 µs, approximately 75 µs following laser-nucleation, during which the cloud has become established and initiated periodic behaviour. Slight upward translation is attributable to the acoustic radiation force of the HIFU insonation, as buoyancy is negligible over the timescale of the observations. Inspection indicates that the quiescent component bubble radius is comparable for both clouds, as expected for acoustic cavitation at a given driving frequency [4,11,17].
Both coalescence of component bubbles during expansion, and fragmentation after collapse, within the clouds is observed. The latter is the mechanism by which the number of constituent bubbles increases, at a rate dependent on the intensity of the HIFU field. Several tens of µs following laser-nucleation, quasi-spherical breathing mode oscillations for the bubbleensembles, closely related to the dynamics of the individual component bubbles, are apparent. The effect is particularly evident towards the latter stages of the movie representation of the high-speed sequences sampled for figures 3(a) and (b), available as supplementary material (from stacks.iop.org/NJP/15/033044/mmedia) 5 . It is well known that oscillating bubbles exert either mutually attractive or repulsive forces via coupling to the radiated acoustic field of other 5 Supplementary material: high-speed sequence in movie format (at full image size acquired by camera). bubbles in the vicinity, or secondary Bjerknes effects [18]. However, due to the proximity of the component bubbles within the clouds generated, figure 3, we attribute the ensemble dynamics to the combined effect of the constituent bubbles, through the physical action of individual bubble oscillations. For example, the collapse of any individual bubble will act to draw its nearest neighbours closer, which when extended to all the bubbles within the population oscillating in phase, generates an overall compression for the cloud.
Figures 3(a)(iv), (viii) and (xii) depict consecutive compressive phases, for the cloud at lower PNP, with approximately one HIFU cycle propagating during the time taken to acquire two high-speed images. This constitutes an ensemble response at f 0 /2, the half-harmonic of the driving frequency, also known as period-doubling [4][5][6][7]. Figure 3(b) is the equivalent high-speed data for a cloud nucleated at higher PNP, whereby a larger bubble-ensemble has developed due to increased levels of fragmentation over the preceding 75 µs. As well as the HIFU intensity-dependent size of the clouds of figures 3(a) and (b), a further notable difference is the additional deflation phases, captured in figures 3(b)(i), (iii), (v), (vii), (ix) and (xi). The full sequence recorded for the cloud at higher PNP, available inmovieformatassupplementarymaterial(seefootnote5),clearlyillustratestheensemble pulsating at more than one frequency.
To quantify the ensemble oscillations, a dark-pixel counting algorithm is implemented to each of the 100 images captured within a high-speed sequence. This effectively yields a 'summed bubble area' variation with time, for every observation of cloud evolution at each PNP investigated. This approach does not explicitly distinguish between ensemble response and constituent bubble dynamics. However, for the high void-fraction clouds being investigated, constituent bubble dynamics and ensemble response are synonymous, as discussed previously. The application of a fast Fourier transform (FFT) to the summed bubble area-time curve from each sequence, is thereby taken to provide the frequency of the ensemble oscillations, at the given HIFU PNP. The resulting high-speed sequence (HSS) spectra for the clouds of figures 3(a) and (b) are presented in figure 4 below, for direct comparison to the frequency content of the acoustic emissions collected from these clouds, specifically. A summary of all experimentally detected acoustic and HSS spectral data is provided in figure 7 below.

Acoustic detection
The acoustic emissions from individual clouds are detected for the duration of high-speed observation, via the PCD device described above. Figure 4(a) is an FFT of a PCD recording of the primary HIFU field at PNP = 1.04 MPa, without the nucleation of cavitation activity (i.e. no laser-pulse incident). The fundamental driving frequency at f 0 = 0.521 MHz is the dominant feature, with a smaller peak at 2 f 0 , the second harmonic. The potential influence of the higher harmonics on the cavitation cloud behaviour is discussed (see the appendix) below.
The green traces of figures 4(b) and (c) are the spectra of the acoustic emissions collected for the clouds of figures 3(a) and (b), respectively. The sub-fundamental peaks (arrowed black) are only detected when cavitation activity is nucleated with a laser-pulse, with structure detail dependent on the PNP of the HIFU field driving the activity. As such, we refer to the frequency of these features as emitted frequencies, f e . The cloud at lower intensity exhibits f e ≈ 260 kHz (7 kHz full-width at half-maximum (FWHM); acoustic data), which corresponds to the f 0 /2 sub-harmonic. respectively. The inset (red, at same frequency scale) represents the ensemble dynamics deduced from a FFT of the dependence of 'summed bubble area' (dark pixel count) on time, throughout the high-speed imaging sequence. An FFT of the RP radius-time curves for a single bubble of selected R 0 (blue dash), under equivalent ultrasonic conditions is also presented, described (see section 4.1). depicted in figure 3(b). The frequency-bifurcation in the ensemble dynamics and acoustic emissions at higher intensities represents a transition into a regime of increased nonlinear response for the cloud-system, a phenomenon which, for single bubbles, has previously received theoretical attention [6,7].
Also depicted in figures 4(b) and (c) are the HSS spectra (red trace inset) obtained from the dark pixel counting algorithm applied to the entire sequence of figures 3(a) and (b), as described. The signal resolution for this analysis is inherently limited by the number of samples available (100 frames per high-speed sequence). Moreover, the frequency of oscillation signal will only become available once the cloud has entered its periodic response phase, typically 20-30 µs following the initial nucleation event. Nonetheless, excellent agreement between the frequency content of the acoustic emissions and ensemble dynamics is apparent. Finally, model spectra (blue dash) from a Rayleigh-Plesset (RP) formulation for single bubble oscillation at selected values of R 0 are also presented, as discussed (see section 4.1) below.

Jetting from peripheral bubbles
The increased nonlinearity of the ensemble dynamics at higher PNPs is underscored by frequent observation of jet, and counter-jet formation [19], from bubbles on the periphery of the clouds, figure 5; a number of examples are arrowed white. Inwardly directed jetting from bubbles that formed at hydrophobic pits, etched in a two-dimensional array on a surface has been reported before [20], albeit under comparatively controlled and idealized conditions. To the best of our knowledge, the observations of figure 5 represent the first at sufficient temporal and spatial resolution to identify jets from bubbles at the periphery of a cloud, that are constituent to it, particularly at a typical HIFU driving frequency. Jetting activity was not been observed at HIFU PNPs < 1.0 MPa.

Single bubble Rayleigh-Plesset model
To investigate the origin of the emitted acoustic frequency content, we implement a RP formulation for a single bubble [4,7,21] given as equation (1). This form of analysis for cloud breathing modes has been undertaken previously for the central region of a 'streamer' in the standing field of an acoustic cleaning bath [7]. A remarkable degree of agreement between the time-varying radius of the cloud, observed at 0.1 Mfps at a driving frequency of 12.96 kHz, and those obtained from an RP formulation was demonstrated. Accordingly, we present our experimental results in parallel with equivalent model predictions, for an R 0 selected such that the features of the frequency spectrum of the single bubble model oscillations match those observed for the bubble-ensembles, through the HSS spectrum and those detected within the acoustic emissions, figures 4(b) and (c): where R is the time-varying radius of the bubble undergoing oscillation and R 0 is as elected quiescent radius, p 0 = 100 kPa is the hydrostatic pressure, p v = 2.33 kPa and κ = 5/3 are the vapour pressure and polytropic exponent of the gas within the bubble. ρ = 10 3 kg m −3 , σ = 72 × 10 -3 N m −1 and η = 0.894 × 10 -3 Pa s are the density, surface tension and liquid viscosity of the host medium, respectively. P(t) represents the HIFU excitation, given the form  (1) is unlikely to deliver spectral features that resemble those experimentally observed. Sample analysis results for variation of host medium viscosity, η and surface tension, σ are available (see section 4.2) below.
The dependence of the model R 0 required to deliver the experimental spectra structure on HIFU PNP, is given in figure 6. The PNP amplitude is both an input parameter to the RP model via the HIFU excitation expression of equation (2), and the experimental factor that determines the rate of fragmentation within the bubble ensemble, and therefore the time-averaged size of the observed clouds.
It is not possible to deduce a quiescent radius for the clouds explicitly from the high-speed sequences. For comparative purposes however, approximations of maximum and minimum cloud radii are also represented, figure 6, using the highest and lowest values of the dark pixel count, averaged over all the high-speed sequences acquired at each PNP investigated. During this process, dark pixels are rearranged into a circle, to homogenize cloud morphologies and effectively assume a void fraction of ∼1, which is reasonable for the ensemble at maximal expansion. The approximation for the minimum cloud radius should not be interpreted literally as some of the collapsed constituent bubbles within the ensemble are likely to be below the imaging resolution for the high-speed camera set-up. Nonetheless, the comparison between the experimental radii approximations and the single bubble model R 0 , which were coupled through the PCD spectra for the acoustic emissions, is compelling, and may be taken to indicate that the frequency content of the acoustic emissions originate from the response of the cloud, acting as a bubble-ensemble. Figure 7 represents an overview of the experimental acoustic and HSS, and RP model, spectral information obtained for each PNP investigated. The bifurcation of f e , both in terms of the emitted acoustic frequencies and the ensemble oscillations, at PNP ≈ 0.78 MPa is clearly visible. The blue dot region represents sub-fundamental RP model oscillation frequencies matched to the experimentally detected f e values, which agrees well with both the bifurcation PNP threshold, and the degree of frequency splitting throughout the bifurcation transition.  Error bars are the standard deviation for each data set, with n 6. The blue dots represent the spectral features above a threshold value derived from the RP model, fitted with selected values of R 0 (see figure 6).

Rayleigh-Plesset robustness analysis
To ensure the sub-fundamental spectral features could not arise from the variation of parameters (including in combination) in the RP model for single bubble dynamics, other than R 0 for a given HIFU PNP amplitude, a robustness analysis is conducted. For brevity, matrices of model spectra are presented for the R 0 and PNPs of interest, through parameter space for surface tension σ , and liquid viscosity η, figures 8(a) and (b), respectively, are presented to demonstrate proof-of-principle. A short discussion on the relevance of the parameters to the observations follows.
The surface tension of a liquid is related to its temperature, such that the range presented corresponds to water from 0 • C (σ = 76 × 10 -3 N m −1 ) to 100 • C (σ = 58 × 10 -3 N m −1 ). Although collapsing cavities are known to generate high core temperatures, including in multibubble configurations [2], the energy is very localized both spatially and temporally, to the location and moment of collapse. We therefore assume room temperature of 25 • C and thus surface tension, σ = 72 × 10 -3 N m −1 , for equation (1).
Increasing the host medium viscosity acts to suppress the amplitude of all spectral features, figure 8(b), as expected. In the extremity, where η = 10 Pa s, the model single bubble predominantly oscillates at the fundamental driving frequency, f 0 = 0.521 MHz, irrespective of R 0 or PNP. The viscosity of tissue is often approximated to that of glycerol, η gl ≈ 1.5 Pa s [22]. However, whole blood has viscosity, η wb ≈ 4 × 10 -3 Pa s [23], a region of parameter space for which the sub-fundamental frequency structure is apparent across the full range of R 0 and PNP reported. As such, the signature acoustic emissions identified may have application for cavitation clouds forming within the vasculature, under HIFU exposure.

Conclusions
We present temporally resolved and directly correlated optical observations and acoustic recordings, of single cavitation clouds developing at a very early stage of evolution in focused ultrasound, for the first time. The frequency of the physical bubble-ensemble oscillations translate directly to frequency content within the acoustic emissions, detectable via hydrophones custom-fabricated for sensitivity over the required bandwidth.
The analysis undertaken does not distinguish between the individual constituent bubble dynamics within the cloud, and the dynamics of the cloud itself. Inspection of figure 3 indicates that constituent bubbles oscillate as part of the ensemble and that the expansion and collapse phases are synonymous for both. The range of quiescent radii inferred from the RP model indicate that the frequency content within the acoustic emissions collected from the clouds originate from a source of radius comparable to that of the cloud, rather than that of the constituent bubbles. Taking the speed of sound in water as 1480 m s −1 , implies a wavelength of λ 0 ≈ 2.7 mm, for the HIFU frequency used in this work. As λ 0 R 0 , the quiescent radius required for the RP model, we conclude that this is a reasonable assumption for the purpose of analysing the acoustic emissions, in terms of scattered primary field. Further work is required to elucidate the role of individual bubbles within the ensemble, and their contribution to the acoustic radiation emitted, particularly through multi-bubble interaction models.
A HIFU PNP amplitude threshold for cloud response transitioning from a stable regime exhibiting f 0 /2, into one of more pronounced nonlinearity with associated frequencies at f 0 /3 and 2f 0 /3 is identified, in terms of the frequencies of the observed ensemble dynamics and in the acoustic emissions detected. The emitted frequencies may be fitted to existing models for bubble dynamics and information regarding the cloud size, relative to the driving frequency and pressure amplitude of insonation, extracted.
This work demonstrates that cavitation clouds can be characterized in terms of signature acoustic emissions, which could potentially be translated for monitoring of cavitation-mediated drug delivery from the vasculature, for focused ultrasound therapy. Tissue represents a much more inhomogeneous and viscoelastic host medium than the one used for this work. However, cavitation activity in blood vessels, from microbubbles delivered intravenously for example, may undergo similar evolution on HIFU exposure. In future work, we will employ the techniques reported to develop cloud manipulation strategies through HIFU duty cycle and intensity modulation, via rapid feedback-control loops detecting at signature frequencies, such as those identified.

Appendix. Assumption of HIFU nonlinearity
The expression used for the HIFU excitation, equation (2), assumes a linear monochromatic wave. It is well known that ultrasound at therapeutic intensities is often nonlinear with potentially strong high frequency harmonic components [24]. In terms of cavitation dynamics, these additional components may result in extraneous oscillations that need to be eliminated as a possible mechanism for the ensemble response observations reported.
Here, we present analysis of the HIFU field for nonlinear components and justify the assumption of linearity for the HIFU expression, by factoring high frequency terms into the RP model at the experimentally determined levels. Figures A.1(a)-(c) are in situ fibreoptic hydrophone recordings of HIFU bursts representative of those used to excite cavitation activity, at PNP = 0.72, 1.04 and 1.29 MPa, representative of the range used in this work. Cursory inspection indicates that the positive pressure amplitude is of approximately the same magnitude as the negative pressure amplitude, commonly taken as an indication of linearity. Figures A.1(d)-(f) are the associated amplitude spectra in the frequency domain, generated by FFT implementation, which reveal slight higher frequency components exist at 2f 0 and 3f 0 , increasing for the larger pressure amplitudes as might be expected.
We demonstrate that the harmonic components do not have a significant effect on the model oscillation dynamics, by modifying equation (2)  where n denotes the harmonic and a n is a scaling factor representing the amplitude of the component. Implementing the RP model with equation (A.1) as the excitation expression, in a HIFU field of PNP = 1.29 MPa (the highest pressure amplitude used, and therefore most nonlinear HIFU generated) for the fundamental frequency f 0 = 0.521 MHz yields the bubbleoscillation spectrum of figure A.2(b). Also included as figure A.2(a), is the equivalent spectrum without the higher frequency harmonic components, as applicable for the linear approximation. Comparison of the spectra indicates that the higher frequency harmonic components of the nonlinear HIFU have no discernible influence on the sub-fundamental peaks in the model bubble oscillations, which match the experimentally detected frequency content of the cloud acoustic emissions. As such, we conclude that the observed frequency splitting in the model oscillation are due to the single bubble itself entering a regime of more pronounced nonlinearity at higher driving pressure amplitudes. We attribute the minimal influence of the harmonics to be due to the relatively small associated amplitudes, and that the higher frequencies are further from resonance with the selected values of R 0 for the model single bubbles.