1 November 2005 Biophysical effects of pulsed lasers in the retina and other tissues containing strongly absorbing particles: shockwave and explosive bubble generation
Author Affiliations +
Damage by pulsed lasers to the retina or other tissues containing strongly absorbing particles may occur through biophysical mechanisms other than simple heating. Shockwaves and bubbles have been observed experimentally, and depending on pulse duration, may be the cause of retinal damage at threshold fluence levels. We perform detailed calculations on the shockwave and bubble generation expected from pulsed lasers. For a variety of different laser pulse durations and fluences, we tabulate the expected strength of the shockwave and size of the bubble that will be generated. We also explain how these results will change for absorbing particles with different physical properties such as absorption coefficient, bulk modulus, or thermal expansion coefficient. This enables the assessment of biological danger, and possible medical benefits, for lasers of a wide range of pulse durations and energies, incident on tissues with absorbing particles with a variety of thermomechanical characteristics.



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 system3, 4, 5, 6 and other biological tissues.7 Absorption of the laser energy generates various responses that include temperature rise8, 9, 10 (thermal), pressure11, 12, 13, 14, 15, 16, 17, 18 (acoustomechanical), and vaporization17, 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 I0 for biological danger to the cell varies as a function of pulse duration. τ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 τL are quite different, the mechanism responsible for threshold damage may change as a function of τL . For long τL , thermal damage requires the lowest I0 and is the threshold damage mechanism. As τL shortens, thermal damage requires less fluence until τL decreases to a characteristic heat conduction time τH , that depends on absorber size and the spacing between absorbers, as explained in Ref. 9. For τL<τH , the threshold fluence remains constant for thermal damage. The damage threshold fluence for explosive vaporization also decreases as τL shortens, but continues to decrease for τL<τH . The characteristic vaporization bubble time τV is equal to the time scale of bubble expansion and for the system treated here is of the order of 100ns . 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 τV , the fluence required for damage from explosive vaporization remains approximately constant. Since τV<τH , for τL between τH and τ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 τL decreases and becomes constant when τL<τc , where τc is a characteristic time for shockwave production. Since τc<τV , for τL between τV and τ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 τL enters the regions between the characterisitic times τH , τV , and τ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 τV and τc , as well as determining how the strength and extent of vaporization and shockwave effects depend on τL and I0 .

Fig. 1

Schematic diagram for damage thresholds as a function of pulse duration. The dashed lines represent suprathreshold damage mechanisms, e.g., for long pulse durations, damage-producing bubbles can occur, however, they occur at fluences above the threshold for thermal damage.





Governing Equations for Absorbing Melanosome

Our model11, 12 consists of a uniform spherical absorber surrounded by a transparent medium. The rate of energy input per unit mass is given by19


where I0 is the incident laser fluence in joules per centimeter squared, a is the radius of the absorbing sphere, τL is the laser pulse duration, ρ0 is the static density of the sphere, and αL is the absorption coefficient of the absorbing sphere. As the particle size shrinks to the wavelength of the incident light, the effective αL may change. The fraction of energy incident on the sphere that is absorbed is given by C(αLa) and is the term in the curley braces, as first defined in Ref. 19. We assume homogeneous energy deposition throughout the spherical absorber. This is valid if the absorber is a good thermal conductor (a24)(cvρλ)<τL , enabling quick internal energy conduction and temperature equilibration, and if αa<1 , such that a considerable amount of the laser energy reaches the portion of the melanosome on the far side from the laser beam. For the ultrashort pulses and melanosome characteristics used in this paper, these conditions are approximately valid. The change in energy deposition as a function of melanosome properties is discussed later in Sec. 5.

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 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,,N . The position of each mass element is denoted by Lagrangian coordinates ui=ui(t) . The numerical value of each ui gives the location of element i in terms of the coordinate r . Each element i is given an initial position ui(t=0)=ri . 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,,n are inside the absorber, and i=n+1,,N are assigned properties of the aqueous medium.

The mathematical dot operation ḟ(t) means a total time derivative, the spatial derivative is taken with respect to r , while the spatial derivative with respect to u at a specific mass point is specifically denoted by u . With this notation, the equation of motion of a point inside the absorbing sphere is


where P is the pressure and ρ=ρ(t) is the time-varying density, which is related to the static density by mass conservation


where u is the radial component of the vector u . In spherical geometry, the equation of motion can be expressed as



If we assume that the absorbing melanosome has a constant bulk modulus B and constant thermal expansion coefficient α , the equation of state (EOS) of the absorber can be written as


where v=1ρ is the specific volume and is related to u by


This enables us to write Eq. 5 as



Energy conservation of the absorber relates the rate of absorption to temperature rise, volume change, and heat lost through conduction to the surrounding medium:


where cv is the specific heat, λ is the thermal conductivity of the absorber, and s is the specific entropy. For the absorber, we do not calculate s , and therefore only use the term on the right hand side of Eq. 8, based on the EOS given by Eq. 7.


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, İ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 pressures22, 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 647K and 218atm . Our treatment enables us to predict the temperature, pressure, and bubble size created by a laser over a wide range of pulse durations (0.1ps<τL<1μs) and fluences (I0<10Jcm2) . 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.


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μm , ρ0=1.35gcm3 , and cv=2.51JgK . We use an absorption coefficient of αL=10,000cm1 . 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 αL . Reliable numbers for the bulk modulus B and bulk thermal coefficient of expansion α 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.4GPa and α=2.98×105K1 . Later, we discuss how the results will change for other values of B and α . The surrounding medium is aqueous with ρ0=1.0gcm3 , and thermal conductivity λ=5.56×103JcmKs . Since for water we are using a tabulation relating P , v , T , and s , we do not use values for cv , B , and α 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 τL=1ps , I0=1.0Jcm2 (suprathreshold for both shockwaves and bubble formation to elucidate the effects) with αL=10,000cm1 . The melanosome undergoes a ΔT=1595K , and the shock front in the near field can be as large as 4kbar , corresponding to a shock front speed of Us=2630ms ( 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.

Fig. 2

Refer to text for laser and absorber properties. (a) Multiple snapshots of the pressure profile outside the melanosome in the aqueous medium at various times after the single laser pulse, including trailing shock fronts due to the remnant ringing of the absorber, (b) decreasing strength of the leading shock front as it propagates out into the medium, the decay is faster than the 1r dependence of an acoustic wave in three dimensions; and (c) log-log plot of (b) to emphasize the changing rate of P decay as a function of P .


Fig. 3

For the same laser pulse as used in Fig. 2, the expansion of the bubble outward from the absorber’s surface at r=a . The bubble is defined by the region of high specific volume v . The time scale for bubble expansion is approximately 50 times larger than that for the pressure wave propagation displayed in Fig. 2. Note that the decreased specific volume in the proximity of the melanosome results from the increased pressure of the vapor due to heat transfer from the hot melanosome.



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 used12 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 1r 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 β=β(P)


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 β and we give a few representative values.


Absorber’s Characteristic Acoustic Time τ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 , α , and ρ 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 τc=2ac , where c=Bρ is the speed of sound in the absorber. This results in


Using graphite’s B=39.4GPa with melanosome’s ρ0=1.35gcm3 , then for an absorbing sphere of a=1μm Eq. 11 gives τc=370ps .

If B and ρ of the absorber are known, then τ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 τ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.


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 initiated32 by heating the liquid water to above 100°C . The region of the aqueous medium farther away, that has a specific volume of approximately 1cm3g , 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μm less than this bubble radius. We also note that over the initial 100ps 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.


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μ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 I0 is almost identical for τL=1ps and 100ps , but drops dramatically for τL=370ps . This τL is also τc for the absorber, as already explained, and defines the stress confinement time. For τL<τc , the pressure effects in the surrounding medium depend on I0 , but are independent of τL . (As explained in Ref. 11, the concept of stress confinement is not valid in the core region of the absorber). For τL> τc , the amplitude of the leading pressure front decreases inversely with τ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 τL<τc . Therefore, pulses with τL<τc will cause much greater pressure induced effects in the surrounding medium than pulses with τL> τc .

Fig. 4

Calculated results for pressure amplitude and maximum bubble radius for laser pulses of different fluence and pulse duration. Values used for properties of the absorbing melanosome are αL=10,000cm1 , B=39.4GPa , and α=2.98×105K1 . (a) The strength (pressure jump) of the expected shock front 1μm outside the absorbing melanosome (r=2a) , and (b) The maximum bubble size, as a function of fluence for different pulse durations.


Figure 4(b) shows the maximum bubble size Rmax as a function of I0 for different τL . The bubble size decreases by only 10% as the pulse duration increases from 1ps to 1ns , which is a manifestation of “cavitation confinement.”12 For τL=100ns , there is a significant decrease in the bubble size, and this defines the characteristic cavitation or vaporization confinement time τv , which is more than two orders of magnitude larger than τ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, τ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 Rmax as a function of I0 calculated from our numerical model with experimental measurements, as well as predictions from two other models. The curve labeled 1ps employs the model of this paper, and is the same 1ps curve plotted in Fig. 4(b). The curve labeled 10ns is of the same model but uses a pulse duration of 10ns . The curve labeled NB is the maximum bubble radius as a function of I0 obtained by fitting the experimental data of Neumann and Brinkmann33 for a pulse duration of 12ns . We found that the curve [Rmax(I0)1μm]=[0.3+44(I00.12)]0.5 , for I0> 0.12 in Jcm2 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 model19 that makes simple assumptions about bubble growth but has the advantage of giving a straightforward analytical expression for Rmax as a function of I0



Fig. 5

Comparison of the maximum bubble radius obtained from model and experiment. The curves labeled 1ps and 10ns employ the numerical model of this paper for pulse durations of 1ps and 10ns , respectively. The curve labeled NB represents a best-fit curve to the experimental data of Neumann and Brinkmann.33 The curve labeled Gerstman is calculated from Eq. 13. The curve labeled RP is calculated from Eq. 14 without the viscous and surface tension terms, and the curve labeled RPVS is calculated from Eq. 14 assuming ηL=ρLυL=500μPas and S=0.07Nm .


We use the following conditions for the parameters in Eq. 12, their choice is more fully explained in Ref. 19. The parameter q=2770Jg is the heat needed to raise one gram of water initially at body temperature of 37°C and P=1atm to water’s critical point temperature of 374°C and pressure of 218atm , ρc=0.315gcm3 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 1Jcm2 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°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 PVγ=constant , with P and V the pressure and volume, respectively, and γ the ratio between the specific heat of the vapor at constant pressure and volume. Also in Eq. 12, cv=2.51JgK is the specific heat of the melanosome, ρ0=1.35gcm3 is the static density of the melanosome, ΔTc=374°C37°C=337°C is the temperature difference necessary to raise the melanosome to the critical point temperature of water, P0=218atm is the critical point pressure, and Pmin=1atm is the ambient pressure, assumed to be the pressure inside the bubble at its maximum extension. The choice of Pmin=1atm was made by Gerstman 19 to balance two opposing effects: outward momentum of the liquid during bubble growth, which would cause an overshoot to Pmin<1atm inside the bubble and a larger Rmax , versus viscous energy loss during expansion that would result in a smaller Rmax . Using these values for the parameters, Eq. 12 for predicting bubble size as a function of fluence becomes


Using C(αLa)=0.703 from Table 1 for αL=10,000cm1 and a=1μm , Eq. 13 predicts that bubble formation of appreciable size requires a minimum I0=0.216Jcm2 , in excellent agreement with the value measured experimentally in Ref. 33, and with our numerical simulations.

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)

The curve labeled RPVS is calculated from the well-known Rayleigh-Plesset34 equation including viscosity and surface tension:


where PB(t) is the time varying pressure inside the bubble; P(t) is the ambient pressure far away from the bubble surface, which we set equal to 1bar and do not vary; ηL=ρLνL=500μPas is the dynamic viscosity; S=0.07Nm is the surface tension; and ρL=1gcm3 is the static density of water. The curve labeled RP uses Eq. 14 but without the terms representing viscosity and surface tension. The dependence of Rmax on I0 is determined by the dependence of the bubble’s initial radius R(t=0) , and pressure PB(t=0) on I0 , as well as assumptions about PB(t) . These values are determined in the same way as for the line labeled Gerstman, as already explained. We found that the solutions to the Rayleigh-Plesset equation vary by less than 1% for reasonable values for the initial bubble wall velocity.

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 controlled34 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 point19 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.


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.


Dependence on αL

The effect of a different absorption coefficient α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 α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 I0 and the fraction of energy absorbed C(αLa) , 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 αL for the absorber is smaller (or greater) than the 10,000cm1 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(αLa) for the fluence are calculated by evaluating C(αLa) in Eq. 1 for any given αL and comparing it to the value obtained using αL=10,000cm1 , i.e., S(αLa)=C(1.0)C(αLa) . For example, assuming a=1μm , the shockwave and bubble produced by a laser of fluence 500mJcm2 and αL=10,000cm1 , will be equivalent to a fluence of only 421mJcm2 if αL=16,000cm1 , or a fluence of 558mJcm2 if αL=8000cm1 . We emphasize that the scaling factor of Table 1 relates I0 to αL at constant absorber radius. If the absorber radius is changed, Table 1 cannot be used directly to scale I0 to calculate shockwaves and bubbles.19


Dependence on B

Since it does not occur as a simple factor as does αL in C(αLa) , 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 τ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, τL=1ps , 10ps , 100ps , 250ps , 1ns , and 10ns , we calculate how the pressure amplitude at u=2a depends on B for I0=100mJcm2 . We see that for a given τ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 τL<τc , the shockwave pressure amplitude in the medium increases with B approximately11 as


where vliq and vabs are the acoustic speed of pressure waves in the liquid medium and absorber, respectively. However, as B increases, τc decreases, as displayed in Eq. 11. For large enough B , τc drops below τL and the shockfront strength then decreases with increasing B , and softens into an acoustic wave. These two competing effects lead to a peak in the shockfront amplitude as a function of B for the constant τL curves in Fig. 6. The graphs in Fig. 6 use I0=100mJcm2 , αL=10,000cm1 , and two different value for α : Fig. 6(a) uses α=2.98×105K1 , and Fig. 6(b) uses α=2.4×104K1 .

Fig. 6

Dependence of the pressure response on the mechanical parameter B . Maximum pressure at r=2μm for (a) α=2.98×105K1 and (b) α=2.4×104K1 . For large α and short pulse durations the curves peak around B80GPa . Values used for laser and absorber properties: I0=100mJcm2 ; αL=10,000cm1 . The labels for the curves in (b) are the same as for the curves (a).


We found that the strength of the predicted shock front at r=2a will be largest for an absorber with B=80GPa if α=2.4×104K1 . For a 10-ps pulse, an absorber with B=80GPa will produce a shock front that is more than 40% stronger than an absorber with B=20GPa . As B is increased above 80GPa , 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 α plays an important role, which we investigate in the next section.


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 α . We therefore investigated the effect of using an absorber with a different coefficient of thermal expansion α . Figure 7(a) shows the dependence of the shock font strength on α , and Fig. 7(b) shows the dependence of the bubble size on α . The values for other parameters used in the calculation for Fig. 7 were I0=880mJcm2 , αL=10,000cm1 , and B=39.4GPa . The dependence of both shockwaves and bubbles on α is notably stronger11 for τL<τc .

Fig. 7

Dependence of the response on the thermal expansion coefficient α for (a) maximum pressure at r=2μm and (b) maximum bubble radius as a function of α for different pulse durations. Values used for laser and absorber properties: I0=0.88Jcm2 , αL=10,000cm1 , and B=39.4GPa .


Another way of showing the strong dependence on α is presented in Fig. 8 , where the fluence required to produce a shock front of 1kbar at r=2a is plotted as a function of α . For Fig. 8 we set τL=0.1ns , αL=10,000cm1 , and B=80GPa . The curve in Fig. 8 displays a 1α dependence, as expected from the relationship between P , α , and I0 , shown in Eq. 15.

Fig. 8

Critical fluence Ic required to obtain a shock front of 1kbar at a position of 1μm outside the surface of the melanosome (r=2μm) as a function of the thermal expansion coefficient of the melanosome α . Values used for other parameters: αL=10,000cm1 , B=80GPa , and τL=0.1ns . For reference we also give a best fit curve showing the 1α dependence of the critical fluence.


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=80GPa and α=2.4×104K1 . Using these values for B and α , we determined the fluence that would produce a shock front of 1000bar , and found that to be I0=120mJcm2 if αL=10,000cm1 . This value of 120mJcm2 decreases slightly for τL<100ps .



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 τc is likely to have an order of 100ps . We also showed that cavitation confinement for bubbles also occurs, and that bubbles become relevant for pulse durations shorter than 100ns . This difference in confinement time scales is a reflection of the difference in the dynamics; vaporization propagates outwards at speeds that are approximately 150 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.


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



Hirsch L. R., Stafford R. J., Bankson J. A., Sershen S. R., Rivera B., Price R. E., Hazle J. D., Halas N. J., West J. L., “Nanoshell-mediated near-infrared thermal therapy of tumors under magnetic resonance guidance,” Proc. Natl. Acad. Sci. U.S.A., 100 (23), 13549 –13554 (2003). 10.1073/pnas.2232479100 0027-8424 Google Scholar


Zharov V. P., Galitovskaya E., Viegas M., “Photothermal guidance for selective photothermolysis with nanoparticles,” Proc. SPIE, 5319 (1), 291 –300 (2004). 10.1117/12.532011 0277-786X Google Scholar


Sliney D. H., “Retinal injury from laser radiation,” Nonlinear Opt., 21 (1), 1 –18 (1999). 1053-3729 Google Scholar


Rockwell B., “Ultrashort laser pulse bioeffects and safety,” J. Laser Appl., 11 (1), 42 –44 (1999). 1042-346X Google Scholar


Thomas R. J., Noojin G. D., Stolarski D. J., Hengst G. T., Toth C. A., Roach W. P., Rockwell B. A., “Retinal damage from femtosecond to nanosecond laser exposure,” Proc. SPIE, 3902 (1), 54 –61 (2000). 0277-786X Google Scholar


Johnson T. E., Keeler N., Dennis J. E., Figueroa C. L., Press H. A., Rockwell R. J., Stuck B. E., Roach W. P., Wartick A. L., “Investigation of laser injuries,” Proc. SPIE, 4953 (1), 61 –69 (2003). 0277-786X Google Scholar


Vogel A., Venugopalan V., “Mechanisms of pulsed laser ablation of biological tissue,” Chem. Rev. (Washington, D.C.), 103 (2), 577 –644 (2003). 10.1021/cr010379n 0009-2665 Google Scholar


Gerstman B. S., Glickman R. D., “Activated Rate Processes and a Specific Biochemical Mechanism for Explaining Delayed Laser Induced Thermal Damage to the Retina,” J. Biomed. Opt., 4 (3), 345 –351 (1999). 10.1117/1.429936 1083-3668 Google Scholar


Thompson C. R., Gerstman B. S., Jacques S. L., Rogers M. E., “Melanin Granule Model for Laser-Induced Thermal Damage in the Retina,” Bull. Math. Biol., 58 (3), 513 –553 (1996). 10.1016/0092-8240(95)00360-6 0092-8240 Google Scholar


Freund D. E., Sliney D. H., “Dependence of retinal model temperature calculations on beam shape and absorption coefficients,” Lasers Life Sci., 8 (3), 229 –247 (1999). 0886-0467 Google Scholar


Sun J. M., Gerstman B. S., “Photoacoustic generation for a spherical absorber with impedance mismatch with the surrounding media,” Phys. Rev. E, 59 (5), 5772 –5789 (1999). 10.1103/PhysRevE.59.5772 1063-651X Google Scholar


Sun J. M., Gerstman B. S., Li B., “Bubble dynamics and shock waves generated by laser absorption of a photoacoustic sphere,” J. Appl. Phys., 88 (5), 2352 –2362 (2000). 10.1063/1.1288507 0021-8979 Google Scholar


Zhigilei L. V., Garrison B. J., “Computer simulation study of damage and ablation of submicron particles from short pulse laser irradiation,” Appl. Surf. Sci., 127–129 (1), 142 –150 (1998). 0169-4332 Google Scholar


Zhigilei L. V., Garrison B. J., “Microscopic simulation of short pulse laser damage of melanin particles,” Proc. SPIE, 3254 (1), 135 –143 (1998). 0277-786X Google Scholar


Vogel A., Busch S., Parlitz U., “Shock wave emission and cavitation bubble generation by picosecond and nanosecond optical breakdown in water,” J. Acoust. Soc. Am., 100 (1), 148 –165 (1996). 10.1121/1.415878 0001-4966 Google Scholar


Paltauf G., Schmidt-Kloiber H., “Internal photomechanical fracture of spatially limited absorbers irradiated by short laser pulses,” Proc. SPIE, 3254 (1), 112 –120 (1998). 10.1117/12.308203 0277-786X Google Scholar


Lin C. P., Kelly M. W., “Cavitation and acoustic emission around laser-heated microparticles,” Appl. Phys. Lett., 72 (22), 2800 –2802 (1998). 10.1063/1.121462 0003-6951 Google Scholar


Faraggi E., Wang S., Gerstman B., “Stress confinement, shock wave formation, and laser-induced damage,” Proc. SPIE, 5695 (1), 209 –215 (2005). 0277-786X Google Scholar


Gerstman B., Thompson C., Jacques S., Rogers M. E., “Laser induced bubble formation in the retina,” Lasers Surg. Med., 18 (1), 10 –21 (1996). 10.1002/(SICI)1096-9101(1996)18:1<10::AID-LSM2>3.0.CO;2-U 0196-8092 Google Scholar


Strauss M., Amendt P. A., London R. A., Maitland D. J., Glinsky M. E., Lin C. P., Kelly M. W., “Computational modeling of stress transient and bubble evolution in short-pulse laser-irradiated melanosome particles,” Proc. SPIE, 2975 (1), 261 –270 (1997). 0277-786X Google Scholar


Haar L., Gallagher J. S., Kell G. S., NBS/NRC Steam Tables, Hemisphere Publishing Corp., New York (1984). Google Scholar


Paltauf G., Reichel E., Schmidt-Kloiber H., “Study of different ablation models by use of high-speed-sampling photography,” Proc. SPIE, 1646 (1), 343 –352 (1992). 0277-786X Google Scholar


Paltauf G., Schmidt-Kloiber H., “Model study to investigate the contribution of spallation to pulsed laser ablation of tissue,” Lasers Surg. Med., 16 (3), 277 –287 (1995). 0196-8092 Google Scholar


Marschall H. B., Morch K. A., Keller A. P., Kjeldsen M., “Cavitation inception by almost spherical solid particles in water,” Phys. Fluids, 15 (2), 545 –553 (2003). 10.1063/1.1535940 1070-6631 Google Scholar


Hirshfelder J. O., Curtiss C. F., Byron Bird R., Molecular Theory of Gases and Liquids, John Wiley, New York (1964). Google Scholar


Ridah S., “Shock waves in water,” J. Appl. Phys., 64 (1), 152 –158 (1988). 10.1063/1.341448 0021-8979 Google Scholar


Glickman R. D., Jacques S. L., Hall R. M., Kumar N., “Revisiting the internal absorption coefficient of the retinal pigment epithelium melanosome,” Proc. SPIE, 4257 (1), 134 –141 (2001). 0277-786X Google Scholar


Lin S.-P., Wang L.-H., Jacques S. L., Tittel F. K., “Measurement of tissue optical properties using oblique incidence optical fiber reflectometry,” Appl. Opt., 36 (1), 136 –143 (1997). 0003-6935 Google Scholar


Sardar D., Mayo M., Glickman R., “Optical characterization of melanin,” J. Biomed. Opt., 6 (4), 404 –411 (2001). 10.1117/1.1411978 1083-3668 Google Scholar


Reynolds W. N., Physical Properties of Graphite, Elsevier Pub. Co., Amsterdam (1968). Google Scholar


Press W. H., Flannery B. P., Teukolsky S. A., Vetterling W. T., Numerical Recipes, 628 Cambridge University Press, New York (1986). Google Scholar


Neumann J., Brinkmann R., “Boiling nucleation on melanosomes and microbeads transiently heated by nanosecond and microsecond laser pulses,” J. Biomed. Opt., 10 (2), 024001 (2005). 1083-3668 Google Scholar


Neumann J., Brinkmann R., “Microbubble dynamics around laser heated microparticles,” Proc. SPIE, 5142 (1), 82 –87 (2003). 10.1117/12.499944 0277-786X Google Scholar


Brennen C. E., Cavitation and Bubble Dynamics, 34 Oxford University Press, New York (1995). Google Scholar
© (2005) Society of Photo-Optical Instrumentation Engineers (SPIE)
Eshel Faraggi, Eshel Faraggi, Bernard S. Gerstman, Bernard S. Gerstman, Jinming Sun, Jinming Sun, "Biophysical effects of pulsed lasers in the retina and other tissues containing strongly absorbing particles: shockwave and explosive bubble generation," Journal of Biomedical Optics 10(6), 064029 (1 November 2005). https://doi.org/10.1117/1.2139970 . Submission:

Back to Top