A Method for Optodynamic Characterization of Erbium Laser Ablation Using Piezoelectric Detection

The Er:YAG laser, with a wavelength of 2.94 μm [1], is a well-established tool in medicine and surgery, particularly in dentistry [2] and dermatology [3]. Its infrared light is absorbed strongly in water and hydroxyapatite, providing effective laser ablation of soft and hard biological tissues [4]. Numerous new medical treatments (particularly in osteotomy) would also benefit from the advantages of Er:YAG laser tissue interaction over the conventional methods (e.g. noncontact intervention, smaller heat-affected zone, absence of mechanical vibration). One of the main technical and scientific challenges yet to be solved is the development of viable and reliable systems for the on-line monitoring of the key parameters, such as the cutting depth and the type of the removed tissue [5] to [8]. Tissue ablation with the Er:YAG laser is driven by microexplosions. Absorbed laser energy is partially released in the form of nonlinear mechanical waves propagated in the ablated tissue and in the air above it. Following this, material ejection occurs. These so-called optodynamic phenomena have received considerable attention in the context of characterization of laser ablation. Various set-ups have been examined for this purpose: spatially resolved techniques, such as schlieren [9] and [10], shadowgraphy [11] to [14] or holography [15], as well as the temporally resolved ones: laser interferometer [16] to [18], laser beam deflection probe [19] and [20] and capacitive or piezoelectric acoustic sensors [7], [8] and [21]. While most of these techniques represent useful research tools within controlled laboratory experiments, only a few of them exhibit the potential to be used for on-line process monitoring in real medical applications. Here, additional technical factors come into prominence: e.g. compactness, affordability, insensitivity to environmental influences, speed of response, etc. In our view, broadband piezoelectric acoustic sensors are the means to make the on-line process monitoring of Er:YAG laser ablation practicable. Fig. 1 illustrates the idea: a piezoelectric sensor is attached to the laser handpiece to detect shock waves in the air above the irradiated tissue. Considering typical conditions (geometry of the handpiece and sensor, focal length of the focusing optics, laser pulse energy, etc.) and previous research results [11] it is reasonable to assume that the shock waveform is nearly hemispherical as it impinges onto the sensor and its amplitude is decreased to the intermediate-to-weak shock level. Piezoelectric shock sensors that respond linearly to this kind of excitation are available. Existing methods of shock wave characterization using acoustic sensors rely on empirically selected signal features, such as acoustic signal energy, signal amplitude and time of flight [7], [21] and [22]. Relationships of these signal features to other influencing factors (e.g. sensor characteristics, orientation or distance from the source) are not known, A Method for Optodynamic Characterization of Erbium Laser Ablation Using Piezoelectric Detection Bosiger, G. – Perhavec, T. – Diaci, J. Georgije Bosiger1,* – Tadej Perhavec1 – Janez Diaci2 1 Fotona d.d. Ljubljana, Slovenia 2 University of Ljubljana, Faculty of Mechanical Engineering, Slovenia


INTRODUCTION
The Er:YAG laser, with a wavelength of 2.94 μm [1], is a well-established tool in medicine and surgery, particularly in dentistry [2] and dermatology [3].Its infrared light is absorbed strongly in water and hydroxyapatite, providing effective laser ablation of soft and hard biological tissues [4].Numerous new medical treatments (particularly in osteotomy) would also benefit from the advantages of Er:YAG laser tissue interaction over the conventional methods (e.g.noncontact intervention, smaller heat-affected zone, absence of mechanical vibration).One of the main technical and scientific challenges yet to be solved is the development of viable and reliable systems for the on-line monitoring of the key parameters, such as the cutting depth and the type of the removed tissue [5] to [8].
Tissue ablation with the Er:YAG laser is driven by microexplosions.Absorbed laser energy is partially released in the form of nonlinear mechanical waves propagated in the ablated tissue and in the air above it.Following this, material ejection occurs.These so-called optodynamic phenomena have received considerable attention in the context of characterization of laser ablation.Various set-ups have been examined for this purpose: spatially resolved techniques, such as schlieren [9] and [10], shadowgraphy [11] to [14] or holography [15], as well as the temporally resolved ones: laser interferometer [16] to [18], laser beam deflection probe [19] and [20] and capacitive or piezoelectric acoustic sensors [7], [8] and [21].While most of these techniques represent useful research tools within controlled laboratory experiments, only a few of them exhibit the potential to be used for on-line process monitoring in real medical applications.Here, additional technical factors come into prominence: e.g.compactness, affordability, insensitivity to environmental influences, speed of response, etc.
In our view, broadband piezoelectric acoustic sensors are the means to make the on-line process monitoring of Er:YAG laser ablation practicable.Fig. 1 illustrates the idea: a piezoelectric sensor is attached to the laser handpiece to detect shock waves in the air above the irradiated tissue.Considering typical conditions (geometry of the handpiece and sensor, focal length of the focusing optics, laser pulse energy, etc.) and previous research results [11] it is reasonable to assume that the shock waveform is nearly hemispherical as it impinges onto the sensor and its amplitude is decreased to the intermediate-to-weak shock level.Piezoelectric shock sensors that respond linearly to this kind of excitation are available.
Existing methods of shock wave characterization using acoustic sensors rely on empirically selected signal features, such as acoustic signal energy, signal amplitude and time of flight [7], [21] and [22].Relationships of these signal features to other influencing factors (e.g.sensor characteristics, orientation or distance from the source) are not known, thus limiting the applicability of these methods to strictly controlled conditions.In this paper we propose another approach that opens the way to on-line process monitoring.We develop a new model of the sensor that takes into account the relative position and orientation of the sensor as well as its mechanical and electrical properties.We verify the model by comparing signals, detected at different sensor distances and orientations relative to the ablated spot, with the theoretical waveforms determined from a numerical solution of the point explosion model [23], taking into account the model of the piezoelectric sensor.Observing excellent agreement between the theoretical and experimental waveforms, we propose a novel method that employs shock-wave energy released during the ablation process as a process characteristic that is almost independent of the exact geometrical properties of the detection set-up.

THEORY
In this section, a new model of a piezoelectric shock sensor is presented.The finite size of the sensor and its mechanical and electrical characteristics are taken into account.Next, a numerical procedure is described that allows the determination of theoretical pressure waves at the sensor surface employing the Taylor-Sedov point explosion model.

Piezoelectric Sensor Model
The theoretical description is simplified on the assumption that the incident shock wave, during the propagation over the piezoelectric transducer surface, can be approximated as a linear spherical acoustic wave [24].For a Taylor-Sedov explosion, this applies well for an intermediate or weak shock over a small distance of propagation.With this assumption, the pressure transient at a given point r s on the transducer can be expressed by the pressure in a centre point r 0 of the transducer, using the equation: where u 0 is the shock propagation speed which is determined using the point explosion model, and where it is expressed as a function of the distance from the source P and the released energy: u 0 = u 0 (r 0 ,E).We also disregard all the effects that result due to a change of the acoustic impedance, assuming that these affect the amplitude and not the waveform.With this assumption we write: where F a (t;r 0 ) is the force that acts on the sensor in its axial direction, and Δp(t,r s ) is the pressure transient at a point on surface A s of the sensor (see Fig. 2).Vector r 0 points to the sensor centre and can be expressed by the horizontal s 0 and vertical h 0 distance.Inserting Eq. ( 1) into (2) and substituting Δp with a spherical impulse disturbance δ, we get Rayleigh's integral [25]: The above equation represents the response of a sensor with a finite aperture to a spherical impulse disturbance.Assuming a cylindrical sensor with a flat aperture, it is possible to express h(t;r 0 ) in terms of the angle Θ s (t;r 0 ) determined by the intersection of the projected spherical wave onto the sensor aperture as shown on Fig. 2: Using the superposition principle, we find that the force acting on a sensor of finite size is proportional to the convolution of a pressure transient Δp(t,r 0 ) in the centre of the transducer and the signal waveform Θ s (t;r 0 ) representing the impulse response of the finite sensor aperture: In order to take into account the electrical and mechanical characteristics, we treat the sensor as a one-dimensional element, where deformation that generates electric current acts only in its axial direction.The basis for this assumption represents the design of the piezoelectric acoustic sensor, where the sensor housing and the insulation suppresses the radial excitations (see Fig. 2).Pure elastic deformations of the sensing element are assumed.With these assumptions, linear static and dynamic characteristics apply.
The Piezoelectric transducer is represented with a current source and capacitor C s in parallel [26].Parallel resistance of the sensing element is usually large and can be neglected.Taking into account the capacitance of the cable C c and the resistive load R L , we get a transfer function for a linear piezoelectric force sensing element, Eq. ( 6), where τ = R L (C s +C c ) is the time constant, d the piezoelectric charge constant and V the measured voltage.H(s) represents the dynamic characteristics of the elastic piezoelectric structure, where natural frequencies can be determined by measuring electrical impedance on an impedance analyzer [26].
By analyzing the transfer function in Eq. ( 6) we can conclude that the sensor's sensitivity in steady state is inversely proportional to the total capacitance.Due to this capacitance, the sensor behaves as a high pass RC filter with time constant τ.The sensor output is in the frequency region between 1/(2πτ) and the first natural frequency of the piezoelectric element is proportional to the force F a .

Shock-Wave Propagation Model
We employ the Taylor-Sedov point explosion model to model propagation of the spherical shock wave [23].
The model assumes that a finite amount of energy E is released instantly in an infinitesimal volume of a perfect gas.Propagation of the blast wave is described by a set of hyperbolic partial differential equations in the Euler form: where U is the vector of the conserved variables, F = F(U) their fluxes and S(U) the geometric source term that results from the transformation of the Euler equations to the spherical coordinates.Subscripts denote partial derivatives with respect to the independent variables; time t and radius r.Primitive variables are the mass density ρ(r,t), fluid speed v(r,t) and pressure p(r,t).Propagation at shock wave-front r = r s (t) is governed by the Rankine-Hugoniot jump conditions.
Transition to dimensionless space ξ = r/r c and time τ = t/t c = c 0 •t/r c coordinates is performed [23] using the characteristic radius r c = (αE / κp 0 ) 1/3 that depends on the released energy E, where c 0 is the sound speed, α = 1.175 and heat capacity ratio κ = 1.4.Primitive variables are normalized with their respective values in the undisturbed gas.The solution thus becomes independent of the released energy and undisturbed gas state, therefore, numerical calculation can be performed once and scaling back into dimensional coordinates adapts it to the given conditions.
In the intermediate and weak shock range, only a numerical solution of the model is obtainable.We use an explicit discrete conservative numerical scheme [27]: where Δt and Δr denote time and spatial steps, respectively.The second order WAF finite volume explicit method is used for the calculation of the fluxes [27]: where the mean average state equals: W denotes the vector of primitive variables and c k the Courant's number.Intercell states are computed with an approximate HLLC Riemann solver [27].The Van Leer limiter ϕ is used to prevent numerical oscillations.Time integration is performed with the forward Euler method.An analytical solution of the strong shock theory serves as an initial condition for the numerical computation.Computation is performed over the dimensionless time interval 1 ≤ τ ≤ 35.The grid is equally spaced with 2 15 points for both independent variables.In the scope of this paper, only pressure is used.According to the model, the pressure transient in a given point is a function of two parameters -the distance from the source r s and the released energy E. Calculated pressure transients, normalized by the shock pressure amplitude, are shown in Fig. 3 for dimensionless distances ξ s : 2 (dashed), 8, 14, 20, 26 and 32 on a time scale retarded by the time of flight of the acoustic wave-front to facilitate waveform comparison.Fig. 3 illustrates a common characteristic of spherical shock waves: the of the compression phase increases with distance ξ s due to supersonic shock wave-front propagation while the duration of the rarefaction phase remains constant.

EXPERIMENT
We use water as the tissue phantom for the purpose of stable and repeatable experimental validation.At the particular Erbium laser wavelength, water trapped or bonded to the tissue plays a key role in tissue ablation.Explosive expansion of laser-heated water generates strong shock waves that propagate in air above the ablated surface, followed by material ejection.The key difference between the ablation of water and the tissue is in the form of the surface on which the resulting shock waves are generated.Water and some soft tissues form a quasi-ideal half-space whereas hard tissues usually form complicated geometry that varies from pulse to pulse during the ablation and affects the spread of the shock wave.The results obtained by ablating the water surface are relevant for the laser surgery in which the shallow holes (craters) are prepared.For laser drilling of deep holes, an appropriate shock-wave propagation model has not yet been developed.In this case the accuracy of the presented method is questionable.

Fig. 4. Experimental set-up: shock wave (SW), piezoelectric sensor (PE), photo-diode (PD), oscilloscope (OSC, pulse generator (PG), personal computer (PC)
The experimental system used to validate the above model is shown in Fig. 4. A free-running Er:YAG laser (Fidelis Plus III, Fotona) is used to irradiate the water surface as it forms a quasi-ideal half-space.The focal point of the laser exit optics is located on the water surface with a spot diameter of 0.9 mm.A signal from a pulse generator triggers the laser system, also setting the pulse duration of the laser flash lamp (45 μs).The supply voltage for laser pumping has been set to 650 V. Resulting laser pulses are short (≈2 μs), causing generation of a single shock wave rather than several that are typical for longer pulses [22].The pulse energy was measured (Nova II, Ophir), where the mean value was 3.14 mJ with a std.deviation of 0.16 mJ.A piezoelectric sensor (CA-1135, Dynasen) with a PZT-5A crystal disc of 1 mm in diameter was translated in parallel and perpendicular to the water surface using two linear stages with micrometer screws.
Responses were measured at three different horizontal distances (7.25, 10.25 and 13.25 mm) and six vertical ones (10, 12, 14, 16, 18 and 20 mm).Six repetitions on each configuration were conducted.The measured signals were sampled using a digital oscilloscope (Agilent DSO6034A) with 300 MHz bandwidth and 12 bit digitization.Signal acquisition was triggered by a signal from an InAs photodiode (J12, Teledyne), mounted behind the back laser mirror.The signals were sent to a PC where they were saved.Room temperature and pressure were measured: T 0 = 298 K, p 0 = 996 mbar, using standard meteorological equipment.

RESULTS AND DISCUSSION
We validate the described model by comparing the theoretical and measured waveforms.To enable this, the released energy is determined for every measured signal by measuring the duration of the compression phase of the measured signal.
Transforming the numerical solution of the point explosion to the dimensional coordinates and using the model of the sensor, we determine a function (Eq.11) that describes the duration of the compression waveform phase t p as a function of characteristic radius r c , given the influencing parameters: relative position of the sensor to ablation spot (determined by s 0 and h 0 ) and the sound speed c 0 which is estimated from the measured room temperature.The function is formulated as a spline with 41 data points (knots) along the characteristic radius r c for values between 0.8 and 2.8 mm, whilst other parameters are held constant.At each data point a theoretical waveform (signal) is found first, from which time t p is then determined.The same procedure is repeated for other sensor positions.
Using measured data for t p we numerically solve the above equation for r c .Released energy for the spherical blast wave E is then determined from the definition of the characteristic radius r c .Duration t p is determined by normalizing the signal with its peak value and searching for the time interval where the normalized amplitude exceeds the threshold value of 0.1 (10%).Signals are normalized because of simplifications of the theoretical model, where absolute signal values are unknown, Eq. ( 2), as the sensor is not calibrated.
The results are shown in Fig. 5. E h is the released energy for the half-space in which the shock wave forms a hemisphere (E h = E/2).No systematic dependency of E h on sensor positions (s 0 , h 0 , defined as shown in Figs. 2 and 4) is observed.It is to be noted that only sensor position (s 0 , h 0 ) varies in this experiment, while all other parameters (especially laser energy) remain the same.Pulse-topulse variations of E h within a sequence of repeated measurements, indicated by the error bars in Fig. 5, are mainly due to pulse-to-pulse variations of the laser pulse energy.We take the observation that E h does not exhibit systematic dependency on sensor position as evidence that the described model correctly describes the key features of the signals and the set-up.It is of interest to note the obtained E h value: its estimated mean value at 0.16 mJ (and std.dev. of 0.02 mJ) implies that about 5% of the incident laser pulse is converted into the energy of the shock wave.
Using estimated released energies E h , normalized theoretical signals are determined and compared to corresponding measured signal waveforms.Fig. 6 shows two limiting examples, characterized by high (53°) and low (20°) angles of incidence between the shock wave-front and the sensor aperture normal.
In both examples, good agreement of the positive (compressive) phase duration is observed between the theoretical and measured waveforms.This is a)    As a further step in model validation, we compare the theoretical and measured times of flight (TOF) of the shock wave-front.The maximum observed relative difference between a measured TOF and its theoretical counterpart is about 1%.Expressing this in terms of distance using the theoretical shock speed, we get a maximum distance uncertainty of about 0.24 mm.We attribute these deviations mainly to the uncertainties of the parameters in Eq. ( 11) and the measured durations t p that are used by the method for the estimation of the released energy E h .
By comparing the theoretical shock waveforms that take into account the sensor model to the respective ones that do not, we find that the former ones systematically exhibit longer duration of the compressive phase.The relative difference between the two increases with the angle of incidence from 5% at the 20° angle to 25% at the 53° angle.This is a result of the prolonged transit time of the shock wave over the finite sensor surface in the case of a larger angle of incidence.This observation presents strong evidence that the finite dimensions of the sensor in given conditions need to be taken into account.
The described approach can be employed to calibrate the sensor.Peak pressures are found from the estimated energies E h using the point explosion model.Then, sensor sensitivity is calculated by dividing the peak voltage of a measured signal by the peak pressure of the corresponding theoretical signal.Analysis of the estimated sensor sensitivities shows that it varies significantly with sensor orientation (s 0 , h 0 ): the estimated pressure peaks in Fig. 6 are 17 and 12 mbar for the upper and lower trace, respectively, and the corresponding sensor sensitivities are 2.4 and 5.7 μV/Pa, respectively.The observed variations of the estimated sensitivities at a fixed sensor position are less than 10%.

CONCLUSION
We describe a method and a set-up that opens the possibilities for on-line process monitoring in Er:YAG laser ablation.We employ a piezoelectric shock sensor to detect shock waves generated in the air above the irradiated surface.Using a comprehensive sensor model and the Taylor-Sedov model of shock wave propagation, we demonstrate excellent agreement between the measured and theoretical waveforms.On this basis we propose that the released shock energy is a process characteristic that is essentially independent of the position and orientation of the sensor.We also demonstrate that this method allows peak shock wave pressure estimation and calibration of the sensor.

ACKNOWLEDGEMENT
The operational part of this research was financed by the European Union and the European Social Fund.

Fig. 2 .
Fig. 2. Geometrical relations used in derivation of the sensor model -placement of the sensor relative to the ablation spot (P): side view (above), top view (below)

Fig. 3 .
Fig. 3. Normalized theoretical waveforms retarded by the acoustic wave-front propagation time