## 1.

## Introduction

The cells of the retinal pigment epithelium naturally contain highly absorbing, micrometer-size melanosomes. Other biological and physical systems can be doped with highly absorbing micrometer- and nanometer-size absorbers, such as gold particles.^{1, 2} Their strong absorption relative to normal tissue makes these particles the likely site for the generation of light induced damage at threshold levels of illumination. In addition to the potential danger, the same biophysical mechanisms can be used for a variety of beneficial medical applications.

Much research has been carried out, both experimentally and theoretically, investigating the interaction of laser light with the visual system^{3, 4, 5, 6} and other biological tissues.^{7} Absorption of the laser energy generates various responses that include temperature rise^{8, 9, 10} (thermal), pressure^{11, 12, 13, 14, 15, 16, 17, 18} (acoustomechanical), and vaporization^{17, 19, 20} (phase transitions). In this paper, we present results for the modeling of the generation of dangerously high pressures and explosive vaporization utilizing a theoretical treatment that incorporates a full empirical equation of state of water. This is the most realistic treatment we have seen and enables us to predict thermomechanical effects, including nonlinear shockwaves and phase transitions that produce bubbles. This treatment is designed to predict the strength of shockwaves and the size of bubbles produced by laser pulses over a wide range of pulse duration and fluence that are absorbed by particles with known physical properties. However, the theoretical model can also be used in “reverse.” Experimental measurements of shockwave strength and bubble size can be used with this analysis to determine mechanical and thermal properties of absorbing particles, such as melanosomes, that are too small for conventional measurements of these properties.

Figure 1
is a schematic showing how the threshold laser fluence level
${I}_{0}$
for biological danger to the cell varies as a function of pulse duration.
${\tau}_{L}$
for each of the three damage mechanisms: temperature rise (thermal), explosive vaporization (bubble generation), and shockwaves. All three mechanisms can operate at all pulse durations. However, since their dependencies on
${\tau}_{L}$
are quite different, the mechanism responsible for threshold damage may change as a function of
${\tau}_{L}$
. For long
${\tau}_{L}$
, thermal damage requires the lowest
${I}_{0}$
and is the threshold damage mechanism. As
${\tau}_{L}$
shortens, thermal damage requires less fluence until
${\tau}_{L}$
decreases to a characteristic heat conduction time
${\tau}_{H}$
, that depends on absorber size and the spacing between absorbers, as explained in Ref. 9. For
${\tau}_{L}<{\tau}_{H}$
, the threshold fluence remains constant for thermal damage. The damage threshold fluence for explosive vaporization also decreases as
${\tau}_{L}$
shortens, but continues to decrease for
${\tau}_{L}<{\tau}_{H}$
. The characteristic vaporization bubble time
${\tau}_{V}$
is equal to the time scale of bubble expansion and for the system treated here is of the order of
$100\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
. For pulse durations that are longer than this time, part of the laser energy will be absorbed after bubble growth has effectively ended, i.e., for a given fluence, the maximum bubble radius, and hence the damage threshold levels, will significantly depend on the pulse duration. Below
${\tau}_{V}$
, the fluence required for damage from explosive vaporization remains approximately constant. Since
${\tau}_{V}<{\tau}_{H}$
, for
${\tau}_{L}$
between
${\tau}_{H}$
and
${\tau}_{V}$
, the fluence required to cause damage from explosive vaporization may drop below the fluence needed to cause thermal damage, and explosive bubble generation may become the threshold damage mechanism.^{19} The fluence required for shockwave production decreases as
${\tau}_{L}$
decreases and becomes constant when
${\tau}_{L}<{\tau}_{c}$
, where
${\tau}_{c}$
is a characteristic time for shockwave production. Since
${\tau}_{c}<{\tau}_{V}$
, for
${\tau}_{L}$
between
${\tau}_{V}$
and
${\tau}_{c}$
, the fluence required to cause damage from shockwaves may drop below the fluence needed to cause bubble damage, and shockwave production becomes the threshold damage mechanism. Thus, the underlying biophysical mechanism for damage at threshold fluences may change as
${\tau}_{L}$
enters the regions between the characterisitic times
${\tau}_{H}$
,
${\tau}_{V}$
, and
${\tau}_{c}$
. In addition to different time scales, each of these mechanisms has different characteristic strengths and length scales. These differences can be important both for assessing danger, and in beneficial medical applications. In this paper, we show how to calculate
${\tau}_{V}$
and
${\tau}_{c}$
, as well as determining how the strength and extent of vaporization and shockwave effects depend on
${\tau}_{L}$
and
${I}_{0}$
.

## 2.

## Model

## 2.1.

### Governing Equations for Absorbing Melanosome

Our model^{11, 12} consists of a uniform spherical absorber surrounded by a transparent medium. The rate of energy input per unit mass is given by^{19}

## 1

$${\stackrel{\u0307}{I}}_{e}=\frac{3{I}_{0}}{4a{\tau}_{L}{\rho}_{0}}\{1-\frac{1}{2{\alpha}_{L}^{2}{a}^{2}}[1-\mathrm{exp}(-2{\alpha}_{L}a)(1+2{\alpha}_{L}a)]\}=\frac{3{I}_{0}}{4a{\tau}_{L}{\rho}_{0}}C\left({\alpha}_{L}a\right),$$In these investigations, we wish to predict the thermomechanical effects, such as shockwaves and bubbles, produced by a laser pulse. This involves the calculation of the motion, pressure, density, temperature, and entropy at all locations within the absorber and in the surrounding medium. The complexity of the system prevents the calculation of analytical solutions and requires us to solve the equations numerically. The position at which these parameters are calculated are given in the Eulerian space-fixed coordinate system
$\mathbf{r}$
. Since we are investigating a spherical absorber, we need only be concerned with the radial coordinate,
$r$
. We use a finite difference algorithm to perform the numerical calculations.^{12} We use a variable grid of mass elements,
$i=1,\dots ,N$
. The position of each mass element is denoted by Lagrangian coordinates
${u}_{i}={u}_{i}\left(t\right)$
. The numerical value of each
${u}_{i}$
gives the location of element
$i$
in terms of the coordinate
$r$
. Each element
$i$
is given an initial position
${u}_{i}(t=0)={r}_{i}$
. Each element is also given an initial density and absorption coefficient, which distinguishes whether the element is in the absorber or outside in the medium. The computer code is written so that elements
$i=1,\dots ,n$
are inside the absorber, and
$i=n+1,\dots ,N$
are assigned properties of the aqueous medium.

The mathematical dot operation $\stackrel{\u0307}{f}\left(t\right)$ means a total time derivative, the spatial derivative $\nabla $ is taken with respect to $\mathbf{r}$ , while the spatial derivative with respect to $\mathbf{u}$ at a specific mass point is specifically denoted by ${\nabla}_{\mathbf{u}}$ . With this notation, the equation of motion of a point inside the absorbing sphere is

where $P$ is the pressure and $\rho =\rho \left(t\right)$ is the time-varying density, which is related to the static density by mass conservationwhere $u$ is the radial component of the vector $\mathbf{u}$ . In spherical geometry, the equation of motion can be expressed asIf we assume that the absorbing melanosome has a constant bulk modulus $B$ and constant thermal expansion coefficient $\alpha $ , the equation of state (EOS) of the absorber can be written as

where $v=1\u2215\rho $ is the specific volume and is related to $\mathbf{u}$ by## 6

$$\stackrel{\u0307}{v}=\frac{\partial v}{\partial t}+\stackrel{\mathrm{\u0307}}{\mathbf{u}}\bullet {\nabla}_{\mathbf{u}}v.$$## 7

$$\frac{\stackrel{\u0307}{v}}{v}={\nabla}_{\mathbf{u}}\bullet \stackrel{\mathrm{\u0307}}{\mathbf{u}}=-\frac{\stackrel{\u0307}{P}}{B}+\alpha \stackrel{\u0307}{T}.$$Energy conservation of the absorber relates the rate of absorption to temperature rise, volume change, and heat lost through conduction to the surrounding medium:

## 8

$${\stackrel{\u0307}{I}}_{e}=T\stackrel{\u0307}{s}-\frac{\lambda}{\rho}{\nabla}_{u}^{2}T={c}_{v}\stackrel{\u0307}{T}+B\alpha T\stackrel{\u0307}{v}-\frac{\lambda}{\rho}{\nabla}_{u}^{2}T,$$## 2.2.

### EOS for Water

The equation of motion, Eq. 4, the EOS, Eq. 7, and the conservation of energy, Eq. 8, constitute the governing equations for the absorbing melanosome. Analogous, but not necessarily similar equations must be obtained for the surrounding medium, which is treated as water. Equation 4 is the same for the absorber and the aqueous medium. Equation 8 has one straightforward modification in that we assume that the medium is transparent and therefore, ${\stackrel{\u0307}{I}}_{e}=0$ .

The most important difference is in the EOS for the aqueous medium. Equation 5 for the absorber is valid only for a material that undergoes no phase changes and experiences changes in its specific volume that are small enough to assume that $B$ remains constant. In Ref. 11, we investigated laser pulses that were weak enough to satisfy these conditions for both the absorber and the aqueous medium. This enabled us to linearize the model and obtain analytical expressions for the acoustic waves that were generated in the medium.

In this paper, we are especially interested in shockwaves and phase transitions (bubbles) in the medium. This requires an EOS that covers a broader range of conditions than that of Eqs. 5, 7, 8, which is used in this paper only for the solid absorbing melanosome. We expand on our treatment of Ref. 12 by employing a variable grid lattice to enable more extensive calculations of the dependencies of shockwaves and bubble formation as a function of laser pulse parameters. In this paper, we are not interested in effects arising from the boundary of the medium. This requires that our lattice end boundary be located far enough away from the melanosome so that no pressure signals reach it during the simulation time. The variable grid method employs a fine grid near the melanosome, where we are most interested, and becomes increasingly coarser at greater distances. This reduces the total number of spatial lattice points that are included in the simulation and hence the computational time becomes almost linear rather than quadratic in the simulated time. We use an EOS formulation for water from the NBS/NRC Steam Tables of Ref. 21. This full EOS treatment includes most of the true thermodynamic behavior, such as phase changes that produce bubbles. It does not permit metastable pressures and produces the stable vapor phase, bypassing intermediate metastable stages.

This stable equilibrium EOS allows for no negative pressures, since tensile stress is unstable with respect to a transition to the vapor phase. Therefore, our calculations will never lead to negative pressures in the aqueous medium. However, this is not an important omission. The medium surrounding a melanosome is an impure aqueous solution and can sustain only relatively small metastable tensile pressures^{22, 23, 24} compared to those generated by the shockwave producing laser pulses that we are interested in. As a consistency test we performed additional calculations, reported in Ref. 12, using the Tait EOS for water,^{25, 26} which was developed to model shockwaves in water and does allow negative pressure conditions, though no phase changes. Our additional work showed that the strongest and therefore most damaging pressure gradients, occurred at the leading edge of a shock front, which always had a positive amplitude. In addition, the leading edge of the shock fronts predicted by the full EOS are almost identical to those predicted using the Tait EOS, as described in Ref. 12. Even though nonlinear phenomena such as shockwave generation are notorious for being difficult to calculate because of numerical instabilities, the agreement between the two different approaches enables us to feel confident of the predictions from the full EOS for positive pressure transients such as shockwaves, as well as for bubble formation. Details of our nonlinear numerical algorithm are given in Ref. 12.

The full NBS/NRC Steam Table EOS relates for stable water the pressure $P$ , the specific volume $v$ , the temperature $T$ , and the specific entropy $s$ , such that a numerical value of any two of these thermodynamic parameters enables calculation of values for the other two. The full EOS is valid for both liquid water and vapor, as well as the fluid phase above the critical temperature and pressure of $647\phantom{\rule{0.3em}{0ex}}\mathrm{K}$ and $218\phantom{\rule{0.3em}{0ex}}\mathrm{atm}$ . Our treatment enables us to predict the temperature, pressure, and bubble size created by a laser over a wide range of pulse durations $(0.1\phantom{\rule{0.3em}{0ex}}\mathrm{ps}<{\tau}_{L}<1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{s})$ and fluences $({I}_{0}<10\phantom{\rule{0.3em}{0ex}}\mathrm{J}\u2215{\mathrm{cm}}^{2})$ . Though the model employed can be used with lasers of arbitrary shape, in the presented calculations we used a top-hat laser temporal profile. We do this for simplicity and to reduce the number of degrees of freedom in the system.

## 3.

## Shockwaves and Bubbles

In this section, we plot the results of exhaustive calculations for the size of shockwaves and bubbles expected for a wide range of laser pulse durations and fluences. These pulse durations and fluences are chosen because they are in the range for which pressure effects and vaporization may be important in determining threshold levels of retinal damage.

In all the calculations reported here we used the following parameters for the melanosome:^{12, 27}
$a=1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
,
${\rho}_{0}=1.35\phantom{\rule{0.3em}{0ex}}\mathrm{g}\u2215{\mathrm{cm}}^{3}$
, and
${c}_{v}=2.51\phantom{\rule{0.3em}{0ex}}\mathrm{J}\u2215\mathrm{gK}$
. We use an absorption coefficient of
${\alpha}_{L}=\mathrm{10,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
. The value of the absorption coefficient used represents an upper limit for any wavelength, hence it corresponds to a lower limit on threshold levels. As there is still some uncertainty with respect to the value of the absorption coefficient due to the difficulty in performing measurements on biological tissue,^{28, 29} later in this paper we present a table that enables a researcher to scale the fluence required to get the same shockwaves and bubbles for other values of
${\alpha}_{L}$
. Reliable numbers for the bulk modulus
$B$
and bulk thermal coefficient of expansion
$\alpha $
are also not available. Our results depend on the values used for these parameters, and independent measurements of them would be extremely valuable, as emphasized in Refs. 11, 12. To continue with the calculations, we adopt graphite as a substitute because of a chemical similarity to melanin,^{30} and use graphite’s
$B=39.4\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$
and
$\alpha =2.98\times {10}^{-5}\phantom{\rule{0.3em}{0ex}}{\mathrm{K}}^{-1}$
. Later, we discuss how the results will change for other values of
$B$
and
$\alpha $
. The surrounding medium is aqueous with
${\rho}_{0}=1.0\phantom{\rule{0.3em}{0ex}}\mathrm{g}\u2215{\mathrm{cm}}^{3}$
, and thermal conductivity
$\lambda =5.56\times {10}^{-3}\phantom{\rule{0.3em}{0ex}}\mathrm{J}\u2215\mathrm{cm}\phantom{\rule{0.2em}{0ex}}\mathrm{K}\phantom{\rule{0.2em}{0ex}}\mathrm{s}$
. Since for water we are using a tabulation relating
$P$
,
$v$
,
$T$
, and
$s$
, we do not use values for
${c}_{v}$
,
$B$
, and
$\alpha $
as is necessary when using the simple EOS, Eq. 5, for the absorbing melanosome.

In Fig. 2 and 3 we use a high fluence to present plots that are useful in interpreting the results that are presented in later figures. We use a laser pulse of duration ${\tau}_{L}=1\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$ , ${I}_{0}=1.0\phantom{\rule{0.3em}{0ex}}\mathrm{J}\u2215{\mathrm{cm}}^{2}$ (suprathreshold for both shockwaves and bubble formation to elucidate the effects) with ${\alpha}_{L}=\mathrm{10,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ . The melanosome undergoes a $\Delta T=1595\phantom{\rule{0.3em}{0ex}}\mathrm{K}$ , and the shock front in the near field can be as large as $4\phantom{\rule{0.3em}{0ex}}\mathrm{kbar}$ , corresponding to a shock front speed of ${U}_{s}=2630\phantom{\rule{0.3em}{0ex}}\mathrm{m}\u2215\mathrm{s}$ ( $\approx 75\%$ above acoustic speed in water). Once the shock front leaves the near-field region, it decays to an acoustic type disturbance. Though not as dangerous as a shock front, a strong acoustic disturbance may cause biological damage, and this possibility remains to be investigated.

## 3.1.

### Shock Front Generation and Propagation

Figure 2(a) contains multiple snapshots of the pressure profiles outside the melanosome in the aqueous medium at various times after the single laser pulse. Note that in the medium, there are multiple trailing shock fronts due to ringing of the melanosome, as explained in Refs. 11, 12. As a shock front expands outward into the medium, the amplitude of the shock front decreases because its energy spreads out over a larger surface area. Also, the sharpness of the pressure gradient converts mechanical energy into heat. (The numerical algorithm used^{12} in this work incorporates a viscosity-like term that depends on the spatial second derivative of the pressure.^{31} It therefore simulates viscosity in that it is important for shock fronts but is of negligible importance for acoustic waves.)

In Fig. 2(b) we plot the decreasing strength of the leading shock front as it propagates out into the medium, away from the melanosome. The shock front decays faster than the $1\u2215r$ dependence of an acoustic wave in three dimensions. The stronger the amplitude of the shock front, the faster it decays with $r$ . This can be represented by

We found that we could nicely fit the pressure decay curve of Fig. 2(b) with the following dependence for $\beta =\beta \left(P\right)$ with $P$ given in bar. In the log-log plot of Fig. 2(c), the slope of the calculated curve at any point is the decay exponent $\beta $ and we give a few representative values.## 3.2.

### Absorber’s Characteristic Acoustic Time ${\tau}_{c}$

Figure 2(a) shows that there are multiple smaller shockfronts that trail the strong leading shockfront. These trailing shockfronts are due to ringing of the absorber.^{11, 12} The difference in
$B$
,
$\alpha $
, and
$\rho $
of the absorber and the thermomechanical properties of the aqueous medium creates a mechanical impedance mismatch at the surface of the absorber. Part of the mechanical energy leaving the absorber is reflected back toward the center. After reaching the center, the energy rebounds and generates another outgoing pressure wave, which is again partially reflected at the surface. This ringing causes a train of decreasing pressure fronts to emanate out into the medium. The time between pressure fronts in the medium is the time for a pressure wave to travel from the surface of the absorber to its center, and then back to the surface, and is due to the thermo-mechanical properties of the absorber. This characteristic absorber time can be expressed as
${\tau}_{c}=2a\u2215c$
, where
$c=\sqrt{B\u2215\rho}$
is the speed of sound in the absorber. This results in

If $B$ and $\rho $ of the absorber are known, then ${\tau}_{c}$ can be predicted from Eq. 11. This logic can also be reversed. If there are microparticles or nanoparticles whose properties are not known, measurement of the time between pressure peaks in the surrounding medium can be used with Eq. 11 to determine $B$ . Later, we show that ${\tau}_{c}$ not only defines the time between pressure peaks, but also is the maximum laser pulse duration that will generate shock fronts in the surrounding medium.

## 3.3.

### Explosive Vaporization: Bubbles

In Fig. 3 we plot a characteristic graph that reflects bubble formation. For various times, we plot the specific volume of the medium as a function of distance from the surface of the melanosome, which is located at
$r=a$
. The region close to the melanosome with the high specific volume represents a shell of vaporization surrounding the melanosome, which is termed a bubble. Note that this vaporized region is initiated^{32} by heating the liquid water to above
$100\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$
. The region of the aqueous medium farther away, that has a specific volume of approximately
$1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{3}\u2215\mathrm{g}$
, is in the liquid state. When we refer to the size of the bubble produced by a laser pulse, we mean the maximum radius of the vaporized shell measured from the center of the melanosome. The thickness of the vapor shell itself is
$1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
less than this bubble radius. We also note that over the initial
$100\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
the vaporized region grows with a speed that is a factor of 50 times smaller than the speed of the pressure waves.^{19} The figure also shows a peak in the specific volume within the bubble. This is due to a pressure gradient inside the bubble that is maintained by heat transfer from the high-temperature melanosome. If there was not continuous heat transfer from the melanosome the pressure would be approximately constant within the bubble (homobaric assumption) and the specific volume would be a radially decreasing function.

## 4.

## Dependence of Shockwaves and Bubbles on Pulse Duration and Fluence

In Fig. 4 we plot the pressure and bubble results for different pulse durations and fluences. Figure 4(a) is a plot of the strength of the shock front at a location $1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ outside the melanosome $(r=2a)$ . This position is chosen because the shockfront has traveled enough distance from the surface of the absorber to fully form, but a small enough distance to still have significant strength. We can see that the shock front amplitude as a function of ${I}_{0}$ is almost identical for ${\tau}_{L}=1\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$ and $100\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$ , but drops dramatically for ${\tau}_{L}=370\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$ . This ${\tau}_{L}$ is also ${\tau}_{c}$ for the absorber, as already explained, and defines the stress confinement time. For ${\tau}_{L}<{\tau}_{c}$ , the pressure effects in the surrounding medium depend on ${I}_{0}$ , but are independent of ${\tau}_{L}$ . (As explained in Ref. 11, the concept of stress confinement is not valid in the core region of the absorber). For ${\tau}_{L}>{\tau}_{c}$ , the amplitude of the leading pressure front decreases inversely with ${\tau}_{L}$ . In addition to the decreasing amplitude, the pressure front smooths to a sinusoidal acoustic wave, with much smaller pressure gradients than the shockfronts generated by pulses with ${\tau}_{L}<{\tau}_{c}$ . Therefore, pulses with ${\tau}_{L}<{\tau}_{c}$ will cause much greater pressure induced effects in the surrounding medium than pulses with ${\tau}_{L}>{\tau}_{c}$ .

Figure 4(b) shows the maximum bubble size
${R}_{\mathrm{max}}$
as a function of
${I}_{0}$
for different
${\tau}_{L}$
. The bubble size decreases by only 10% as the pulse duration increases from
$1\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
to
$1\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
, which is a manifestation of “cavitation confinement.”^{12} For
${\tau}_{L}=100\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
, there is a significant decrease in the bubble size, and this defines the characteristic cavitation or vaporization confinement time
${\tau}_{v}$
, which is more than two orders of magnitude larger than
${\tau}_{c}$
. Vaporization resulting from laser absorption is mainly due to heat flow into the surrounding medium and is hot cavitation, though cold cavitation bubbles that are generated by a sudden decrease in pressure are incorporated in the EOS that we employ for the medium. Therefore,
${\tau}_{v}$
is determined by heat flow. In the next section, we show that the amplitude of both the shock fronts and bubbles depend on both the thermal and the mechanical properties of the absorber.

In Fig. 5
, we compare
${R}_{\mathrm{max}}$
as a function of
${I}_{0}$
calculated from our numerical model with experimental measurements, as well as predictions from two other models. The curve labeled
$1\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
employs the model of this paper, and is the same
$1\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
curve plotted in Fig. 4(b). The curve labeled
$10\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
is of the same model but uses a pulse duration of
$10\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
. The curve labeled NB is the maximum bubble radius as a function of
${I}_{0}$
obtained by fitting the experimental data of Neumann and Brinkmann^{33} for a pulse duration of
$12\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
. We found that the curve
$[{R}_{\mathrm{max}}\left({I}_{0}\right)\u22151\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}]={[0.3+44({I}_{0}-0.12)]}^{0.5}$
, for
${I}_{0}>0.12$
in
$\mathrm{J}\u2215{\mathrm{cm}}^{2}$
fits best the experimental data, and we use this function to extrapolate the experimental results of their Fig. 4 to the fluence scale of our Fig. 5. The curve labeled Gerstman is calculated from a model^{19} that makes simple assumptions about bubble growth but has the advantage of giving a straightforward analytical expression for
${R}_{\mathrm{max}}$
as a function of
${I}_{0}$

## 12

$$\frac{{R}_{\mathrm{max}}}{a}={\{1+\frac{1}{q{\rho}_{c}}[\frac{3C\left({\alpha}_{L}a\right){I}_{0}}{4a}-{c}_{v}{\rho}_{0}\Delta {T}_{c}]{\left(\frac{{P}_{0}}{{P}_{\mathrm{min}}}\right)}^{1\u2215\gamma}\}}^{1\u22153}.$$We use the following conditions for the parameters in Eq. 12, their choice is more fully explained in Ref. 19. The parameter $q=2770\phantom{\rule{0.3em}{0ex}}\mathrm{J}\u2215\mathrm{g}$ is the heat needed to raise one gram of water initially at body temperature of $37\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$ and $P=1\phantom{\rule{0.3em}{0ex}}\mathrm{atm}$ to water’s critical point temperature of $374\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$ and pressure of $218\phantom{\rule{0.3em}{0ex}}\mathrm{atm}$ , ${\rho}_{c}=0.315\phantom{\rule{0.3em}{0ex}}\mathrm{g}\u2215{\mathrm{cm}}^{3}$ is the critical point density of water. Since under short pulse laser illumination the melanosome is expected to reach temperatures of thousands of degrees Celsius (e.g., at $1\phantom{\rule{0.3em}{0ex}}\mathrm{J}\u2215{\mathrm{cm}}^{2}$ retinal fluence), the initial temperature and radius of the bubble would be determined by the rate of energy transfer to the melanosome. For laser pulses such that the laser transfers all its energy on time scales shorter than the characteristic bubble formation time, which is the most damaging case, the water in contact with the melanosome will be heated up to the maximum allowed thermodynamic state of liquid-vapor coexistence, the critical point. Numerically we can estimate the bubble formation time by considering that it takes a time scale of a few hundred picoseconds for the vapor shell around the melanosome to double its specific volume.

As was done by Gerstman, we assume that the laser energy is used to raise the temperature of the melanosome to
$374\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$
, and to bring a certain volume of water to the critical point [cf. Eq. 11 in Gerstman ^{19}]. The amount of water brought to the critical point is used to calculate the initial bubble radius. Because of the short time scales involved we further assume that the bubble will grow under adiabatic conditions and use the expression
$P{V}^{\gamma}=\text{constant}$
, with
$P$
and
$V$
the pressure and volume, respectively, and
$\gamma $
the ratio between the specific heat of the vapor at constant pressure and volume. Also in Eq. 12,
${c}_{v}=2.51\phantom{\rule{0.3em}{0ex}}\mathrm{J}\u2215\mathrm{g}\phantom{\rule{0.2em}{0ex}}\mathrm{K}$
is the specific heat of the melanosome,
${\rho}_{0}=1.35\phantom{\rule{0.3em}{0ex}}\mathrm{g}\u2215{\mathrm{cm}}^{3}$
is the static density of the melanosome,
$\Delta {T}_{c}=374\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}-37\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}=337\phantom{\rule{0.2em}{0ex}}\xb0\mathrm{C}$
is the temperature difference necessary to raise the melanosome to the critical point temperature of water,
${P}_{0}=218\phantom{\rule{0.3em}{0ex}}\mathrm{atm}$
is the critical point pressure, and
${P}_{\mathrm{min}}=1\phantom{\rule{0.3em}{0ex}}\mathrm{atm}$
is the ambient pressure, assumed to be the pressure inside the bubble at its maximum extension. The choice of
${P}_{\mathrm{min}}=1\phantom{\rule{0.3em}{0ex}}\mathrm{atm}$
was made by Gerstman ^{19} to balance two opposing effects: outward momentum of the liquid during bubble growth, which would cause an overshoot to
${P}_{\mathrm{min}}<1\phantom{\rule{0.3em}{0ex}}\mathrm{atm}$
inside the bubble and a larger
${R}_{\mathrm{max}}$
, versus viscous energy loss during expansion that would result in a smaller
${R}_{\mathrm{max}}$
. Using these values for the parameters, Eq. 12 for predicting bubble size as a function of fluence becomes

## 13

$$\frac{{R}_{\mathrm{max}}}{a}={\{1+488[\frac{C\left({\alpha}_{L}a\right){I}_{0}(\mathrm{J}\u2215{\mathrm{cm}}^{2})}{a\left(\mu \mathrm{m}\right)}-0.152]\}}^{1\u22153}.$$## Table 1

As a function of αLa ; Percent of incident energy absorbed C(αLa) , and scaling factor S(αLa)=C(1.0)∕C(αLa) for fluence to produce an identical shock front and bubble as for any of the results presented in this study if αLa is changed from 1.0 ( αL=10,000cm−1 , a=1μm ).

αLa | C(αLa) | S(αLa) |
---|---|---|

0.1 | 0.124 | 5.67 |

0.2 | 0.231 | 3.05 |

0.3 | 0.323 | 2.18 |

0.4 | 0.402 | 1.75 |

0.5 | 0.472 | 1.49 |

0.6 | 0.531 | 1.32 |

0.7 | 0.584 | 1.20 |

0.8 | 0.629 | 1.12 |

0.9 | 0.668 | 1.05 |

1.0 | 0.703 | 1.00 |

1.2 | 0.760 | 0.92 |

1.4 | 0.804 | 0.88 |

1.6 | 0.838 | 0.84 |

1.8 | 0.865 | 0.81 |

2.0 | 0.886 | 0.79 |

The curve labeled RPVS is calculated from the well-known Rayleigh-Plesset^{34} equation including viscosity and surface tension:

## 14

$$\frac{{P}_{B}\left(t\right)-{P}_{\infty}\left(t\right)}{{\rho}_{L}}=R\stackrel{\u0308}{R}+\frac{3}{2}{\stackrel{\u0307}{R}}^{2}+\frac{4{\nu}_{L}}{R}\stackrel{\u0307}{R}+\frac{2S}{{\rho}_{L}R},$$The difference between the Gerstman prediction and the RP line is due to the difference in the dynamics of the bubble growth. We see from Fig. 5 that the present numerical model and the Gerstman model both do a good job in explaining the experimental data, with the present numerical model better at lower fluences and the Gerstman model better at higher fluences. The RP line greatly overestimates the size of the bubble. The RPVS line does a better job, but also overestimates the bubble size. This overestimation by the RPVS model may be evidence that the growth of large bubbles is thermally controlled^{34} by heat loss, which the RPVS model does not take into account. The better agreement between the Gerstman model and the experiments is possibly due to not including inertial effects, which compensates for not allowing heat loss during bubble growth. The present numerical model explicitly includes heat loss. Another possible explanation for the overestimation of maximum bubble size in the Gerstman and RPVS models is that bubble nucleation occurs at temperatures and pressures below the critical point^{19} for pulses longer than the characteristic bubble formation time. A lower nucleation temperature would lead to a lower initial vapor pressure and subsequently to smaller bubble sizes.

## 5.

## Dependence of Shockwaves and Bubbles on Melanosome Properties

As discussed at the beginning of the previous section, important physical properties of melanosomes are not well known. To obtain the results presented in the previous sections, we chose specific values representing these properties. Since once they are measured, the values that we chose may turn out not to be those of melanosome, we also investigated how our results will change if the values of some of these parameters are changed. These results are also useful for investigations of interesting microparticles and nanoparticles other than melanosomes.

## 5.1.

### Dependence on ${\alpha}_{L}$

The effect of a different absorption coefficient
${\alpha}_{L}$
is straightforward to take into account, and the preceding results can be used with simple scaling. For an absorber with a constant radius, a different value of
${\alpha}_{L}$
merely means that a different fraction of the incident energy will be absorbed. The response of the system depends on the rate of energy absorbed as shown in Eq. 8, which depends on the product of the fluence
${I}_{0}$
and the fraction of energy absorbed
$C\left({\alpha}_{L}a\right)$
, as shown in Eq. 1. The response of the system does not distinguish between low fluence with strong absorption or high fluence with weak absorption. Therefore, if
${\alpha}_{L}$
for the absorber is smaller (or greater) than the
$\mathrm{10,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
used in this paper, the results presented above can be used if the fluences quoted above are multiplied by the scaling factors given in Table 1 so that the rate of energy absorption remains the same. These scaling factors,
$S\left({\alpha}_{L}a\right)$
for the fluence are calculated by evaluating
$C\left({\alpha}_{L}a\right)$
in Eq. 1 for any given
${\alpha}_{L}$
and comparing it to the value obtained using
${\alpha}_{L}=\mathrm{10,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, i.e.,
$S\left({\alpha}_{L}a\right)=C\left(1.0\right)\u2215C\left({\alpha}_{L}a\right)$
. For example, assuming
$a=1\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
, the shockwave and bubble produced by a laser of fluence
$500\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$
and
${\alpha}_{L}=\mathrm{10,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, will be equivalent to a fluence of only
$421\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$
if
${\alpha}_{L}=\mathrm{16,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, or a fluence of
$558\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$
if
${\alpha}_{L}=8000\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
. We emphasize that the scaling factor of Table 1 relates
${I}_{0}$
to
${\alpha}_{L}$
at constant absorber radius. If the absorber radius is changed, Table 1 cannot be used directly to scale
${I}_{0}$
to calculate shockwaves and bubbles.^{19}

## 5.2.

### Dependence on $B$

Since it does not occur as a simple factor as does
${\alpha}_{L}$
in
$C\left({\alpha}_{L}a\right)$
, the effect of
$B$
is not as simple to determine. For a specific absorber, such as the melanosomes in the retina, the prediction of the shockfront strength requires knowledge of the value of
$B$
. Alternatively, if absorbing particles are inserted into a cell to generate shockwaves for biomedical applications, the
$B$
of the particle material will affect the size of the shockfront as well as other characteristics of the pressure response such as
${\tau}_{c}$
. Thus, understanding the dependencies on
$B$
is important for guiding the choice of particle composition. The shock front strength depends on
$B$
in a nonlinear fashion, as seen in Fig. 6
. For six different pulse durations,
${\tau}_{L}=1\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
,
$10\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
,
$100\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
,
$250\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$
,
$1\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
, and
$10\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$
, we calculate how the pressure amplitude at
$u=2a$
depends on
$B$
for
${I}_{0}=100\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$
. We see that for a given
${\tau}_{L}$
, the value of
$B$
of the absorbing particle controls the strength of the pressure wave that is generated. There are several competing effects that depend on
$B$
. For
${\tau}_{L}<{\tau}_{c}$
, the shockwave pressure amplitude in the medium increases with
$B$
approximately^{11} as

## 15

$${P}_{\mathrm{liq}}\left(u\right)=\frac{\alpha B}{{c}_{v}}\left(\frac{{v}_{\mathrm{liq}}}{{v}_{\mathrm{liq}}+{v}_{\mathrm{abs}}}\right){\left(\frac{a}{u}\right)}^{\beta}C\left({\alpha}_{L}a\right){I}_{0},$$We found that the strength of the predicted shock front at $r=2a$ will be largest for an absorber with $B=80\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$ if $\alpha =2.4\times {10}^{-4}\phantom{\rule{0.3em}{0ex}}{\mathrm{K}}^{-1}$ . For a $10\text{-}\mathrm{ps}$ pulse, an absorber with $B=80\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$ will produce a shock front that is more than 40% stronger than an absorber with $B=20\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$ . As $B$ is increased above $80\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$ , the dropoff in the size of the expected shock front as a function of $B$ depends on the length of the laser pulse. We also see that the value of the absorber’s $\alpha $ plays an important role, which we investigate in the next section.

## 5.3.

### Dependence on the Thermal Expansion Coefficient

The strength of the shock front and the size of the bubble both depend on the value used for
$\alpha $
. We therefore investigated the effect of using an absorber with a different coefficient of thermal expansion
$\alpha $
. Figure 7(a)
shows the dependence of the shock font strength on
$\alpha $
, and Fig. 7(b) shows the dependence of the bubble size on
$\alpha $
. The values for other parameters used in the calculation for Fig. 7 were
${I}_{0}=880\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$
,
${\alpha}_{L}=\mathrm{10,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, and
$B=39.4\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$
. The dependence of both shockwaves and bubbles on
$\alpha $
is notably stronger^{11} for
${\tau}_{L}<{\tau}_{c}$
.

Another way of showing the strong dependence on $\alpha $ is presented in Fig. 8 , where the fluence required to produce a shock front of $1\phantom{\rule{0.3em}{0ex}}\mathrm{kbar}$ at $r=2a$ is plotted as a function of $\alpha $ . For Fig. 8 we set ${\tau}_{L}=0.1\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ , ${\alpha}_{L}=\mathrm{10,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , and $B=80\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$ . The curve in Fig. 8 displays a $1\u2215\alpha $ dependence, as expected from the relationship between $P$ , $\alpha $ , and ${I}_{0}$ , shown in Eq. 15.

Because of the uncertainty in the values of the physical properties of a melanosome, we determined a set of values that would produce the strongest shock front for a given fluence. We also required the values of the melanosome parameters to be reasonable based on existing evidence. We found that those values are $B=80\phantom{\rule{0.3em}{0ex}}\mathrm{GPa}$ and $\alpha =2.4\times {10}^{-4}\phantom{\rule{0.3em}{0ex}}{\mathrm{K}}^{-1}$ . Using these values for $B$ and $\alpha $ , we determined the fluence that would produce a shock front of $1000\phantom{\rule{0.3em}{0ex}}\text{bar}$ , and found that to be ${I}_{0}=120\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$ if ${\alpha}_{L}=\mathrm{10,000}\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ . This value of $120\phantom{\rule{0.3em}{0ex}}\mathrm{mJ}\u2215{\mathrm{cm}}^{2}$ decreases slightly for ${\tau}_{L}<100\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$ .

## 6.

## Summary

We have presented detailed results from calculations that show the expected size of shockwaves and bubbles generated by lasers pulses of various durations and fluences incident on a spherically absorbing melanosome immersed in an aqueous medium. The nonlinearity of the process requires careful attention to the numerical algorithm employed. The strength of the shock front drops rapidly as it expands outward, and eventualy decays to acoustic strength and speed. Whether or not these strong acoustic waves can still cause damage was not investigated.

We showed that stress confinement is valid for shockwaves in the surrounding medium when pulse durations are shortened below a characteristic absorber time. For melanosomes, this ${\tau}_{c}$ is likely to have an order of $100\phantom{\rule{0.3em}{0ex}}\mathrm{ps}$ . We also showed that cavitation confinement for bubbles also occurs, and that bubbles become relevant for pulse durations shorter than $100\phantom{\rule{0.3em}{0ex}}\mathrm{ns}$ . This difference in confinement time scales is a reflection of the difference in the dynamics; vaporization propagates outwards at speeds that are approximately $1\u221550$ that of pressure waves.

A table was presented that gave a factor by which the fluence should be scaled to get the same pressure and bubble results for absorbers with different coefficients of absorption. We also showed that the size of the shock front that is generated depends on the physical properties of the absorbing melanosome, such as its bulk modulus and thermal coefficient of expansion. The importance of these properties makes it crucial for them to be measured independently.

The results of this paper can be used to guide biomedical procedures involving the retina or other tissues with strongly absorbing particles, and experiments that investigate shockwave and bubble formation from spherical absorbers in biological tissues or other materials. The results can also be used as a basis for determining the biological molecules or structures in a cell that are most important when an incident laser pulse causes damage.

Finally, note that the previous discussion bears relevance to understanding the transition between bulk and nonbulk properties in materials. From Eq. 11 we see that for a constant bulk modulus, the critical acoustic time depends linearly on the radius of the absorber. However, as the size of the absorber is reduced, we expect to eventually reach a regime where finite size effects will influence the bulk modulus of the material. Hence, as the size of the absorber is made smaller, and as one approaches the nonbulk regime, where finite size effects become important, one would observe a deviation from the linear relationship expressed in Eq. 11, when measuring the oscillation time of the absorber.

## Acknowledgments

The authors would like to thank the Air Force Office of Scientific Research (AFOSR) for funding through Grant No. F49620-03-1-0221.