Photoacoustic (PA) imaging offers several attractive features as a biomedical imaging modality, including excellent spatial resolution and functional information such as tissue oxygenation. A key limitation, however, is the contrast to noise ratio that can be obtained from tissue depths greater than 1–2 mm. Microbubbles coated with an optically absorbing shell have been proposed as a possible contrast agent for PA imaging, offering greater signal amplification and improved biocompatibility compared to metallic nanoparticles. A theoretical description of the dynamics of a coated microbubble subject to laser irradiation has been developed previously. The aim of this study was to test the predictions of the model. Two different types of oil-coated microbubbles were fabricated and then exposed to both pulsed and continuous wave (CW) laser irradiation. Their response was characterized using ultra high-speed imaging. Although there was considerable variability across the population, good agreement was found between the experimental results and theoretical predictions in terms of the frequency and amplitude of microbubble oscillation following pulsed excitation. Under CW irradiation, highly nonlinear behavior was observed which may be of considerable interest for developing different PA imaging techniques with greatly improved contrast enhancement.

## I. INTRODUCTION

Over the past decade, photoacoustic (PA) imaging has emerged as a new modality combining the safety and portability of ultrasound imaging with the specificity of optical imaging and offering excellent spatial resolution.^{1,2} PA imaging exploits the PA effect^{3} whereby absorption of modulated electromagnetic radiation leads to heating followed by expansion and contraction of the absorbing material. This results in the generation of a pressure wave and these acoustic emissions can be detected and used to reconstruct an image. The utility of PA imaging has been demonstrated in applications ranging from blood oxygenation mapping to functional brain imaging^{4,5} and is rapidly being translated into pre-clinical use.

The imaging depth, and consequently also the contrast, in PA imaging are quite severely limited by the attenuation of both the optical and acoustic signals in tissue. To address this challenge and also to improve specificity for molecular imaging applications, a number of PA contrast agents have been developed.^{6–8} These include dyes such as methylene blue or indocyanine green, which give enhanced absorption at specific wavelengths, and also metallic nanoparticles. The latter are of particular interest, as tuning the plasmon resonance characteristics to a particular optical wavelength can lead to effective absorption cross sections several orders of magnitude larger than those typically obtained from a dye.^{9} Concerns about the potential toxicity of nanoparticles have, however, hindered their development for clinical use.

A possible alternative for contrast enhancement is to use microbubbles having a gas core of 1–2 μm in diameter, which can be transient or stabilized by a surfactant or polymer coating. These are well established as contrast agents for ultrasound imaging on account of their excellent acoustic scattering cross section and nonlinear behavior.^{10,11} Phase change microdroplets, which turn into bubbles upon laser irradiation, have been shown to produce a tenfold increase in the amplitude of acoustic emissions compared to those produced by nanoparticles.^{12,13} One disadvantage, however, is the large amount of energy required to vaporize the liquid. As PA imaging advances toward molecular diagnostic^{14} and therapeutic imaging applications,^{15,16} there is a need for more sensitive and multifunctional contrast agents. Such particles should have a specific PA signature and ideally be visible under other modalities such as ultrasound. Recent work has included the development of several PA microbubbles, i.e., stabilized microbubbles coated with an absorbing material that respond to laser excitation.^{17–25} Contrast agents have also recently been used in combination with a modulated continuous wave (CW) laser that offers extended possibilities in frequency domain PA imaging.^{26–30} Here, we propose a PA microbubble agent, following an analogous principle, that is responsive to both pulsed and CW PA imaging. We show that, by coating microbubbles with a thin layer of an appropriately selected optically absorbing liquid, the heat deposited during optical irradiation can initiate oscillation of the microbubbles. This can then lead to high amplitude acoustic emissions and requires less energy than the vaporization of liquid droplets. Microbubbles were fabricated in order to exhibit this behavior. Their response was characterized using a combination of ultra high-speed video microscopy and monitoring of acoustic emissions in order to validate the theoretical modeling. The results not only show that such excitation is indeed feasible in practice, but also that close parallels can be drawn between laser-driven and ultrasound (US)-driven bubbles.

## II. THEORY

The parameters appearing in the equations can be found in Table I and the corresponding physical parameters used for the derivation are displayed in Fig. 1.

Symbol . | Description . | DCM value . | Toluene value . | Units . |
---|---|---|---|---|

cp _{o} | Oil heat capacity | 1190 | 1685 | J/K/kg |

λ _{o} | Oil thermal conductivity | 0.14 | 0.14 | W/m/K |

ρ _{o} | Oil density | 1327 | 867 | kg/m^{3} |

D_{o} | Oil diffusivity | 8.9 × 10 ^{−8} | 1.0 × 10 ^{−8} | m^{2}/s |

μ _{o} | Oil viscosity | 0.43 | 0.59 | mPa s |

$ \sigma o , w $ | Oil $/$ water interfacial tension | 28.3 | 37.1 | mN/m |

σ _{o} | Oil surface tension | 26.5 | 28.4 | mN/m |

Symbol | Description | Value | Units | |

cp _{w} | Water heat capacity | 4184 | J/K/kg | |

λ _{w} | Water thermal conductivity | 0.61 | W/m/K | |

ρ _{w} | Water density | 1000 | kg/m ^{3} | |

D _{w} | Water thermal diffusivity | 1.5 × 10 ^{−7} | m ^{2} /s | |

μ _{w} | Water viscosity | — | Pa s | |

R _{i} | Inner radius | — | m | |

R _{e} | Outer radius | — | m | |

x | $ R i R e $ | — | — | |

y | $ R i , eq R i , 0 $ | — | — | |

V $ \u2009 w 0.1 $ | Heated water volume | — | m ^{3} | |

V _{oil} | Oil volume | — | m ^{3} | |

B $ \u2009 a , a v $ | Average power deposited | — | W/m ^{3} | |

B _{a} | Power deposited | — | W/m ^{3} | |

P $ \u2009 \u221e $ | Atmospheric pressure | 10^{5} | Pa | |

P _{g} | Gas pressure | — | Pa | |

T _{g} | Gas temperature | — | K | |

T_{0} | Room temperature | 293 | K |

Symbol . | Description . | DCM value . | Toluene value . | Units . |
---|---|---|---|---|

cp _{o} | Oil heat capacity | 1190 | 1685 | J/K/kg |

λ _{o} | Oil thermal conductivity | 0.14 | 0.14 | W/m/K |

ρ _{o} | Oil density | 1327 | 867 | kg/m^{3} |

D_{o} | Oil diffusivity | 8.9 × 10 ^{−8} | 1.0 × 10 ^{−8} | m^{2}/s |

μ _{o} | Oil viscosity | 0.43 | 0.59 | mPa s |

$ \sigma o , w $ | Oil $/$ water interfacial tension | 28.3 | 37.1 | mN/m |

σ _{o} | Oil surface tension | 26.5 | 28.4 | mN/m |

Symbol | Description | Value | Units | |

cp _{w} | Water heat capacity | 4184 | J/K/kg | |

λ _{w} | Water thermal conductivity | 0.61 | W/m/K | |

ρ _{w} | Water density | 1000 | kg/m ^{3} | |

D _{w} | Water thermal diffusivity | 1.5 × 10 ^{−7} | m ^{2} /s | |

μ _{w} | Water viscosity | — | Pa s | |

R _{i} | Inner radius | — | m | |

R _{e} | Outer radius | — | m | |

x | $ R i R e $ | — | — | |

y | $ R i , eq R i , 0 $ | — | — | |

V $ \u2009 w 0.1 $ | Heated water volume | — | m ^{3} | |

V _{oil} | Oil volume | — | m ^{3} | |

B $ \u2009 a , a v $ | Average power deposited | — | W/m ^{3} | |

B _{a} | Power deposited | — | W/m ^{3} | |

P $ \u2009 \u221e $ | Atmospheric pressure | 10^{5} | Pa | |

P _{g} | Gas pressure | — | Pa | |

T _{g} | Gas temperature | — | K | |

T_{0} | Room temperature | 293 | K |

### A. Linear theory for the microbubble oscillations

Everywhere, the subscripts *w* and *o* refer to water and oil, respectively. The microbubble oscillations may be shown to obey a Rayleigh-Plesset type equation of the form

where *P _{g}* is determined by the heat transfer into and out of the system, and assumes the following expression for the pulsed laser case:

where *F _{a}* is the thermal energy per unit volume deposited by the laser pulse (in $ J / m 3 $). The more complex interactions relevant for the CW laser irradiation case are described by

where $ B a ( t ) $ is the heat density deposited by the laser ( $ W / m 3 $). Inserting Eq. (3) into Eq. (1) results in

The expression of the coefficients follows directly from Eqs. (1) and (3). Their detailed expression is given in Eqs. (A18) and (A19) in the Appendix. The key features of the microbubble behavior can be further described by linearizing Eq. (4) for small bubble oscillation amplitudes and by Fourier transformation to the frequency domain

Here, *ω* is the angular frequency of the laser excitation. The third-order transfer function [Eq. (5)] reduces to the product of a resonator transfer function (of order 2) and an integrator $ 1 / j \omega $. In the more simple case of pulsed illumination, the integrator disappears since the driving consists in a quasi-instantaneous energy deposition rather than a time modulated power deposition. The gain of the transfer function of Eq. (5) is expressed as

and, similarly,

for the case of pulsed irradiation. Here, *R _{i}* is the internal bubble radius and

*V*

_{oil}is the volume of the optically absorbing liquid surrounding the microbubble. The only other requirement on the encapsulating liquid is that it is immiscible in water. This liquid is hereafter referred to as “oil.”

*R*

_{0}is the initial bubble gas core radius. Equations (6) and (7) show that the overall gain of the transfer function is highest for low density and low heat capacity oils. The major determinant of the oscillation amplitude at resonance, however, is the damping coefficient [directly obtained from the canonical transfer function Eq. (5)] that was found to depend dramatically on the viscosity of the oil chosen

where $ x = R i / R e $ is the ratio of the inner to outer bubble radii. Thus, the damping coefficient decreases with increasing bubble size, with decreasing layer thickness (in competition with the laser energy absorption), and, primarily, with decreasing oil viscosity. This relation emphasizes the necessity for using low viscosity oils.

### B. Undamped natural frequency

The undamped natural frequency of the system follows directly from the expression of the linear resonance curve [Eq. (5)] as follows:

A crucial feature of this expression is the dominance of *R _{i}*, the instantaneous or “hot” bubble radius, in determining the natural frequency. Also, the influence of the oil density, which is present in the denominator, can shift significantly the natural frequency, whereas the surface tension only has a secondary effect.

### C. Scaling function

One challenge in comparing the results from the theoretical model to the experimental data is the difficulty in fully characterizing the system. For example, the variations in the thickness of the oil layer within individual microbubbles is extremely difficult to measure. Some influences, however, such as differences in the absorption of individual microbubbles (i.e., encapsulation efficiency) and potential inhomogeneities in the irradiating laser beam due to its finite size can be compensated for by expressing the amplitude at resonance

and looking at its dependency on thermal dilatation, which can be measured experimentally and individually for each bubble. This results in a scaling function that compensates for the influence of variations in most of the parameters

where the variable $ y = R i , eq / R i , 0 $ is a measure of the thermal dilatation of the microbubble that can be obtained from the optical recordings. This scaling function therefore enables comparison between the recorded individual bubble responses and their predicted behavior. The remaining influences in the bubbles responses such as the influence of an unknown oil layer thickness that varies from bubble to bubble will not allow for a prediction of a single line but of an area wherein the individual bubble responses are expected to fall. This area does not depend on their absorbance nor on the conditions of irradiation. Higher control over the microbubble fabrication process and properties will therefore directly result in narrowing of this area and further refining the understanding of the experimental results.

## III. MATERIALS AND METHODS

### A. Materials

Albumin from bovine serum (BSA), oil red, Nile Red, dichloromethane (DCM), and toluene were purchased from Sigma-Aldrich (Bornem, Belgium). Perfluoropentane (PFP, 99%) was obtained from Fluoromed (Round Rock, TX, USA). All chemicals were of reagent grade and used without further purification.

### B. Microbubble preparation

In this pilot study, we coated the microbubbles with a selection of two oils, DCM and toluene, that can be readily incorporated into the microbubbles. The properties of the oils are given in Table I. To prepare a saturated dye solution, 50 mg of oil red was dissolved in 1 mL of DCM. The solution was stirred for 10 min and excess oil red removed by allowing sufficient time for separation from the DCM. Then, 10 *μ*L of PFP was mixed with 200 *μ*L of the oil red solution. The mixture was added to an aqueous solution of BSA (1 mg/mL) and sonicated using a probe sonicator (Bandelin Sonorex HD3200, Berlin, Germany) with the tip held in the center of the vessel for 5 s and at the air/liquid interface for 15 s. The prepared bubbles were extracted using a syringe and dispersed in 20 mL of deionized water. To make BSA/toluene bubbles, the same volume of toluene was used instead of DCM, and the process was repeated. The microbubbles were not centrifuged because for these experiments the microbubble radius was a study parameter and so it was desirable to have wide size distribution from which individual microbubbles of particular sizes could be selected.

### C. Scaling

The scaling function given in Sec. II C was applied to the theoretical prediction after correction by a factor $ 1 / 1.7 $ to account for the non-fully stable thermal state of the microbubble after 200 *μ*s. This number originates from considering heat transfers at a finite time in the experiments, whereas a fully stable regime is assumed in the theory.

### D. Confocal microscopy

In order to estimate the thickness of the oil layer that encapsulates the microbubbles, a fluorescently labeled batch was produced. Nile Red dye was chosen for this purpose. Confocal microscopy images were taken of ∼100 microbubbles for each formulation. The thickness of the oil layer was extracted by measuring the 1/e width of the Gaussian fit of the fluorescence intensity surrounding the microbubbles [Fig. 2(b)] using a script written in MATLAB. The thickness of the oil layer for the toluene-coated microbubbles was estimated at 1.77 μm with a standard deviation of 0.37 μm and that of the DCM-coated microbubbles was estimated at 1.92 μm with a standard deviation of 0.32 μm. The bubble radius ranged from ∼1 μm to ∼20 μm and no obvious dependency of the oil layer thickness in the bubble size was found.

### E. Experimental setup

The combined optical and acoustical behavior of the microbubbles was studied by introducing a diluted suspension into an Opticell™ (Nunc, Thermo Fisher Scientific, Wiesbaden, Germany) cell culture device, consisting of two optically and acoustically transparent membranes in a rigid frame. The Opticell was immersed in a water bath at room temperature and laser illumination was provided using a frequency-doubled Nd:YAG pulsed laser (Evergreen, Durham, CT, 150 mJ), which was used for pulsed excitation (pulse width Δ*t* = 8 ns, *λ* = 532 nm) and a CW diode-pumped solid-state laser (2 W, *λ* = 532 nm, Changchun New Industries, Changchun, China) which was used for CW excitation. The laser was focused onto a 40 μm diameter spot within the Opticell through a water-immersion objective [LUMPLFL, 60×, numerical aperture (NA) = 0.9, Olympus, Tokyo, Japan]. The Brandaris 128 ultra high-speed camera^{36,37} was used to take optical recordings (128 images) of individual microbubbles through the same objective at frame rates of ∼10 × 10^{6} frames per second. Illumination was provided by a Xenon flash light (EG and G, FX-1163, Perkin Elmer Optoelectronics, Salem, MA).

For the pulsed laser irradiation experiments, the laser was first fired only in the second of the series of 6 movies (sequence of 128 images) recorded by the Brandaris camera, so that the first movie would provide a reference. The pulsed laser was fired for a total of three times at low energy setting (∼50 mJ/cm^{2}). The laser output before entering the microscope was set close to the minimum (∼3 mJ, calibrated using a Coherent FieldMax_{II} laser energy meter, Santa Clara, CA) and further reduced using neutral density (ND) filters (ND 1.65).

The CW laser output was modulated in time at a chosen frequency by an acousto-optic modulator (AA-optoelectronics, Orsay, France). The CW laser was turned on 233 *μ*s before the start of the high-speed recording and the acousto-optic modulator (AOM) was triggered to transmit the laser light onto the sample 200 *μ*s before the start of the Brandaris recording to allow for the bubbles to reach a quasi-stable thermal state. Each bubble was irradiated five times consecutively with a 100 ms time interval. Toluene- and DCM-coated microbubbles were irradiated with a 0.75 MHz modulated laser beam in order to study the effect of variations in the physical properties of the coating. The DCM-coated bubbles were irradiated by a 1.0 MHz and 0.75 MHz modulated laser in order to investigate the effect of the frequency change.

The acoustic signals emitted by the vibrating microbubbles were collected by a 1 MHz (90% bandwidth) focused broadband transducer (C302, Panametric Olympus, Waltham, MA) located at the bottom of the water tank, then fed into a pre-amplifier and stored using a digital oscilloscope (DPO3034, Tektronix, Son en Breugel, Noord Brabant, The Netherlands). The optical focus and the acoustic focus were co-aligned using a 0.2 mm diameter needle hydrophone (Precision Acoustics, Dorchester, UK). Here, we study bubbles with radii between 2 and 6 μm, which is a relevant size for imaging contrast agents. Such bubbles are also easily imaged under a microscope. The expected resonance frequency for this size range is close to 1 MHz, and this determined our choice of transducer center frequency. Using the bubble size as a study parameter further required the use of a broadband transducer. In the CW laser experiment, the large bandwidth of the transducer allowed for the simultaneous extraction of the harmonic and subharmonic signals. We chose a focused transducer for its advantage in sensitivity.

### F. Data processing

#### 1. Pulsed laser irradiation

The radius-time curves describing the response of the microbubbles to the pulsed laser irradiation were extracted from the ultrahigh-speed movies using a custom written MATLAB script (version R2012a, The MathWorks, Natick, MA) implementing a combination of contrast correction, thresholding, and area detection. This analysis method showed good stability despite the variations in illumination intensity across the recording. A second script was used to extract the properties defining the microbubbles oscillatory behavior, through fitting of the experimental data to the following equation describing the impulse response of a damped oscillator:^{38,39}

Here, *A* is the amplitude of the response, *z* is the damping coefficient of the oscillator, and *ω*_{0} is its eigenfrequency. *t* is the time and *R*_{0} is the initial bubble radius. The script implements an error minimization based on a discrete variation of the parameters *A*, *z*, and *ω*_{0}. Discrete steps of $ 10 \u2212 2 $ for the damping coefficient, 10 nm for the amplitude and 10 kHz for the frequency were chosen. This approach allows for a more robust analysis while ensuring the required precision.

#### 2. CW laser irradiation data

The radius-time curves of the bubble responses from the high-speed recordings were extracted using the same MATLAB script as for the pulsed excitation. The frequency content and corresponding amplitude of the oscillations were then extracted using fast Fourier transform (FFT) analysis. The phase of the oscillations was determined by fitting a sine wave to the microbubble oscillations. Both operations were performed in MATLAB.

The pressure emission curves recorded by the 1 MHz transducer and corresponding to the first CW laser irradiation were all analyzed using the fast Fourier transform function in MATLAB also for the harmonics and subharmonics. The transducer was calibrated in emission by sending a chirp wave to a needle hydrophone (0.2 mm, Precision Acoustics) placed at the focus of the transducer. Then a perfect reflector (stainless steel) was placed at the focus of the transducer to measure the receive characteristics using the same chirp wave.^{40} Finally, the conversion of the time to a bubble radius was performed with the help of a finite difference simulation of the heat transfer using the method described in the Appendix and in Ref. 12. During the CW laser irradiation, the bubble heats up and grows in time following a function $ f ( 1 / \tau ) $ with $ \tau = R i eq 2 / D w $ with *τ* the typical time scale, *D _{w}* the thermal diffusivity of water, and

*R*

_{i}_{eq}the equilibrium radius of the bubble, which is also equal to the thickness of the thermal boundary layer at equilibrium.

^{12}The function

*f*is given by the finite differences simulation. The conversion is finally achieved by using the two reference points given by the cold initial radius of the microbubbles measured in the first high-speed recording (laser off) and the average hot radius of the microbubble after 200

*μ*s of heating measured in the second high-speed recording, i.e., at the first laser irradiation.

## IV. RESULTS

### A. Pulsed irradiation

#### 1. Experimental bubble response

The aim of the pulsed laser experiments was to examine the impulse response of the microbubbles and to extract their oscillatory behavior. Two examples of impulse responses of DCM and toluene-coated microbubbles are displayed in Figs. 3(a) (Ref. 41) and 3(b),^{41} respectively, derived from the ultra high-speed optical recordings. The laser pulse of 8 ns hits the microbubble at *t* = 0. During the pulse, the oil quickly heats up and transfers the thermal energy to the microbubble gas core, thereby initiating the microbubble oscillations. Example frames from the optical recording of a DCM-coated microbubble are shown in Fig. 3(c).

The same experiment was reproduced to obtain 66 individual responses for the toluene-coated microbubbles and 57 for the DCM-coated microbubbles. The corresponding undamped natural frequencies *ω*_{0} for each bubble were extracted by fitting Eq. (12) to the radius-time curve. The results are shown in Figs. 3(d) and 3(e) with a root-mean-square error of 0.27 MHz and 0.17 MHz, respectively. The theoretical undamped natural frequencies were calculated from Eq. (9) and are also plotted in Figs. 3(d) and 3(e).

The choice of the simple impulse response equation [Eq. (12)] is justified, since the time scales of heat deposition, cooling, and oscillations are decorrelated. In practice, this decorrelation implies that the core temperature is oscillating around a hot temperature instantly produced by the laser impulse. As a response to this impulse excitation, the bubble oscillates while cooling down on a much slower time scale. The theoretically computed microbubble undamped natural frequencies agree very well with the measured values for both microbubble populations. The difference between the undamped natural frequencies of both microbubble populations can also clearly be seen in Figs. 3(d) and 3(e).

The toluene-coated microbubbles exhibited a damping coefficient of $ z = 0.148 \u2212 0.1 + 0.15 $ while for the DCM-coated microbubbles it was significantly higher $ z = 0.247 \u2212 0.15 + 0.15 $, which is expectedly much higher than the theoretical value for a bubble free of BSA. Such a result is also well known for acoustically driven bubbles. For the measured range of bubble sizes and oil thicknesses, the oil-coated bubble free of BSA has a damping coefficient from Eq. (12) $ z = 0.0175 \xb1 0.0055 $ for the toluene-coated microbubbles and $ z = 0.0119 \xb1 0.0036 $ for the DCM-coated microbubbles. It is also likely that partial phase change of these relatively volatile oils impacts on the damping coefficient.

### B. CW laser irradiation

#### 1. Experimental bubble response

An example of an optically recorded response of a DCM-coated microbubble to CW laser irradiation is shown in Figs. 4(a)–4(c). The first recording was taken before the laser was on thereby displaying the size of the “cold” bubble.

As indicated in Fig. 4(a) (Ref. 41) and in selected frames, Figs. 4(b) and 4(c), the average size of the microbubble when exposed to the laser light increases significantly. Clear oscillations around this average radius are also visible in Fig. 4(a). These oscillations are a direct consequence of the modulation of the laser intensity. The bubble was further irradiated for 100 *μ*s, i.e., for a total of 300 *μ*s, during which the vibrating microbubble emits an ultrasound wave that was recorded by the receiving transducer positioned at the bottom of the water tank. An example of the acoustic emission from a DCM-coated microbubble irradiated by a 1 MHz modulated laser is displayed in Fig. 4(d). The acoustic trace has been filtered in the Fourier domain to remove the high frequency noise arising from the electronics.

#### 2. Single microbubble resonance

Figures 4(b) and 4(c) show the increase in the microbubble size between two time points. In fact, the microbubbles are growing continuously during the 300 *μ*s during which the laser is on. Thus, in some cases, the microbubble will grow from below the resonant size to above the resonant size thereby describing its own resonance curve. This feature cannot be optically captured as no existing optical system can record the phenomenon with sufficiently high time resolution over such a long duration. The ultrasound transducer, however, can capture these events via the pressure emitted by the microbubble during the entire duration of the laser exposure. Knowing the size of the bubble at two time points from the optical recordings, these acoustic emission curves can be processed in order to retrieve the full resonance curve of a single bubble from a single CW laser exposure. The details are given in Sec. III F.

A sample of such resonance curves is given in Fig. 5 for DCM-coated microbubbles irradiated by a CW laser modulated at frequencies of 1 MHz [Figs. 5(a) and 5(b)] and 0.75 MHz (Ref. 41) [Figs. 5(c) and 5(d)] and toluene-coated microbubbles irradiated by a CW laser modulated at 0.75 MHz [Figs. 5(e) and 5(f)]. The horizontal axis gives the radius of the bubble normalized by its resonant radius. It is clearly visible in Fig. 5 that the resonance is correlated with a strong shift in the phase of the received pressure wave, as expected for a resonating system. It is also clear that individual microbubbles can display significantly different behavior from one another and the resonance frequency of each can vary by as much as 15% from the predicted theoretical value. The mean resonant size measured from the resonance curves was 3.42 ± 0.23 μm, 4.13 ± 0.27 μm, and 4.31 ± 0.53 μm for the DCM bubble excited at 1.0 MHz, the DCM bubbles excited at 0.75 MHz and the toluene bubbles excited at 0.75 MHz, respectively. The theoretical values given in the same order are 2.89 μm, 3.81 μm, and 4.03 μm. Note that this deviation from the theoretical predictions was not observed when the bubble is irradiated with a pulsed laser. Also as observed from the impulse response of the microbubbles, the damping is significantly different for each microbubble as is evident from the various widths of the individual resonance peaks.

#### 3. Microbubble population response

In order to better understand the behavior of the microbubble population as a whole, we now look at the response of each bubble at a single time point given by the optical recording. By doing so, we look at the global behavior of such laser-driven microbubbles. The influence of the variation in absorbance from one microbubble to another can be compensated for using Eq. (11). Using the range of damping coefficients measured from the impulse response of the microbubbles, a theoretical range for the microbubble responses can be calculated. This range is represented by the red area in Fig. 6(a). On the same graph, the response of the microbubbles is represented by blue error bars. They represent the standard deviation per bin of ten bubbles. Overall, the microbubbles show a stronger response than predicted. This is once more probably due to a partial vaporization/condensation of the oil layer, which absorbs part of the pressure variations in the core of the microbubble and thus allows for a larger oscillation amplitude. Despite this discrepancy, the resonant radius and the relative phase shift shown in Fig. 6(b) are in reasonable agreement with the expected behavior.

#### 4. Harmonic and subharmonic generation

The resonance phenomenon shown in Fig. 6 is based on the response of the microbubbles at the fundamental frequency, i.e., that of the laser modulation. Just as in the case of acoustic excitation, the microbubbles also generate at harmonic and subharmonic frequencies due to the intrinsic nonlinearity of their volumetric oscillations. The nonlinearities are known to be further strengthened by the presence of the coating (surfactant, protein, or phospholipid, for example) necessary to stabilize the microbubbles long enough to use them in practice.^{42} Figure 7 shows the amplitude of the first harmonic and of the subharmonic generated by a microbubble as a function of its size. The same bubbles as those depicted in Fig. 5 are used here and the horizontal axis is scaled to the resonant radii of the individual microbubbles.

From the results it appears that the harmonic generation mostly occurs around resonance [Figs. 7(a), 7(c), and 7(e)] with an amplitude up to −10 dB as compared to the fundamental, which makes them practically usable. Furthermore, only a few bubbles did not generate a significant second harmonic. Interestingly, the generation of subharmonic components was not highest from bubbles of resonant size but in the case of DCM-coated microbubbles for radii $ \u223c 10 % $ bigger for the DCM-coated microbubbles [Figs. 7(b) and 7(d)]. In some cases, the subharmonic was actually larger than that at the fundamental frequency, with a maximum of +8 dB. A positive subharmonic generation on a dB scale practically implies that the microbubble experiences a period doubling as observed in a few optical recordings. This effect was not observed for the toluene-coated microbubbles (possibly due to the smaller number of bubbles growing across their resonant size) and could be a result of the higher volatility of the DCM as compared to the toluene.

## V. DISCUSSION

In the experiments described in this study, the exposure lasts for 200 *μ*s before the first recording is taken. Also from the finite difference simulation, we know that a resonant bubble at 1 MHz reaches 73% of the final thermal equilibrium radius after 100 *μ*s and 79% in 200 *μ*s, which represents only a small difference. After a short time ( $ < 100 \u2009 \u2009 \mu s $) the system temperature rises very slowly (due to the spherical configuration), and the bubble grows to 95% of its final equilibrium size in tens of milliseconds. This problem can be addressed theoretically more rigorously by setting the equilibrium radius to be time dependent, and considering its variations to be slow compared to the time scale of the oscillations. Strictly speaking, however, this treatment should be combined with a treatment of the partial vaporization of the fluids that is for now the main imprecision of the proposed model. Further development of the theoretical model must include vaporization and molecular diffusion, which are both relevant at longer irradiation time scales. It is important to note here that the differences observed between the theoretical resonant radius and the measured resonant radius can partly result from the addition of the imprecisions arising when simulating the bubble heating without Rayleigh-Plesset dynamics. This is expected to have a minor effect on the measurement, as compared to the effect of the partial vaporization of the oil.

With regard to the materials used, the liquids chosen for the purposes of these experiments are referred to as oils. It is to be noted that these are not necessarily oils in the chemical sense, and are only required to be immiscible with water. It is also clear that DCM and toluene are not biocompatible liquids, but they offer the advantage of being well characterized and easy to handle. The positive results obtained in this paper provide strong motivation to develop formulations consisting of more biocompatible oils. These oils should also be less volatile. In terms of safety, the microbubbles were seen to be highly responsive at low laser energy: The energy used for the pulsed excitation was ∼50 mJ/cm^{2} well within the range considered safe for medical use^{43} within the biological window (∼1000 nm). The use of CW lasers offers many advantages in terms of costs, ease of use, signal-to-noise ratio, and allows for making use of some rather unique features as demonstrated above. It also requires, however, higher levels of energy deposition by the laser. We estimate the intensity of the CW laser used for these experiments to be between 50 and 500 kW/cm^{2} which, integrated over the irradiation time is too high for safe use. Furthermore, after 200 *μ*s of laser irradiated, <10% of the input energy is stored in the bubble itself, which limits the efficacy of the system. Nonetheless, achieving the measured signal amplitudes with these relatively low powers during the concept phase is a good sign that this technology has the potential to meet the biomedical requirements for the use of CW lasers.

The theory in Sec. II A and the Appendix does not account for the details of the water/oil interface dynamics that is influenced by the stabilizing BSA shell. The model could be improved by including the coating effects. It is, however, more relevant to develop this modification for more common coating materials, such as phospholipids, as opposed to those used in the present test system. Also, the stabilization of these bubbles was assured by the presence of BSA. It is possible that heat induced a denaturation of the BSA while the oscillations would prevent a stiff reticulation. This effect was not investigated and not accounted for since it has a limited interest for this proof of concept. It would, however, become important in a next step where the obtained damping coefficient is explained in detail.

The damping of the oscillations sustained by these BSA-stabilized microbubbles is larger than that predicted for a surfactant-free bubble. It is, however, interesting to note that the presence of a low viscosity oil layer is sufficient, in principle, to significantly decrease viscous damping. This conclusion also applies to the case of ultrasound-driven bubbles, which offers interesting leads for the design of even more efficient acoustic (and PA) contrast agents.

The CW laser in this study was modulated using a sine wave for practical simplicity. In theory, a square wave is expected to induce a slightly higher response and to also generate different harmonic response peaks. A more simple sine modulation containing a single harmonic frequency seemed, therefore, more suited for a proof of concept study where unexpected effects have to be pointed out, but this should be explored further in future work.

## VI. CONCLUSIONS

In this study, we have fabricated two types of oil-coated microbubbles, one with dichloro-methane and a second type with toluene. In both cases, the oil entraps a dye that absorbs the laser light. The bubbles were first irradiated with a pulsed laser and good agreement between the measured natural frequency and the theoretical predictions was found. The non-dimensional damping coefficient was also measured and shown to be higher than that expected for an oil-coated microbubble. There was also considerable variation from one bubble to another. The microbubbles could successfully be driven by a modulated CW laser irradiation and displayed the expected resonance behavior associated with strong and nonlinear acoustic emissions. The microbubble response was seen to vary with repeated laser pulses, with up to a 60% change in response amplitude being recorded for individual microbubbles, which is attributed to a modification of the dye/oil structure within the coating of the bubble. Under CW laser irradiation individual microbubbles also displayed considerable variation in behavior, mostly in their resonant size, with some bubbles growing from sub-resonant to super-resonant size during exposure. Despite this variability, there was still good agreement with the theoretical predictions for the overall behavior of the microbubble population. The microbubbles exhibited stronger nonlinear behavior than expected, manifested in the generation of acoustic emissions at harmonic and subharmonic frequencies. In particular, a sub-population of bubbles ∼10% larger than the resonant size generated significant subharmonics, possibly due to period doubling behavior. These surprising results may be of significant interest for contrast-enhanced PA imaging.

## ACKNOWLEDGMENTS

This project was made possible by the funding of NanoNextNL, a micro and nanotechnology consortium of the Government of the Netherlands and 130 partners project, and the UK Engineering and Physical Sciences Research Council (Grant No. EP/I021795/1). We also warmly thank Gert-Wim Bruggert, Martin Bos, and Bas Benschop for their ongoing and effective technical support.

### APPENDIX: MATHEMATICAL DERIVATION OF THE THEORETICAL MODEL

##### 1. Momentum equation

The Navier-Stokes equation for an incompressible, Newtonian fluid is as follows:

Body forces will be negligible and a spherical symmetry case is investigated leading to

In the simulation, the bubbles will have an oscillation amplitude in the order of micrometers and the frequency will be in the order of MHz. Speeds will therefore be approximately $ 1 \u2009 \u2009 m \u2009 \u2009 s \u2212 1 $ and thus much lower than the speed of sound. For this reason, incompressibility of the liquid is assumed leading to

where $ R \u0307 $ is $ d R / d t $. With this we find

This equation can be written for both the oil layer and the water. When integrating from *r* = *A* to *r* = *B* the term $ \mu \u2207 2 v $ drops out and this gives

Taking the inner bubble radius *R _{i}* for

*R*and integrating over the water domain

Integrating over the oil domain

can be rewritten as

##### 2. Normal component stress tensor

###### a. Over the oil gas interface

where $ \sigma o \xaf \xaf $ is the strain tensor and $ e r \u2192 $ denotes that it is in the *r* direction. $ \delta P 1 $ is the difference in pressure over the oil–gas interface.

where $ u \u2032 $ is the velocity derivative to the radius. *σ* is the surface tension. Knowing $ v = R \u0307 i R i 2 / r 2 $ we also know $ v \u2032 ( R i ) = \u2212 2 R \u0307 i R i 2 / R i 3 $. Therefore:

###### b. Over the oil water interface

Knowing $ v = R \u0307 i R i 2 / r 2 $ we also know $ v \u2032 ( R e ) = \u2212 2 R \u0307 i R i 2 / R e 3 $. Rewriting the equation above then gives

resulting in

##### 3. Combining

We know that $ P g \u2212 P \u221e = P ( R i \u2212 ) \u2212 P \u221e $ because the pressure at the inside of the inner radius of the bubble is by definition in the gas and, therefore, $ P g = P ( R i \u2212 ) $. We can rewrite by adding and subtracting similar terms

Part 1 of Eq. (A6) is defined in Eq. (A4), part 2 is defined in Eq. (A3), part 3 is defined in Eq. (A5), and part 4 is defined in Eq. (A1). Thus, the complete equation is

Rewriting gives

And rewriting further

To reach the modified Rayleigh-Plesset equation

where the viscosity of water *μ _{w}* is temperature dependent following the relation

^{45}

where *T* is the temperature of the water at the water–oil interface. This new Rayleigh-Plesset (RP) equation simplifies to the classic RP equation for a bubble with only one liquid around it when the properties of water and oil are chosen identical.

##### 4. Small variations around equilibrium

In this part, small variations are added to the static solution in order to obtain a simple model describing the simulation result. The static solution assumes the temperature in the gas to be homogeneous. For this to be true for a modulated laser signal, the diffusion time of the heat in the gas should be smaller than the half period of the laser modulation. This can be contained in the following equations:

we make the assumption *a priori* (verified *a posteriori*) that the microbubble resonance frequency is in the same range as its acoustic resonance frequency and using a Minnaert approximation, the temperature in the gas can be considered constant for bubbles smaller than

A diffusion distance estimation can be done for the water to find approximately $ 0.1 \u2009 \u2009 \mu m $ for $ 1 \u2009 \u2009 \u2009 MHz $ frequency.

The change in temperature over time can be described as follows:

with *ρ _{o}* the density of the oil, $ V w 0.1 $ the volume of the first $ 0.1 \u2009 \u2009 \mu m $ of water,

*cp*the heat capacity at constant pressure of the oil,

_{o}*B*

_{0}the maximum power of the laser, and

*B*the absorbed laser power ( $ W \u2009 \u2009 m \u2212 3 $),

with *B _{a}* being the amplitude of the oscillations in the power and $ \omega = 2 \pi f $. Using the initial equilibrium solution

*I* is the average laser power. Filling in $ C 1 o $ and $ C 2 o $ and rewriting gives

We know that

where *T _{g}* can be filled in and

*P*is the total pressure as used in the modified Rayleigh-Plesset equation [Eq. (A10)]

where *P _{g}* is the found

*P*and $ P \u221e $ is the pressure at infinity. Expressing this in the new variables

*P*

_{0}as the pressure at the beginning of the static solution and

*P*as the total pressure in the gas, and filling in

_{g}*T*and the found pressure, this gives

_{g}Organizing for $ R \u0307 , \u2009 \u2009 R \u0307 2 $, and $ R \u0308 $ and, as an approximation, taking all *R _{i}* and

*R*to be $ R i , eq $ and $ R e , eq $ in case they are multiplied by $ R \u0307 , \u2009 \u2009 R \u0307 2 $, or $ R \u0308 $.

_{e}In order to find an equation that does not contain an integral, everything is derived to time.

in which case $ R \u0307 e $ is assumed to be approximately $ R \u0307 i $ and $ T gas , eq + T room $ is now called $ T g 2 $. The term $ 2 \gamma R \u0308 i R \u0307 i $ is of higher order and is therefore neglected. *R _{i}* is expected to act as an harmonic oscillator and will therefore have the shape of

which has the shape of a transfer function

with *G* being the gain, *O* the output, *I* the input, *z* the damping, and *ω*_{0} the angular eigenfrequency. One thing that can be noted here is that this transfer function is of third order where a standard RP equation would be of second order. The expected phase difference in our case is therefore *π* at resonance instead of $ \pi / 2 $ such as in the normal RP equation,

From Eq. (A15) we can find

Therefore, *ζ* can be simplified as

where the equilibrium pressure $ P g , eq $ is the atmospheric pressure plus the Laplace pressure jump over both interfaces

*ω*_{0} is not a function of time so all *R _{i}* and

*R*are now $ R i , eq $ and $ R e , eq $,

_{e}which altogether is an expression for the angular eigenfrequency as a function of $ R i , eq $ and $ R e , eq $. This shows the eigenfrequency is inversely related to the bubble size, but also shows that the oil layer thickness plays a role. The denominator under the square root shows an inertial shift of the resonance curve: Because oil and water have different densities, the thickness of the oil layer influences the mass to be displaced and, therefore, the resonance frequency.

Now to find an expression for the damping. According to Eq. (A27),

## References

^{6}frames per second.