Laser-driven resonance of dye-doped oil-coated microbubbles A theoretical and numerical study

Microbubbles are used to enhance the contrast in ultrasound imaging. When coated with an optically absorbing material, these bubbles can also provide contrast in photoacoustic imaging. This multimodal aspect is of pronounced interest to the ﬁeld of medical imaging. The aim of this paper is to provide a theoretical framework to describe the physical phenomena underlying the photoacoustic response. This article presents a model for a spherical gas microbubble suspended in an aqueous environment and coated with an oil layer containing an optically absorbing dye. The model includes heat transfer between the gas core and the surrounding liquids. This framework is suitable for the investigation of both continuous wave and pulsed laser excitation. This work utilizes a combination of ﬁnite difference simulations and numerical integration to determine the dependancy on the physical properties, including composition and thickness of the oil layer on the microbubble response. A normalization scheme for a linearized version of the model was derived to facilitate comparison with experimental measurements. The results show that viscosity and thickness of the oil layer determine whether or not microbubble resonance can be excited. This work also examines the use of non-sinusoidal excitation


I. INTRODUCTION
Ultrasound imaging is a safe, fast, and relatively cheap imaging modality that offers high resolution images deep inside the human body. However, ultrasound imaging lacks the specificity of other techniques such as photoacoustic imaging. Contrast in photoacoustic (PA) imaging results from variations in the absorption of pulsed or modulated light and the subsequent propagation of sound in tissue. 1 The amplitude of the emitted acoustic signal is unique for every tissue type which makes PA imaging a very attractive clinical imaging modality. A major limitation of PA imaging, however, is its limited penetration depth that restricts the technology to superficial or catheter-based tissue imaging. Diffusive scattering and absorption of both the incident light and subsequent acoustic emissions prevent adequate signalintensities from being obtained beyond a depth of a few mm. 1 A possible solution to this problem is to use contrast agents, in the form of dyes or nanoparticle suspensions, to increase optical absorption at the site of interest and thereby increase the amplitude of the acoustic emissions from that region. 2,3 Metallic nanoparticles in particular have been shown to offer considerable improvement in photoacoustic contrast. 4 By exploiting the plasmon resonance phenomena, their optical absorption at a given wavelength can be much greater than that available from, e.g., hemoglobin. 5 Nanoparticle agents have been designed with multiple shapes and sizes to tune the absorption wavelength and to add specific functionalities. [5][6][7] The biological safety of nanoparticle agents, however, remains uncertain, motivating the scientific community to look for alternatives. It has recently been proposed that even greater contrast enhancement in photoacoustic imaging could be achieved through the use of volatile droplets whose vaporization is triggered by light. 2,8,9 However, the necessary phase change consumes energy under the form of latent heat, that cannot then participate in the acoustic generation, which makes the use of the available energy inherently suboptimal. Stable microbubbles modified with the addition of an optically absorbing coating have been proposed as a potential solution to this problem. 10,11 This is an attractive approach since gas bubbles are well established as contrast agents for ultrasound imaging. [12][13][14][15] In this case, the microbubble oscillations are stimulated by the heating and cooling of the coating and subsequently of the gas core upon optical irradiation.
Previous work has examined the case of a microbubble exposed to a laser pulse. 16 The purpose of that study, was to compare the relative efficiency of different types of PA contrast agents according to their geometrical arrangement of the various components. The effects of heat transfer between the microbubble coating, gas core and surrounding liquid were neglected, for simplicity. In this paper, we show that these effects are in fact important for PA generation by light absorbing microbubbles. We propose a revised theoretical description of a spherical bubble consisting of a gas core surrounded by an optically absorbing layer suspended in an aqueous environment. The absorbing layer considered here consists of an oil in which a dye of specific optical properties can be dissolved. 17 We show that through an appropriate selection of materials, a strong microbubble resonance can be excited. We also study the influence of the different microbubble and laser light exposure parameters upon the amplitude and frequency spectrum of the acoustic emissions and demonstrate the potential for utilizing a harmonic imaging technique to achieve further improvement in imaging sensitivity.

A. Physical problem and analytical derivation
In this paragraph, we give a summary of the theoretical derivation containing the main steps of the reasoning and of the derivation. All details can be found in Appendix A. The microbubble system (Fig. 1) consists of three domains: a gas core, an oil layer and the surrounding liquid and is assumed to remain spherically symmetric. In each domain, three interrelated physical processes are to be evaluated. First, the thermal diffusion-convection problem will define the heat transfer and the instantaneous temperature in each domain. In spherical coordinates, it obeys the relation, where T denotes the temperature field, r the radial coordinate andT ¼ rT. The density of the fluid is q j , cp j is the specific heat capacity, D j is the heat diffusivity ½D j ¼ k j =ðq j cp j Þ, B j is the thermal power deposition density (W/m 3 ) and t is the time variable. The subscript letter j refers to one of the three domains. Second, the gas equation of state determines the relation between the temperature and the pressure in the gas core. Considering low Mach numbers (Ma < 0.01) the pressure can be considered homogeneous in the gas core. Here we chose the ideal gas law as the state function for the gas core in the considered temperature range (between 20 C and 100 C) and pressures (around ambient pressure 10 5 Pa) thus, where P g is the gas pressure, R i is the bubble radius and in our case also the inner oil radius, K is the universal gas constant, T g is the gas temperature and n is the number of moles of the gas that is assumed constant. We neglect vaporization and molecular diffusion phenomena in this derivation. Finally, writing the momentum equation in both the oil and the water phase will determine the dynamic behavior of the system by relating the motion of the fluids to the inner gas pressure. Because of the low Mach number in the water and the oil (Ma < 5 Â 10 À3 ), the radial momentum equation can be derived by integration of the Navier-Stokes equation in potential flow in each domain. The boundary conditions at the interfaces are obtained by balancing the normal stress tensors on each side of the interface. For the gas/oil interface that provides For the oil/water interface, we obtain Here i refers to the gas/oil interface and e refers to the oil/water interface. R is the radius of the interface (gas/oil interface and oil/water) and the superscripts þ and À refer to the outer and inner side of the interface, respectively. l w is the dynamic viscosity of water and its temperature dependence can be written as 18 l w ¼ 2:414 Â 10 À5 Â 10 247:8=ðTÀ140Þ for the water with T the temperature of the water at the water-oil interface. r wo is the oil/water interfacial tension and r o is the gas/oil interfacial tension. The subscripts o, w, and g refer to the oil, the water and the gas, respectively.
The boundary conditions are inserted into the integrated momentum equation, leading to where P 1 is the pressure far away from the bubble. In principle, this derivation is identical to that of the Rayleigh-Plesset equation, 19 but now including multiple domains. The above equations can be discretized and used in a finite difference model (FDM).

B. Governing differential equation
A number of simplifications can be made to the description of the physical system presented above to derive an analytical solution. First, the microbubble resonance frequency is assumed to lie in the same range as its acoustical resonance frequency [f r (MHz) R 0 (lm) $ 3.3 (Ref. 20) verified a posteriori]. During fast thermal processes (in the transient regime), the thermal boundary layer in the gas will obey d ¼ ffiffiffiffiffiffiffiffiffi ffi pD g t p where D g is the thermal diffusivity of the gas and t is the time. The temperature in the gas can then be considered constant when the establishment of the thermal boundary layer is faster than the variations in deposited heat by the modulated laser beam. This holds for bubbles smaller than R ¼ ffiffiffiffiffiffiffiffi ffi pD g 6:6 r % 11 lm: This limit comprises the range of bubble sizes relevant for medical use. Thus, we consider both the pressure and temperature to be homogeneous in the gas core. Therefore the bubble is considered to oscillate around its equilibrium state given by the solution of the static heat diffusion equation. One then easily obtains R i;eq ¼ R i;0 T g;eq P 0 T 0 P 1 þ 2r go R i;eq þ 2r ow R e;eq 2 6 4 3 7 5

1=3
: Here k is the thermal conductivity, the subscript eq refers to the equilibrium state and B a,av (in W/m 3 ) is the average thermal power deposited by the laser. Over the course of each oscillation, the temperature has to change not only in the oil layer but also in a layer of water corresponding to the thermal diffusion radius over half of the excitation period: with f las the laser driving frequency. Here, we are primarily interested in frequencies in the MHz range, corresponding to commonly used ultrasound imaging frequencies. Higher frequencies offer higher resolution whereas lower frequencies offer deeper penetration. To a first approximation, the thermal diffusion radius is then equal to r d % ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pD w =2f las p $ 0:1 lm with a corresponding volume V w0.1 . The heat capacity of the gas, i.e., the thermal energy that can be stored in the gas is negligible as compared to the specific heat of the oil and the water and can therefore be omitted. The heat diffusion in the gas is rapid when the condition given by Eq. (6) is satisfied. Thus, the temperature variation in the gas can be estimated by considering the change in enthalpy of the system following a change in the temperature, with V oil the volume of the oil layer. This temperature variation can then be inserted into the equation of state of the gas together with the equation for the equilibrium radius [Eq. (8)] to obtain the gas pressure as a function of the system parameters and as a function of the initial and equilibrium bubble radii. In turn, the gas pressure can be replaced in the modified Rayleigh-Plesset equation [Eq. (5)] to give

C. Linearization and resonance behavior
We can now differentiate Eq. (10) to suppress the integral and approximate R i to R i,eq in the non-linear products giving where If the purely non-linear term 2c € R i _ R i is neglected, Eq. (11) then becomes Here x is the angular frequency and the underline refers to the complex notation. A consequence of the term 1/jx (that has a Àp/2 phase) in Eq. (12) is a Àp phase at resonance instead of Àp/2 for an acoustically driven bubble. This is a consequence of the necessary integration of the heat deposition as the energy is the quantity driving the system. One notices that the direct influence of the laser intensity appears only in the gain of this transfer function. Equation (12) leads to the undamped natural frequency of the system, where x ¼ R i;eq =R e;eq is a non-dimensional variable describing the influence of the oil layer thickness. x also includes a secondary influence of the laser induced heating via the dilatation of the bubble denoted by the subscript eq: R i;eq > R i;0 implies R i;eq R e;eq < R i;0 R e;0 : The undamped natural frequency of the laser-driven microbubble [Eq. (13)] has a form similar to that of an acoustically driven bubble and is (mostly) inversely proportional to the equilibrium bubble radius. This corresponds to the similar nature of the volumetric oscillations in both cases. For comparison, if the bubble is acoustically driven, R i,eq ¼ R i,0 , R e,eq ¼ R e,0 , P g,eq ¼ P g,0 and the term f becomes f ¼ 3jðP g;0 =R i;0 Þ where j is the polytropic exponent of the gas. j is then also the corrective factor of the description of the thermal behavior of the gas between the acoustically and thermally driven bubbles. Neglecting the surface tension terms for the larger bubbles, the resonance frequency for both driving modes will differ by the quantity ffiffiffi j p . For the smaller bubbles where the surface tension dominates, both expressions become identical. As a remark, both resonance frequencies for an acoustically driven and a thermally driven bubble become identical when j ¼ 1, i.e., in the isothermal case.
This equilibrium radius term R i,eq can in fact be generalized to be the average bubble radius at any time following the reasoning presented above. The expression is therefore also applicable to the transient thermal state.
The calculated resonance frequency given in Eq. (13) also has a strong dependency on the density of the oil and interfacial tensions as one would expect from the mechanical nature of the system as was also shown by Church. 21 From the same equation one can obtain the theoretical damping coefficient of the microbubble, that is, Thus, the damping coefficient depends on the oil and water viscosities, oil and water densities, as well as the bubble size.

D. Pulsed laser excitation
The previous sections describe the case of a modulated continuous wave laser exposure. In practice, pulsed lasers are the preferred excitation tools in photoacoustics. It is therefore worth investigating the case of pulsed light excitation. The derivation of the differential equation for the motion of the microbubble in this case also follows from Eq. (5). However, the different timescales are now well separated. First, the heating of the oil can be treated as instantaneous as the duration of the pulse is typically a few nanoseconds. Then, the gas heats up by diffusion over a timescale given by s g ¼ R 0 2 =D g where R 0 is the bubble radius and D g the gas thermal diffusivity. s g is typically a few hundreds of nanoseconds, which, following the argument of Eq. (6) is much faster than the bubble motion timescale. Finally, the heat diffuses away from the bubble over a timescale s w ¼ R 0 2 =D w where D w is the thermal diffusivity of the water. R 0 is also the terminal thickness of the thermal boundary layer in the static case. s w here typically reaches tens of microseconds. An estimate of the impact of heat diffusion on a short timescale can be obtained from energy conservation, such that where the small amplitude approximation allows the bubble surface area to be taken as the resting surface area S ¼ 4pR 2 i0 , where the approximate expression for the development of the thermal boundary layer for short times, ffiffiffiffiffiffiffiffiffiffi pD w t p , is used and where the influence of the thickness of the oil layer on the global heat diffusion is neglected. The temperature in the gas then becomes where H represents the Heaviside step function and T hot the temperature of the oil (and gas by extension) just after the laser pulse. T hot can further be approximated to T hot % F a / q o cp o þ T 0 where F a is the thermal energy deposited by the laser per unit volume. Subsequently, the first term of Eq. (10) simply becomes Although simpler, this term is actually similar to that found for the continuous wave (CW) laser and the resulting equation shows the same characteristics in terms of resonance frequency and damping, without the need for the additional differentiation performed prior to obtaining Eq. (11).

E. Amplitude response and parameter space
The parameter space of the problem has four dimensions: the thermal power deposition density, which is proportional to the laser intensity and the oil absorption coefficient, the oil layer thickness, the laser modulation frequency and the initial bubble radius. The amplitude of the microbubble response at resonance can be derived from Eq. (12) for Similar to x, the variable y ¼ R i;eq =R i;0 is a measure of the thermal dilatation of the microbubble. e w ¼ 0.1 lm is the thickness of the water layer subject to thermal diffusion during the irradiation cycle. Basically, Eq. (19) is a normalization function that must be calculated from experimentally measured quantities. The presence of the variables x and y prohibits a simple scaling as far as the laser intensity and oil layer thickness are concerned. On the other hand, varying the frequency alone and using Eq. (13) enables us to reduce the resonant oscillation amplitude description in Eq. (18) to a x À5=2 scaling law for a thin oil layer and a x À3 scaling law for a thick oil layer. Thus, a thermally driven microbubble experiences a decrease in the maximum response according to a À5/2 power law with increasing excitation frequency for a thin oil layer. As a consequence the parameter space can be divided in two subspaces. The first is the modulation frequency and the second captures the influence of the oil layer thickness and the laser intensity described using Eq. (19).

III. FINITE DIFFERENCE MODEL (FDM)
In Sec. II, we derived a simple analytical expression that describes, after a number of approximations, the thermal and mechanical behavior of the microbubble. The validity of these approximations, using the same initial set of equations, can be verified by means of a more refined simulation. We therefore design a finite difference simulation in which the three physical equation describing the problem, i.e., the heat diffusion equation, the momentum conservation in the water, and the thermodynamic behavior of the gas core. In short, each equation is simulated as such and communicate with the others at each time step. The mass of each grid volume is defined to be constant and the grid is recalculated every time step to hold on to this definition. Initially, before the laser is turned on and before the bubble starts to oscillate, the grid is defined with a regular size interval of 31.25 nm up to a radius twice the typical bubble size. Beyond this radius the grid gets 15% bigger for each step outward. Without an excessive amount of grid points, the outer point, kept at room temperature, is at a distance well over 3000 times the typical bubble radius making the thermal drift due to a finite simulation volume negligible. The gas is considered to obey the ideal gas law. The pressure is thus defined by where P k V k =T k can be written as where P is the pressure, V is the volume. The subscript 0 stands for the initial value before the laser is turned on and k is the index on the grid spacing vector p and the radius vector r. The speed of the bubble wall _ R is assumed to be much smaller than the speed of sound in air or water so that the pressure in the bubble is considered homogeneous. Together with mass conservation, this leads to a condition on the radius of the grid points in the bubble, where r k is the radius in meters. To speed up the calculation, a Taylor expansion is carried out to the order 3. The oil and the water are assumed to be incompressible, and therefore, the volume in each grid volume in the oil and water domains remains constant, The heat convection diffusion equation writes where I is in units of W/m 3 . The heat equation can be approximated with a central difference scheme in space and a forward finite difference scheme in time, which gives a condition on the temporal evolution of the temperature (Fig. 2). The momentum equation defined in Eq. (5) was also discretized in time using finite differences. Both the velocity and the acceleration are calculated at each time step. More details of the derivation, including the discretization steps, are given in Appendix B.

A. ODE integration method
The ordinary differential equation (ODE) integration for both the non-linear and linear bubble dynamics equations were performed in MATLAB (version R2012a, The Mathworks, Natick, MA). For the non-linear simulations, the integration was performed using the ODE 113 solver as the problem is non-stiff for the considered range of parameters. Simulations were performed to cover the full threedimensional parameter space using both a sine and a square wave modulation for the laser.

B. Finite difference model
The scheme as describe in Sec. III was coded as such and used to simulate the responses of the microbubbles upon laser irradiation using Fortran and is hereafter referred to as finite difference model/simulations (FDM).

C. Parameter estimate
The four key parameters relevant for biomedical imaging applications are the thermal power deposited by the laser, the oil layer thickness, the microbubble initial radius, and the laser modulation frequency. The microbubble size ranges from 1 to 10 lm, corresponding to both the medically relevant bubble size range and to the upper limit of validity of the proposed theory, as demonstrated below. In the present study, the oil layer thickness was varied from 0.1 lm to 3 lm. The power deposition is, to a first approximation, equal to I abs with abs the absorption coefficient of the oil. A value of 2500 m À1 was chosen for the absorption that can be reached by dissolution of a dye such as oil red in an organic oil. The laser intensity was varied from 0.25 Â 10 10 W/m 2 , which is necessary for obtaining a significant bubble response, to 2.5 Â 10 10 W/m 2 that can still be easily reached in practice.

A. Undamped natural frequency
The undamped natural frequency, Eq. (13), was calculated for triacetin oil that has been used in previous studies 17 to coat microbubbles; and heptane oil that is a commonly used linear chain organic oil. The results are displayed in Fig. 3. As expected the natural frequencies for the range of bubble sizes considered lie within the MHz range corresponding to that used in standard ultrasound and photoacoustic imaging systems. There is a noticeable difference in natural frequency between the heptane-coated and triacetin-coated microbubbles owing to the different densities of the two oils. In the same figure, we also plot the resonance frequency of a free gas bubble driven by ultrasound. The undamped natural frequency of a laserdriven oil-coated microbubble differs from that of the acoustically driven bubble due to both the different properties of the microbubbles (oil density, interfacial tensions, In order to fully characterize a dynamic system it is essential to determine the damping coefficient. In the present case the damping coefficient is given by Eq. (14) and depends primarily on the oil and water viscosities, as well as the thickness of the oil layer. Common oils have a specific density of approximately 0.7, but the viscosity can vary over two orders of magnitude, from as low as 386 lPa s for heptane to as high as 17 mPa s for triacetin for example. Figure 4 shows the variation in the damping coefficient with radius for microbubbles coated with heptane and triacetin. The difference in viscosity between these two oils leads to a difference in damping of more than an order of magnitude. A damping coefficient as high as 0.5 for a 3 lm triacetin-coated microbubble drastically limits the benefit of the mechanical resonance of these microbubble as compared to the heptane-coated bubbles. In practice, there will be additional dissipation due to the necessary use of stabilizing agents thereby decreasing even further the oscillation amplitude. Thus, using a low viscosity oil is a necessary condition. Since the damping coefficient is also a function of the oil layer thickness, there is a further difference between high viscosity oils (>1 mPa s) and low viscosity oils (<1 mPa s): for heptane, the damping decreases when increasing the oil thickness whereas for triacetin, the damping increases with increasing oil thickness. This is captured by Eq. (14) and shown in Fig. 4.

C. Low viscosity oil in the parameter space
We have shown that using a low viscosity oil is crucial for obtaining a strong bubble response. For the remainder of the study, we will therefore consider heptane-coated bubbles only. The shape of the excitation waveform when using ultrasound is limited to quasi-sine waves due to the limitations of the relatively narrowband transducer technology. The laser intensity on the other hand can be modulated using any arbitrary waveform with frequency components up to hundreds of MHz, which is the maximum frequency of the currently available acousto-optic modulators. In contrast to ultrasound, light penetration in not further limited by the modulation frequency. Preliminary simulations were therefore performed using sine and square wave modulated waveforms. The latter showed a slightly higher efficiency as the oscillation amplitude is slightly larger for the same energy deposition.
The simulated response of heptane-coated microbubbles to a square wave modulated laser as a function of the initial bubble radius for different laser intensities, oil layer thickness, and modulation frequency is shown in Fig. 5. show the corresponding phase difference between the laser excitation and the microbubble response. The response simulated with both models is very similar both in amplitude and in quality factor, showing that the proposed theory is representative of the physical problem simulated in the FDM. The only significant difference between the finite difference and the non-linear model is the resonant radius of the initial bubble. This discrepancy originates from the fact that an equilibrium assumption is made in the non-linear theory whereas the FDM is applied over a period of 100 ls. The fully  developed thermal regime in the FDM is however only reached after several milliseconds of exposure due to the relatively slow thermal processes. A laser exposure duration of 100 ls places the system in a regime where the thermal processes (with the exception of the high-frequency modulation) evolve on a much shorter timescale as compared to the bubble oscillations. In this quasi-stable regime, the gas core temperature can lag by 20% to 30% from the equilibrium temperature as demonstrated in Appendix C. Figures 5(a) and 5(d) depict the variation of the microbubble response when increasing the oil layer thickness. The oscillation amplitude decreases for a thicker oil layer despite the decreasing damping coefficient. It can also be seen from Figs. 5(a) and 5(d) that the response amplitude passes through a maximum for a heptane oil layer thickness of 1.1 lm. Physically, for a thin oil layer, the energy deposited in the oil is quickly transferred to the gas and the surrounding water and thus the temperature in the system changes quickly enough to follow the excitation. In this regime, the microbubble response is limited by the energy deposited by the laser that is to first order proportional to the volume of absorbing oil. Above a critical oil layer thickness, the time required to change the temperature in the oil is no longer negligible compared to the laser modulation period. Thus, the temperature in the oil fails to follow the variation in the heat deposited by the laser. This phenomenon of thermal inertia decreases the amplitude of the temperature variation in the gas, and therefore decreases the response amplitude. This effect can be seen mathematically in the term a and its dependency on the oil volume and transiently heated water volume. Asymptotically, the maximum of the resonance curve can be written as in the case of a thin oil layer and for a thick oil layer. Figures 5(b) and 5(e) show the variations in the microbubble response as a function of the heat deposited by the laser. As expected, the microbubble oscillations become stronger with increasing laser intensity. Interestingly, and unlike the temperature variation, the response does not increase linearly with the intensity. The laser intensity in Eq. (19) is represented by the variable y that is a measure of the thermal dilatation of the microbubble. Nevertheless, in Eq. (19), the laser intensity is presented as a proportional term that describes the influence of the density of heat deposition and an inverse term that corresponds to the change in bubble size with increasing temperature. Figures 5(c) and 5(f) present the simulation results for different modulation frequencies. The response of the laserdriven microbubbles at resonance increases strongly when decreasing the laser modulation frequency. Physically, the temperature variations in the gas core and therefore the driving pressure will vary with the energy deposited during half the period of the laser excitation and will therefore be larger when the modulation frequency becomes lower. Finally, Fig. 5(g)-5(i) show the variation of the phase difference between the microbubble oscillations and the laser excitation. As expected from the linearized equations, the phase varies from Àp/2 to À3p/2 and crosses Àp at the natural frequency. Thus, the bubble oscillates in anti phase with the laser excitation.

D. Scaled resonance curves
In order to apply these theoretical findings to an experimental case, one must consider the practical difficulties involved in producing stable microbubbles with the same oil thickness and to expose them to the same laser intensity. It is therefore desirable to normalize the microbubble resonance curves using parameters that are experimentally accessible. This can be achieved using the normalization functions given in Sec. II E and Eq. (19), derived from the linearized equation Eq. (12). We also know from Eq. (13) that the resonant bubble size can be described as a function of the "hot" bubble radius. This hot bubble radius can thus be used instead of the initial bubble radius and then be normalized to the resonant radius using Eq. (13).
The resonance curves simulated from the non-linear theory for varying intensity and oil thickness and normalized using Eq. (19) are plotted in Fig. 6(a) and those for a varying laser modulation frequency normalized using the power laws derived in Sec. II E are plotted in Fig. 6(b). The normalization function from Eq. (19) is effective for scaling the microbubble responses simulated using the non-linear theory. A scaling according to a À5/2 power law gives a similar result for varying modulation frequencies but with a larger deviation for the lowest frequency (500 kHz). The resonance curves simulated by the FDM and scaled with the same equations are plotted in Fig. 7. For the FDM, the amplitude of the scaled curves for the oil layer thickness and laser intensity [ Fig. 7(a) and 7(b)] present a larger amplitude and deviation, which is mostly due to the differences between the thermal equilibrium radius used in the theory and the quasi equilibrium average radius simulated by the FDM. The proposed normalization applied to the results of the FDM reduces a sevenfold variation in the response amplitude to a small error margin. Figure 7(c) shows the FDM simulation for different excitation frequencies normalized by the expected À5/2 power law.

E. Harmonics and subharmonics
Beyond the strength of the fundamental resonance, one feature of microbubble oscillations has become increasingly important and been widely investigated for ultrasound imaging: the harmonic and subharmonic pressure wave generation. Figures 8(a) and 8(b) show the harmonic and subharmonic microbubble oscillations, respectively, relative to the fundamental response on a dB scale for square wave excitation and for different heat deposition densities. Figures 8(c) and 8(d) depict the same quantities for a sine wave laser modulation. Both the square and sine wave laser modulation generate harmonics from bubbles around the resonant size but the square wave also generates significant harmonics from much smaller bubbles. We attribute this to the harmonic composition of the square wave itself that includes a higher frequency component at three times the fundamental frequency. The generation of harmonics thus depends strongly upon the choice of the laser modulation waveform.

VI. DISCUSSION
As discussed above, the FDM simulations were run over a period of 100 ls, allowing for the model to reach a quasi steady-state that is nonetheless significantly different from the perfect equilibrium state used as a reference in the theory. Simulations were also run over a longer timescale with non-modulated laser driving and converged toward the expected thermal equilibrium. A duration of 100 ls or less is, however, more relevant for practical application of these bubbles, as enough time should be allowed for imaging and signal integration/accumulation whilst avoiding excessive heat deposition in the tissue. The theoretical model could be modified to match this quasi steady state but the timescale to choose would then depend on the experiment or application.
The observed increase in signal amplitude could potentially increase the tissue depth from which photoacoustic images can be obtained in two ways. First, the use of a sensitive contrast agent would increase the signal to noise ratio. Second, the use of a CW laser offers temporal integration possibility, together with an optimal use of the bubble resonance effect. This last approach is very novel and now needs to be translated into (pre)clinical practice.
The normalization function in Eq. (19) derived from the linear theory fails to satisfactorily scale the thinner oil layers (%100 nm) and these are therefore not presented in Fig. 6(a) and Fig. 7(a). We justify this omission by the much lower predicted response for such thin oil layers [ Fig. 5(a) and 5(c)] that therefore present a more limited interest in terms of normalization.
Both the proposed theory and the FDM were considered in the context of incompressible potential flow. In fact, an oscillating bubble emits an acoustic wave, therefore producing a weakly compressible flow carrying energy away from the microbubble. This effect has been addressed in earlier studies, such as those by Keller et al. 22 and Prosperetti et al. 23 From these models, an approximate form of the losses by radiation can be included as a supplementary pressure term in Eq. (10) for low Mach numbers, where P g is the gas pressure. Using the state function of the gas and Eq. (9), Eq. (22) becomes Here the dominant damping term is on the right side. The addition of this term in the theory and corresponding simulations produced only a negligible effect on the total damping. Reradiation is therefore negligible compared to viscous dissipation, even in the low viscosity case of heptane oil-coated bubbles.

VII. CONCLUSIONS
On the basis of theoretical considerations, we have demonstrated the possibility of using the mechanical resonance of microbubbles in the MHz frequency range, typically used in ultrasound imaging, to increase the signal amplitude in photoacoustic imaging. We have developed a theory supported by finite difference model simulations that clarifies the underlying phenomena and predicts the natural frequencies and resonance characteristics of laser driven microbubbles. The resonance frequency of the laser-driven microbubbles is different from the acoustically driven bubble by a factor ffiffiffi j p , j being the polytropic exponent of the gas. The natural frequencies of these microbubbles depend on the density of the oil coating, as was already shown for acoustically driven bubbles. The proposed theory also indicates the crucial importance of choosing an oil with a low viscosity, i.e., similar to that of water in order to obtain a sufficiently high quality factor for the resonance. It appears that a suitable selection of the oil properties can improve the overall quality of the system, even above that of acoustically driven bubbles. We have also extracted from the linearized theory normalization functions within the investigated parameter space that will allow for the scaling of datasets using experimentally accessible parameters to reduce significantly the error due to variations in oil thickness and laser intensity. The proposed theoretical considerations are valid for continuous wave exposure and can also be applied to the more typical case of pulsed excitation photoacoustics. Similarly, the response of optically absorbing microbubbles that do not contain an oil layer can be extrapolated from the proposed work by appropriate selection of the equivalent parameter set (density, thickness and viscosity) for the absorbing layer. Finally, we have shown that laser-driven microbubbles are expected to exhibit similar behavior to that of acoustically driven microbubbles in terms of harmonic and subharmonic generation, thereby opening new possibilities for harmonic photoacoustic imaging that could potentially enable a much improved contrast to tissue ratio.

ACKNOWLEDGMENTS
This project was made possible by the funding of the Dutch national NanoNextNL program, a micro and nanotechnology consortium of the Government of the Netherlands and 130 partners and the UK Engineering and Physical Sciences Research Council (Grant No. EP/ I021795/1). We also warmly thank Rodolfo Ostilla Monico for his valuable advices concerning the finite difference simulation.

APPENDIX A: MATHEMATICAL DERIVATION OF THE THEORETICAL MODEL
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 q @v @t þ v @v @r ¼ À @P @r þ lr 2 v: In the simulation, the bubbles will have an oscillation amplitude of the order of lm and the frequency will be in the order of MHz. Speeds will therefore be approximately 1 m/s and thus much lower than the speed of sound. For this reason, incompressibility of the liquid is assumed leading to where _ R is dR=dt. 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 lr 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, ; (A2) rewritten as where r o is the strain tensor andẽ r denotes that it is in the r direction. dP 1 is the difference in pressure over the oil-gas interface.
where u 0 is the velocity derivative to the radius. r is the surface tension.
Over the oil water interface r w Áẽ r À r o Áẽ r ¼ dP 2 ;

Rewriting the equation above then gives
Resulting in

Combining
We know that P g À P 1 ¼ PðR À i Þ À P 1 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 Þ. We can rewrite by adding and subtracting similar terms, Part 1 of (A6) is defined in (A4), part 2 is defined in (A3), part 3 is defined in (A5) and part 4 is defined in (A1). Thus, the complete equation is Rewriting gives and rewriting further To reach the modified Rayleigh-Plesset equation, where the viscosity of water l w is temperature dependent following the following relation: 18 l w ¼ 2:414 Â 10 À5 Â 10 247:8=ðTÀ140Þ ; where T is the temperature of the water at the water-oil interface. This new RP equation reduces to the classic RP equation for a bubble with only one liquid around it when the properties of water and oil are chosen to be identical.

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 R ¼ pD g 6:6 % 11lm: A diffusion distance estimation can be made for the water to find app. 0.1 lm for 1 MHz frequency.
The change in temperature over time can be described as follows: with q o the density of the oil, V w0.1 the volume of the first 0.1 lm of water, cp o the heat capacity at constant pressure of the oil and B a (t) the absorbed laser power (W/m 3 ): Using the initial equilibrium solution, B a,av is the average laser power deposited (W/m 3 ). Filling in C 1o and C 2o and rewriting gives We know that where T g can be substituted 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 1 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 g as the total pressure in the gas, and filling in T g and the found pressure, this gives Organizing for _ R; _ R 2 , and € R and, as an approximation, taking all R i and R e to be R i,eq and R e,eq in case they are multi- In order to find an equation that does not contain an integral, everything is derived to time: In which case _ R e is assumed to be approximately _ R i and T gas,eq þ T room is now called T g2 . The term 2c € R i _ R 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 x 0 the angular eigen frequency. 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 p at resonance instead of p/2 such as in the normal RP equation: From Eq. (A15), we can find T g2 R 3 i;0 R 4 i;eq ¼ T 0 P g;eq P 0 R i;eq : (A29) Therefore, f can be simplified where the equilibrium pressure P g,eq is the atmospheric pressure plus the Laplace pressure jump over both interfaces, i;eq þ 2r ow R i;eq R e;eq ! ; (A31) (A32) x 0 is not a function of time so all R i and R e are now R i,eq and R e,eq : 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), The simulation is based on a time dependent mesh. The mass of each grid volume is defined to be constant and the grid is recalculated every time step to hold on to this definition. Initially, before the laser is turned on and before the bubble starts to oscillate, the grid is defined with a regular size interval of 31.25 nm up to a radius of 6 lm which is twice the typical bubble radius. Beyond this radius the grid gets 15% bigger each step outward. This way, without having an excessive amount of grid points, the most outer grid point, kept at room temperature is more than 3000 times the typical bubble radius making negligible the thermal drift due to a finite simulation volume.

The ideal gas law
In order to find an expression for the pressure, we look at the ideal gas law, where P is pressure, V is volume, m is mass,l is the molar mass in kilograms per mole, R g is the ideal gas constant. The simulation is defined such that the mass in each grid-volume remains constant in time: The speed of the bubble wall _ R is much smaller than the speed of sound in air or water. Therefore the pressure in the bubble is considered homogeneous. With this information it can be shown that the pressure is defined by where P V k =T k can be written as with R e being the radius of the bubble at the oil water interface. Filling in for what was found in Appendix B 5, Rearranging gives with k being the gridpoint on the boundary between water and oil. Similarly, the boundary condition between the oil and the gas is with R i being the radius of the bubble at the gas-oil interface. Rearranging this gives