Emission tomography of laser induced plasmas with large acceptance angle apertures
aa r X i v : . [ phy s i c s . p l a s m - ph ] J u l Emission tomography of laser induced plasmas with largeacceptance angle apertures
S. V. Shabanov , and I. B. Gornushkin Department of Mathematics, University of Florida, Gainesville, FL 32611, USA BAM Federal Institute for Materials Research and Testing, Richard-Willst¨atter-Strasse11, 12489 Berlin, Germany
Abstract
It is proposed to use apertures with large acceptance angles to reduce the inte-gration time when studying the emissivity of laser induced plasmas by means of theAbel inversion method. The spatial resolution lost due to contributions of angledlines of sight to the intensity data collected along the plasma plume diameter isrestored by a special numerical data processing. The procedure is meant for thelaser induced plasma diagnostics and tomography when the integration time neededto achieve a reasonable signal to noise ratio exceeds a characteristic time scale ofthe plasma state variations which is short especially at early stages of the plasmaevolution.
Laser-induced breakdown spectroscopy (LIBS) is widely used for quantitative chemicalanalysis of an elemental content of various materials [1, 2]. A typical experimental setup(see, e.g., [3, 4, 5]) is depicted in Fig. 1. A plasma plume is created by a laser ablationat the distance of two focal lengths from a lens. The lens projects the plasma plumeimage onto the image plane where a detecting device is placed (e.g., a spectrometerslit). The plasma plume expands, cools, and radiates. All photons propagating alonga particular line of sight parallel to a vector e through a point at the distance y fromthe plasma plume symmetry axis can be collected at the point (0 , − y, f ), e.g., by aspectrometer slit. The intensity I Ω ( y ) measured at the point (0 , − y, f ) depends on thesolid angle Ω in which the light is collected. If Ω is small, then the intensity I Ω ( y )can be approximately viewed as the intensity I ( y ) in (nearly) parallel lines of sight. Inthis approximation, the measured intensity is related to the plasma emissivity by theAbel integral equation that can be solved to recover the emissivity [6]. For this reason,spectrometers used in these measurements must have a high f − number in order to restrictthe solid angle of the acceptance cone. The latter implies, in particular, a small amplitudeof the measured signal per unit time. In turn, a long integration time is required toachieve a reasonable signal to noise ratio. If the plasma plume remained stationary during1he integration time, the reconstructed emissivity (as well as the plasma temperaturedistribution and local densities of electrons, ions, and atoms) would have correspondedto the true (instantaneous) plasma parameters. It appears, however, the plasma plumecannot be viewed as stationary during a typical integration time of 1 µs or higher. Thisis especially true for early stages of the plasma plume evolution [7, 8].Here the following possibilities are discussed to reduce the integration time. If theplasma plume has the spherical symmetry (or, at least it can be approximated as such),then it is shown that any device with an arbitrary large acceptance cone can be used tocollect the data. In particular, spectrometers with low f − numbers can be used. The lossof the spatial resolution due to contributions of angled lines of sight into the data function I Ω ( y ) can be eliminated by a proposed numerical data processing such that the intensity I ( y ) in parallel lines of sight is accurately recovered and can be used in the subsequentAbel inversion to obtain the plasma emissivity. The data processing is particularly simpleif the device acceptance cone is slightly restricted by a large aperture of a special shape(see Section 3) placed in the lens focal plane (as depicted in Fig. 1). But with minormodifications the reconstruction algorithm is shown to work for any acceptance cone ofthe measuring device (or any aperture).If the plasma plume has the axial symmetry , then a measuring device may also havean arbitrary large acceptance cone, but the lines of sights that are not perpendicular tothe plasma plume symmetry axis must be blocked by an aperture in the focal plane whosedimension along the symmetry axis is sufficiently small. In other words, only lines of sightthat lie within a thin slice of the plasma plume (normal to the symmetry axis) contributeto the measured intensity. This is similar to the conventional experimental setup (e.g., aspectrometer with a high f − number). The difference is that photons propagating alongall angled lines of sight within the plasma plume slice can be collected to amplify themeasured signal. A simple numerical data processing algorithm is proposed to accuratelyrecover the intensity in parallel lines of sight.So, a rational behind the use of larger apertures is a reduction of the integration time.This is especially important when studying such a dynamical object as laser-inducedplasmas. The plasma diagnostics and tomography based on the Abel inversion becomesinaccurate when the integration times exceeds a characteristic time scale of the plasmastate variations which is short especially at early stages of the plasma evolution (e.g., afew ns after the ablation process is over). The method proposed below would allow fora sufficient signal to noise ratio with integration times in the range of 10’s to 100’s ns atthe price of loosing a spectral resolution. It also allows for the use of spectrometers withlow f − numbers, which reduces the integration time, but not as much as in the case whenthe spectral resolution is not needed. 2 u,v)(cid:13) u(cid:13)y(cid:13)y(cid:13)0(cid:13) 2f(cid:13) 4f(cid:13)A(cid:13) Image(cid:13)plane(cid:13) z(cid:13)Focal(cid:13)plane(cid:13)Plasma(cid:13) e(cid:13) e(cid:13) Lens(cid:13) -y(cid:13)
Figure 1:
A typical experimental setup for LIBS plasma measurements. A plasma plume iscreated by a laser ablation at the distance of two focal length from a lens. The light intensityis measured in the image plane. An aperture is placed in the lens focal plane. It restricts theacceptance cone at each point where the intensity is measured (provided the acceptance cone ofthe light collecting device is wider than that of the aperture). Photons propagating along anyparticular line of sight parallel to e = e θ,ϕ that lies within the solid angle of the acceptance coneand goes though the point (0 , y,
0) in the plasma plume are collected at the point (0 , − y, f ).The measurements give the intensity I Ω ( y ) as a function of y . For large acceptance cones, I Ω ( y )cannot be viewed as the intensity in parallel lines of sight and, hence, cannot be directly usedin the Abel inversion to reconstruct the plasma emissivity. Suppose that a plasma plume is created on the optical axis of a lens that is at a distanceof two focal lengths from the plume (see Fig. 1). Then an image of the plasma plume iscreated in the plane normal to the optical axis that is at the distance of two focal lengthsbut on the other side of the lens. Let the plasma emissivity be given as a function ofthe position vector r , ε = ε ( r ). The question to be studied is the amount of plasmaradiation that can be collected at a point in the image plane. The coordinate systemis set so that the optical axis coincides with the z − axis. Consider an infinitesimal areaelement dσ = dxdz at the point (0 , y,
0) (the origin is inside the plasma plume) normalto the z − axis. Let e θ,ϕ = (cos ϕ sin θ, sin ϕ sin θ, cos θ ) be the unit vector that definesthe direction of a line of sight through the point (0 , y, θ and ϕ is transparent. The line has the angle θ with the optical axis. Variationsof ϕ correspond to lines that lie on the cone of the angle θ whose the axis passes through(0 , y,
0) and is parallel to the optical axis. By letting θ and ϕ range over the intervals[0 , π/
2] and [0 , π ], respectively, all lines of sight through the point (0 , y,
0) in the plasmaplume are parameterized. Without any restriction on the lens diameter, all such linesintersect at the point (0 , − y, f ) in the image plane, where f is the lens focal length.A lens has a finite size and yet it is subject to optical aberration effects. For thisreason a diaphragm should be placed in the lens focal plane (see Fig. 1). The diaphragmaperture restricts the ”emission” solid angle within which the light can be collected at3he point (0 , − y, f ). As noted before, spectrometers with high f − numbers are used tocollect the light (the spectrometer slit is placed in the image plane parallel to the y − axis). For a high f − number, the spectrometer acceptance solid angle is smaller than the”emission” solid angle. So, not all the plasma radiation that arrives at the observationpoint is collected. If the acceptance solid angle is small, then the measured intensity isproportional to the intensity in (nearly) parallel lines of sight going through the plasmaplume diameter (see below) that is used in the Abel inversion method.Let us put aside the question about how exactly the light is collected at the point(0 , − y, f ) and assume that all the light coming through the aperture can somehow becollected. Alternatively, one can think of any measuring device that has a large acceptancecone. As is clear from Fig. 1, any acceptance cone is equivalent to some aperture placedin the lens focal plane. So, measuring devices and apertures restricting the ”emission”cone can be discussed on the same footing.The energy collected per unit time, unit frequency, and unit area that was emittedinto a solid angle d Ω in the direction of e θ,ϕ reads [9] I Ω ( y ) = Z Ω d Ω cos θ Z ∞−∞ ds ε (cid:16) r θ,ϕ ( s, y ) (cid:17) ≡ Z Ω d Ω cos θ I ( θ, ϕ, y ) , (1)where the dependence on the frequency and time is not explicitly shown (it is not relevantfor what follows), r θ,ϕ ( s, y ) = (0 , y,
0) + e θ,ϕ s is the parametric equation of the line ofsight through the point (0 , y,
0) and parallel to the unit vector e θ,ϕ , s is the arclengthparameter, d Ω = dθdϕ sin θ , and the range for the angles θ and ϕ is determined by theshape of the aperture in the focal plane. The function I ( θ, ϕ, y ) is the light intensity alongthe line of sight parallel to e θ,ϕ through the point (0 , y,
0) in the plasma plume. To findthe dependence of I Ω on the aperture geometry, the integration over the solid angle istransformed to a double integral over the aperture by a change of variables. Let ( u, v ) berectangular coordinates in the focal plane such that the origin lies on the lens optical axisand the u − axis is parallel to the y − axis (see Fig. 1). If dA = dudv is the area element atthe aperture point ( u, v ), then d Ω = dA cos θf + u + v and the relation between the angles and rectangular coordinates in the focal plane is givenby v = f tan θ cos ϕ , u = f tan θ sin ϕ . (2)Equation (1) can be written in the form I Ω ( y ) = Z A dA f ( f + u + v ) I ( θ, ϕ, y ) , (3)where the angles θ and ϕ in the arguments of I are expressed via u and v by invertingrelations (2). The integration in (3) is carried out over the aperture area denoted by A .4onsider first the case of a spherically symmetric plasma plume (the case of the axialsymmetry is discussed in Section 4), i.e., ε ( r ) = ε ( r ) where r = | r | . The distance r θ,ϕ ( y )between the line r = r θ,ϕ ( s, y ) and the origin (the plasma plume center) is given by themagnitude of the vector that is the cross product (0 , y, × e θ,ϕ . A simple calculationshows that r θ,ϕ ( y ) = | (0 , y, × e θ,ϕ | = y (1 − sin θ sin ϕ ) / for y ≥
0. If I (0 , , y ) = I ( y )is the intensity along the line of sight parallel to the optical axis, then it follows from theassumption of the spherical symmetry of the plasma plume that I ( θ, ϕ, y ) = I ( r θ,ϕ ( y )) , r θ,ϕ ( y ) = y q − sin θ sin ϕ (4)because the intensity along any line of sight through the point (0 , y,
0) depends only onthe distance of that line to the plasma plume center. By making use of relations (2) toexpress r θ,ϕ as a function of u and v , Eq. (3) can be written in the form I Ω ( y ) = Z A dA f ( f + u + v ) I (cid:16) yµ ( u, v ) (cid:17) , µ ( u, v ) = s f + v f + u + v . (5)If A is an infinitesimal aperture of the area ∆ A centered on the optical axis, then I Ω ( y ) ≈ ∆Ω I ( y ) where ∆Ω = ∆ A/f . In this case, the data function can be used as the intensityalong parallel lines of sight I ( y ) ≈ I Ω ( y ) / ∆Ω that is related to the plasma emissivity ε ( r )by the the Abel integral equation I ( y ) = Z ∞−∞ ds ε (cid:16)q y + s (cid:17) = 2 Z ∞ y dr r ε ( r ) √ r − y . (6)It can be solved for ε ( r ). The procedure is known as the Abel inversion [6]: ε ( r ) = − π Z ∞ r dy I ′ ( y ) √ y − r = Z ∞ dk kJ ( kr ) ˜ I ( k ) , ˜ I ( k ) = Z ∞−∞ dye iky I ( y ) , where ˜ I ( k ) is the Fourier transform of I ( y ) and J ν is the Bessel function of index ν = 0.Substituting (6) into the right side of Eq. (5), an analog of the Abel integral equation(6) is obtained for an aperture with a finite acceptance angle. If Eq. (5) can be solvedfor I ( y ) for a given I Ω ( y ), then a solution to the generalized Abel equation is found bymeans of the Abel inversion for I ( y ). In Sections 5 a numerical algorithm is proposed forsolving Eq. (5) in the case of an aperture of a special shape. The case of an aperture ofa general shape is presented in Section 6. A numerical example of is given in Section 7. It turns out that there is a particular shape of the aperture for which Eq. (5) is substan-tially simplified and becomes convenient for a qualitative analysis of the relation betweenthe intensity in parallel lines of sight I ( y ) and the data function I Ω ( y ) for large apertures.5 (cid:13)v(cid:13) u(cid:13)v(cid:13) u(cid:13) a(cid:13) A(cid:13) nk(cid:13) -(cid:13)
A(cid:13) nk(cid:13) +(cid:13) a(cid:13) u(cid:13) a(cid:13) u'(cid:13) a(cid:13)
Figure 2:
Left panel:
A hyperbolic aperture that allows for a factorization in the sum overangled lines of sight. It is bounded by the lines v = ± ∆ a and hyperbolas u = ± tan θ a p v + f .Here u a = f tan θ a = f p − µ a /µ a and u ′ a = tan θ a p f + ∆ a . For ∆ a ∼ f , the aperture cannotbe approximated by a rectangle. A rectangular approximation is accurate if ∆ a ≪ f . It canbe used for studies of axially symmetric plasmas as explained in the text. Right panel:
Anillustration to the reconstruction algorithm for a general aperture. Curves 1, 2, and 3 are thelevel curves (hyperbolas) µ ( u, v ) = µ − nk , µ ( u, v ) = µ + nk , and µ ( u, v ) = µ − nk +1 , respectively, thatare used to partition the region occupied by the aperture. In each partition region A ± nk , thefunction I n ( µ ) is approximated by a linear function as is explained in the text. Here µ a is theminimal value of the function µ ( u, v ). Let A be symmetric relative to the reflection u → − u and let A + be the half of A that lies in the half-plane u ≥
0. Then R A = 2 R A + in the right side of Eq. (5). Considerthe change of variables in the double integral (5): ( u, v ) → ( µ, v ) where the new variable µ = µ ( u, v ) is defined by the second relation in (5). Then dA = dudv = dvdµJ where theJacobian is J = ∂u/∂µ = ( ∂µ/∂u ) − . The coordinate line µ = const is a hyperbola inthe ( u, v ) − plane: u = √ − µ µ q v + f . (7)It intersects the u − axis at u = f √ − µ /µ and has slant asymptotes u = | v |√ − µ /µ at large | v | . In particular, the line u = 0 (or the v − axis) is the limiting hyperbola as µ →
1. Note that the u − intercept of the hyperbola (7) tends to zero, while the slope ofthe slant asymptotes becomes infinite, i.e., the hyperbola tends to the vertical straightline through u = 0. Let A + be bounded by horizontal lines v = ± ∆ a , the vertical line u = 0, and by the hyperbola (7) with µ = µ a < − ∆ a ≤ v ≤ ∆ a and µ a ≤ µ ≤
1. A straightforward calculation shows that the doubleintegral (5) is factorized in the new variables thanks to the ”hyperbolic” shape of theaperture: I Ω ( y ) = 4 Z ∆ a dv √ f + v Z µ a dµ µ √ − µ I ( µy ) = 4∆ a q f + ∆ a Z µ a dµ µ √ − µ I ( µy ) . (8)6ut µ = cos θ so that u = f tan θ a ≡ u a defines the u − intercept of the hyperbola (7) with µ = µ a that bounds the aperture area A + . Put also tan ϕ a = ∆ a /f . Then I Ω ( y ) = 4 sin ϕ a Z θ a dθ cos θ I ( y cos θ ) . (9)This equation is simple enough to estimate the signal amplification due to an aperturewith a large acceptance angle and compare it with signal measurements in reported LIBSplasma tomography experiments [4, 5].In the conventional LIBS plasma tomography experiments, the slit of a spectrometer isplaced into the image plane parallel to y − axis. However, not all the light that enters intothe spectrometer slit at the point (0 , − y, f ) is collected but rather its only portion thatilluminates a spectrometer mirror (that further directs the light onto the spectral grating).If the aperture in the focal plane is not too small or not present at all, this portion isdetermined only by the f − number of the spectrometer so that the acceptance solid anglein Eq. (1) is defined by 0 ≤ ϕ ≤ π and 0 ≤ θ ≤ θ s where tan θ s = 1 / (2 f s ) and f s isthe spectrometer f − number. Spectrometers in the reported LIBS plasma tomographyexperiments typically have relatively large f − numbers, f s ∼
5. So, the acceptance angleis small, tan θ s ≈ θ s ∼ − . It is then sufficient to calculate the integral (1) in theleading order of θ s , while neglecting terms of order θ s . In this approximation I ( θ, ϕ, y ) = I ( y ) + O ( θ ) (cf. (4)) and the evaluation of the integral (1) yields: I Ω ( y ) = πθ s I ( y ) + O ( θ s ) . The smallness of the factor πθ s in the rate of the energy flux per unit time is exactlythe reason to have a long integration time to get a sufficient signal to noise ratio, whendetermining I ( y ) used in the Abel inversion. Suppose now that all the light that comesto the point (0 , − y, f ) through the hyperbolic aperture is collected. Assuming, for sakeof simplicity, I ( y cos θ ) to be a constant I , the evaluation of the integral in (9) yields: I Ω ≈ ϕ a (2 θ a − sin(2 θ a )) I .
So the signal amplification (as compared to the spectrometer case) can be estimated bya factor of 4 sin ϕ a (2 θ a − sin(2 θ a )) / ( πθ s ) ∼
50 if θ s ∼ − and ϕ a ∼ θ a ∼ π/ a ∼ u a ∼ f ).To increase the angles ϕ a and θ a , spectrometers with small f − numbers can be used.An effective aperture corresponding to this case is a circular disk of the radius f / (2 f s ).If the hyperbolic aperture is fit into this disk, then Eq. (9) applies directly, and I ( y ) isrecovered by the numerical data processing proposed in Section 5. Otherwise, a circularaperture should be treated by the algorithm of Section 6.If the spectral resolution is not needed, then a gated ICCD camera can be placed inthe image plane. The pixel row of the camera is set along the plasma image diameter. Thelight collected by each pixel provides the spatially resolved data I Ω ( y ). If the aperture has7he hyperbolic shape, Eq. (9) applies directly to recover I ( y ). Otherwise, the algorithmof Section 6 must be used. To study the plasma emissivity in a narrow frequency windowof interest, a narrow band (or interference) filter can be positioned in front of the detector(the pixel array). This method can be particularly useful for diagnostics of laser inducedplasmas shortly after the ablation process ( ∼ ns ). At this stage the plasma exhibitsa fast dynamics and emits mostly a continuum spectrum radiation. For axially symmetric plasma plumes, photons propagating along lines of sight that arenot perpendicular to the plume symmetry axis must be prevented from being collected bya detector so that only a plasma plume slice normal to the symmetry axis is visible to ameasuring device. For spectrometers with high f − numbers, this is achieved automaticallydue to their narrow acceptance cones. But a narrow acceptance cone also eliminates mostof lines of sight in the visible slice itself that are angled to the optical axis. If R is aplasma plume radius, then the lines of sight through (0 , y,
0) in the acceptance cone ofthe angle θ s can go through plasma points separated at most by the distance 2 R tan θ s .This determines the spatial resolution along the y − axis (as well as along the symmetryaxis). For θ s ∼ − , the resolution is such that the data can be collected at about 10positions across the plasma plume diameter without a significant overlap of the plasmaregions from which the light is collected. So, when spectrometers with lower f − numbersor pixel rows of an ICCD camera (as suggested above) are used to increase the signal, thespatial resolution become even worse along both the directions, the y − and symmetryaxes.The following compromise can be made. The spatial resolution in the direction parallelto the symmetry axis is achieved by placing a slit-like aperture extended along the y − axisin the focal plane. If a slit of height 2∆ a is positioned in the focal plane ( − ∆ a ≤ v ≤ ∆ a in the above notations), then only photons traveling along the lines of sight that are inthe plasma ”slice” of height 2 R tan ϕ a = 2 R ∆ a /f can reach the image plane. The shapeof a sufficiently narrow rectangular aperture can well be approximated by the hyperbolicaperture with small ∆ a because the hyperbola (7) has the vertical tangent line at theintercept point u = u a (see the caption of Fig. 2). Thus, Eq. (8) or (9) remains applicablein this case, but the angle θ a is determined by the measuring device f − number, i.e.,tan θ a = 1 / (2 f s ) (the approximation of a small θ s is invalid for small f s ), and ϕ a is smallenough to provide for the spatial resolution along the plasma plume symmetry axis.Two-dimensional detectors used in LIBS experiments [4, 5] in combination with aspectrometer typically have about a few hundred pixels for the spatial resolution along thespectrometer slit, while only 10 to 20 spatially resolved data points are taken (each datapoint corresponds to an average intensity over a cluster of pixels). The algorithm providedbelow allows for an accurate reconstruction of the spatial resolution of just a single pixelsize. So, it may not be deprived of sense to apply the proposed algorithm to improve the8patial resolution even in spectrometers with higher f − numbers. This, however, dependsvery much on the gradients of thermodynamic parameters of the observed plasma. Forhigh gradients, this is definitely the case. Let the light be collected by a linear array of pixel detectors placed along the y − axis inthe image plane. Each pixel has a size ∆. The pixel centers are located at y = y n = n ∆, n = 0 , , ..., N (by the symmetry of the observed plasma, it is sufficient to reconstruct I ( y ) (the intensity in parallel lines of sight) only for non-negative y , I ( − y ) = I ( y )). Inwhat follows the factor 4 sin ϕ a in (9) will be omitted. The energy of all photons collectedby a pixel n per unit time and unit frequency is given by I Ω n = Z y n +1 / y n − / dyI Ω ( y ) = Z µ a dµ µ √ − µ Z y n +1 / y n − / dyI ( µy ) = Z µ a dµ µ √ − µ I n ( µ ) (10) I n ( µ ) = Z µy n +1 / µy n − / du I ( u ) , (11)where y n ± / = ∆( n ± /
2) are the boundary points of the n th pixel, and the change ofthe integration variable u = µy has been carried out. If only the photons propagatingalong parallel lines of sight were registered by the n th pixel, then the data would haveconsisted of values I n = Z y n +1 / y n − / du I ( u ) = I n (1) . (12)Exactly these values must be used in the Abel inversion to reconstruct the emissivity.The idea is to find a suitable interpolation for the functions I n ( µ ) using the ”sampling”points I n . Then the integral with respect to µ in (10) can be evaluated thus turning Eq.(10) into a system of algebraic equations for I n that can be solved for a given data set I Ω n .Let P n ( µ ) denote the interval [ µy n − / , µy n +1 / ] so that P n (1) = P n is the intervaloccupied by the n th pixel. Depending on the values of µ and n , the interval P n ( µ ) mayhave the following positions relative to the intervals occupied by physical pixels: P n ( µ ) overlaps with both P n − k +1 and P n − k , or P n ( µ ) is contained in P n − k (13)for k = 1 , , ..., N n where N n = [(1 − µ a ) n ] + 1 and [ x ] denotes the largest integer thatdoes not exceed x . The integer N n is defined by the condition y n − N n < µ a y n ≤ y n − N n +1 .For a fixed value of n it is not difficult to find the ranges of µ in both the cases (13): µ + nk ≤ µ ≤ µ − nk , or µ − nk +1 ≤ µ ≤ µ + nk , (14) µ ± nk = n − k + 1 / n ± / , k = 1 , , ..., N n , µ − nN n +1 = µ a . (15)For example, for k = 1, the first case in (13) holds as long as µy n +1 / remains in P n , i.e., µy n +1 / ≥ y n − / or µ ≥ µ + n . As µ becomes less than µ + n , the whole interval P n ( µ )9emains in P n − until the left boundary of P n ( µ ) crosses the left boundary of P n − , i.e., µy n − / ≥ y n − / or µ ≥ µ − n , and so on. A special consideration is required for the case k = N n because, depending on the value of µ a , either both the cases (13), or just the firstone, or none of them may occur. This issue will be clarified shortly. For time being, thevalues µ ± nk are defined by (15), i.e., µ a is assumed not to exceed µ nN + n (both the cases (13)are possible for k = N n ).Consider two integrals R ba dx f ( x ) and R b ′ a ′ dx f ( x ) for some f ( x ) where a ≤ a ′ < b ′ ≤ b .How can the second integral be approximated by the first one? The following approxima-tion is proposed: Z b ′ a ′ dx f ( x ) ≈ b ′ − a ′ b − a Z ba dx f ( x ) . (16)The accuracy of this approximation is discussed in Appendix. For a small interval length b − a and smooth f , this is a good approximation as is shown in Appendix. If µ is suchthat the interval P n ( µ ) is contained in P n − k , then the approximation (16) applies directly: I n ( µ ) = µ ∆∆ I n − k = µ I n − k , µ − nk +1 ≤ µ ≤ µ + nk . (17)If µ is such that P n ( µ ) overlaps with P n − k and P n − k +1 , then I n ( µ ) is approximated by alinear combination of I n − k +1 and I n − k with the coefficients being the lengths of fractionsof P n ( µ ) in P n − k +1 and P n − k , respectively, in units of the pixel size ∆: I n ( µ ) = (cid:16) µa n − a n − k (cid:17) I n − k +1 + (cid:16) µ (1 − a n ) + a n − k (cid:17) I n − k , µ + nk ≤ µ ≤ µ − nk , (18)where a n = n + 1 /
2. The integral over the interval [ µ a ,
1] in (10) is split into the sum ofintegrals over the intervals in which the approximations (17) and (18) hold I Ω n = N n X k =1 Z µ − nk µ + nk + Z µ + nk µ − nk +1 ! dµ µ √ − µ I n ( µ ) . (19)The integrals are easily evaluated and the desired system of linear equations for I n isobtained. There is a simple recurrence algorithm to solve it. Before describing it, theissue about the integration limits in the term k = N n must be addressed.Consider the last three integrals in the sum (19): Z µ + nNn µ − nNn − + Z µ − nNn µ + nNn + Z µ + nNn µ a , (20)in accord with the definition (15). By examining the possible overlaps of P n ( µ ) with P n − N n and P n − N n +1 depending on the value of µ a , the following three possibilities arefound to occur. If µ a ≥ µ − nN n where the value of µ − nN n is computed by the rule (15),then the last two integrals in (20) cannot occur in the sum (19) and must be removed.If µ + nN n ≤ µ a < µ − nN n where the values of µ ± nN n are calculated by the rule (15), then thevery last integral in (20) cannot be present and must be removed from the sum (19).10nly if µ a < µ nN n all three integrals (20) are present in the sum (20). For algorithmimplementation purposes, it is convenient to retain the sum (19) as it is. The above threecases can be accounted for if, after calculating the values of µ ± nk by the rule (15), thevalues of µ ± nN n are redefined by the rule: µ − nN n → max( µ a , µ − nN n ) , µ + nN n → max( µ a , µ + nN n ) . It is straightforward to verify that the unwanted integrals in (20) (and, hence, in (19))always vanish as their upper and lower limits coincide.Define the functions F ( α, β ) = Z βα dµ µ √ − µ = √ − α − q − β ,F ( α, β ) = Z βα dµ µ √ − µ = 12 (cid:16) cos − ( α ) − cos − ( β ) (cid:17) + 12 (cid:16) α √ − α − β q − β (cid:17) . They are used to evaluate the integrals in the sum (19). For the central pixel ( n = 0),one finds I = I Ω0 F ( µ a , . (21)For n = 1 , , ..., N , Eq. (19) yields the recurrence relation I n = I Ω n − I Fn C n , (22)where C nk = F ( µ + nk , µ − nk ) a n − F ( µ + nk , µ − nk ) a n − k , k = 1 , , ..., N n ,I Fn = N n X k =2 C nk I n − k +1 + N n X k =1 h F ( µ − nk +1 , µ − nk ) − C nk i I n − k . Note that I Fn depends only on I , I ,..., I n − . So by executing the recurrence relation (22)with the initial condition (21), all the values I n are found. They can then be used todetermine the emissivity by a suitable Abel inversion algorithm. There is a vast literatureon the numerical Abel inversion (see, e.g., [10]), and, for this reason, it will not be discussedhere. In Section 7, it is demonstrated that for emissivity functions representative for LIBSplasmas, the proposed algorithm provides a very accurate reconstruction of I n . Thus,the accuracy in the emissivity reconstruction would be determined by the accuracy of aparticular Abel inversion algorithm. It is not difficult to extend the proposed algorithm to an aperture of a generic shape.Suppose that A is symmetric under the reflection u → − u and A + is the half of A for11hich u ≥
0. A generalization to the non-symmetric case is trivial. Integrating Eq. (5)over the interval occupied by the n th pixel, a generalization of (10) is obtained I Ω n = 2 Z A + dA f ( f + u + v ) I n ( µ ( u, v )) , (23)where the function µ ( u, v ) is defined by the second equation in (5). The area A + can bepartitioned into regions that lies between the level curves (hyperbolas) of the function µ ( u, v ): µ ( u, v ) = µ + nk , µ ( u, v ) = µ − nk , where µ ± nk are defined in (15) and µ a = min A + µ ( u, v ). The procedure is shown in theright panel of Fig. 2. Let A − nk and A + nk be the regions in which µ + nk ≤ µ ( u, v ) ≤ µ − nk and µ − nk +1 ≤ µ ( u, v ) ≤ µ + nk , respectively. Then a generalization of Eq. (19) to the case of ageneric aperture reads I Ω n = 2 N n X k =1 Z A − nk + Z A + nk ! dA f ( f + u + v ) I n ( µ ( u, v )) . (24)The approximations (17) and (18) are used for I n ( µ ) in the regions A + nk and A − nk , respec-tively to evaluate the double integrals in (24). Depending of the shape of A , it can bedone either analytically or numerically. Put F m ( B ) = 2 Z B dA f ( f + u + v ) µ m ( u, v ) , m = 0 , , for a planar region B . The initial condition (21) and the recurrence relation (22) become I = I Ω0 F ( A + ) , I n = I Ω n − I Fn C n , n = 1 , , ..., N , (25)where C nk = a n F ( A − nk ) − a n − k F ( A − nk ) , a n = n + 1 / ,I Fn = N n X k =2 C nk I n − k +1 + N n X k =1 (cid:16) F ( A nk ) − C nk (cid:17) I n − k , and A nk is the union of A + nk and A − nk . The reconstruction algorithm has been tested with synthetic data. The emissivity func-tion is taken in the form ε ( r ) = ε exp( − r / (2 σ ))( a − cos( br )) where the values of theparameters (arbitrary units) are ε = 1, σ = 3, a = 1 .
4, and b = 0 .
5. The emissivityprofile is typical for the LIBS plasmas observed in experiments [5]. The intensity I ( y ) in12
10 0 10−0.200.20.40.60.811.2 y −10 0 10−0.200.20.40.60.811.2 y
Figure 3:
Numerical results for a hyperbolic aperture with µ a = cos θ a , θ a = 51 . o . Left panel:
The solid (black) curve is the graph of the intensity I ( y ) in parallel lines of sight for the emissivity ε ( r ) described in the text. The dashed (blue) curve is the observed intensity I Ω ( y ) calculatedfrom I ( y ) as specified in Eq. (8). Both the curves are normalized on their maximal values. Theamplification of the signal is not shown explicitly. It is determined as explained in the text.Solid (red) dots are reconstructed values of I n from I Ω n . The number of pixels is 60 ( N = 31). Right panel:
The same as in the left panel but the number of pixels is 30 ( N = 16). A smallinaccuracy of the reconstructed values of I n is observed near the maxima of the intensity I ( y )which occurs through errors in the approximations (17) and (18). parallel lines of sight is calculated by means of Eq. (6). Then the integral intensity overall lines of sight within the acceptance solid angle defined by the hyperbolic aperture with µ a = cos θ a , θ a = 51 . o is found by Eq. (9) (the amplification factor 4 sin ϕ a is omittedas it is irrelevant for the reconstruction algorithm). The synthetic data are calculated bymeans of Eq. (10) when 60 and 30 pixels are placed across the plasma plume diameter.Then the reconstruction algorithm is applied to find the values I n . The results are shownin Fig. 3. The left and right panels represent the cases with 60 and 30 pixels, respectively.The solid (black) and dashed (blue) curves are the graphs of I ( y ) and I Ω ( y ), respectively.Both the curves are normalized to their maximal values so that the amplification factordue to the vertical size ∆ a of the hyperbolic aperture cannot be seen. This normalizationwould not make sense for a generic aperture because there is no factorization in the sumover angled lines of sight for apertures of a general shape (see Section 6). Yet, the factor-ization due to the hyperbolic aperture leads to a mild deviation of I Ω ( y ) from I ( y ) whenthey are normalized on their maximal values. The solid (red) dots show the reconstructedvalues of I n . The reconstruction is very accurate for the case of 60 pixels. When thepixel size increases, the approximations (17) and (18) becomes less accurate as proved inAppendix. A reconstruction error increases near the maxima of I ( y ) when the pixel sizeis doubled (the right panel of Fig. 3). For a typical pixel detector, the pixel size is about0 . mm . So, for a typical LIBS plasma plume diameter of 2 mm , there are about 100pixel collecting the data. 13 . Appendix Here the accuracy of the approximation (16) is assessed. For a continuous function f ( x )there exist points a ≤ x ∗ ≤ b and a ′ ≤ x ′∗ ≤ b ′ such that I = Z ba dxf ( x ) = ( b − a ) f ( x ∗ ) , I ′ = Z b ′ a ′ dxf ( x ) = ( b ′ − a ′ ) f ( x ′∗ ) . It follows from these relations that I ′ = ∆ ′ ∆ I + ∆ ′ (cid:16) f ( x ′∗ ) − f ( x ∗ ) (cid:17) , where ∆ = b − a and ∆ ′ = b ′ − a ′ . By the mean value theorem, the last term in thisequality can be transformed so that I ′ = ∆ ′ ∆ I + ∆ ′ ( x ′∗ − x ∗ ) f ′ ( ξ ) ≡ ∆ ′ ∆ I + γ for some ξ between x ′∗ and x ∗ . Therefore the error γ is bounded by | γ | ≤ M ∆ ′ ∆ , M = max [ a,b ] | f ′ ( x ) | . If ∆ →
0, while ∆ ′ / ∆ remains finite, the error decreases as ∆ . Acknowledgments
The authors acknowledge stimulating discussions with Prof. U. Panne (BAM) and D.Shelby (UF, Chemistry). S.V.S. is grateful to Prof. U. Panne for his continued supportand thanks Department IV of BAM for a kind hospitality extended to him during his visit.The work of I.B.G. is supported in part by the DFG-NSF grant GO 1848/1-1 (Germany)and NI 185/38-1 (USA).
References [1] I.D. Winefordner, I.B. Gornushkin, T. Correll, E. Gibb, B.W. Smith, N. Omenetto,J. Anal. At. Spectrom. 19, 1061 (2004)[2] C. Pasquini, J. Cortez, L.M.C. Silva, F.B. Gonzaga, J. Braz. Chem. Soc. 18, 463(2007)[3] J. W. Olesik and G. M. Hieftje, Anal. Chem. 57, 2049 (1985)[4] R. ´Alvarez, A. Roberto, M. C. Quintero, Spectrochim. Acta B 57, 1665 (2002)145] J. A. Aguilera, C. Arag´on, and J. Bengoechea, Appl. Optics 42, 5938 (2003)[6] A. C. Eckrbeth,
Laser diagnostic for combustion temperature and species , (Gordon &Breach, Amsterdam, 1996);H. R. Griem,
Principles of plasma spectroscopy (Cambridge University Press, Cam-bridge, 1997)[7] A. Casavola, G. Colonna, and M. Capitelli, Appl. Surf. Sci. 208-209, 559 (2003)[8] M. Capitelli, F. Capitelli, and A. Eletskii, Spectrochim. Acta B 56, 567 (2001)[9] S. Chandrasekhar,
Radiative transfer (Dover Publications, New York, 1960)[10] W. Lochte-Holtgreven, in:
Plasma diagnostics , ed. W. Lochte-Holtgreven (North-Holland Pub. Co., Amsterdam, 1968), p. 135;J.D. Algeo and M.B. Denton, Appl. Spectrosc. 35, 35 (1981);W. R. Wing and R. V. Neidigh, Am. J. Phys. 39, 760 (1971);R.N. Bracewell,
The Fourier transform and its applications (McGraw-Hill Book Co.,New York, 1965);D.R. Keefer, L. Smith, and S.J. Sudhasanan,