Multilayer target PIXE spectral simulation (X,X) secondary fluorescence correction algorithm

Abstract

Simulation of particle induced X-ray emission (PIXE) spectra is not a recent subject. Still, when samples are not homogeneous, problems emerge even in the simplest case of layered samples. If it is necessary to consider the presence of the same chemical element in more than one physically distinct layer the number of available simulation codes is very small. In addition, although X-ray emission spectra from PIXE experiments are much less prone to significant secondary fluorescence issues than their X-ray fluorescence spectrometry (XRF) counterpart, cases do emerge where secondary fluorescence calculations are necessary to ensure good PIXE spectral simulations, even if corrections are small. The case of secondary fluorescence induced by primary X-rays in thick homogeneous samples was solved long ago by various authors. In the case of non-homogenous targets, the problem becomes much more complex and, although also addressed long ago, a general solution cannot be found in the standard accessible literature on the PIXE technique. In the present work we revise a secondary fluorescence correction method presented in 1996 to handle homogeneous targets and extend it to be applicable to multilayered targets. Its implementation in the DT2 code allows simulation of PIXE spectra taking into account this type of matrix effect correction in complex multilayer targets. Fluorescence between different physical layers, the possibility of the presence of one chemical element in more than one layer, and the potential “illusional” presence of a chemical element in a given layer due to secondary fluorescence effects, when its real concentration in that layer is null, are dealt with. This is the first of what is intended to be a series of three papers. In this part I work, the model is presented for the case of secondary X-rays induced by primary X-rays produced by particle collisions. Applications and potentially demanding experimental conditions will be dealt with in part II, and the case of secondary X-rays induced by primary radiation from non-radiative transitions of fast electrons will be addressed in part III.

---

1 Introduction

Quantitative work on particle induced X-ray emission (PIXE)¹ can be performed using a very simple approximation if the targets are thin enough so that the ion beam particles impacting the target do not lose any significant amount of energy while crossing it and the characteristic X-rays of sample elements are well separated in the spectra. In this case, if standard samples have been previously analysed under the same conditions, integrating the characteristic X-ray peaks, or even just using their height, will provide quantitative data without the need for much complex processing.²

Still, in many cases the situation is not so simple. If the target is not thin enough, the target X-ray yield must be determined by integration of the yield function along the ion beam particle path in the target, and it can even happen that enhancement of X-ray emission relative to the yield expected from particle induced ionizations takes place. In “standard” cases, as mentioned by Folkmann in 1974,³ it is important to consider the fluorescence processes that result from the absorption of primary X-rays (the X-rays induced directly by particle collisions), in particular those cases that result from the absorption of the primary characteristic X-rays in the sample material. This absorption is named self-absorption, and the fluorescence processes are usually referred to as secondary fluorescence, which is probably the most important phenomenon leading to this enhancement.

Being quite significant when studying some types (e.g.: metal alloys) of thick targets (targets that are thick enough to completely stop the incident ion beam), the X-ray yield enhancement effect due to secondary fluorescence was addressed long ago by several authors and solved for the case of homogeneous thick samples. In the case of PIXE work, Reuter et al. in 1975,⁴ Ahlberg in 1977 (ref. 5) and Richter and Wätjen in 1981 (ref. 6) presented analytical solutions to the problem. Van Oystaeyen and Demortier in 1983 (ref. 7) developed a Monte Carlo method; Campbell et al. in 1989 (ref. 8) calculated the need for tertiary corrections and Ryan et al., at the beginning of the 1990s,^9,10 implemented calculation processes in GeoPIXE to deal with thin layers and inclusions in complex geological samples.

The secondary fluorescence effect in PIXE is similar to what is observed in X-ray fluorescence spectrometry (XRF), and therefore some of these methods resemble and reflect the 1960s work of Shiraiwa and Fujino,¹¹ even though the primary yield determination in the case of PIXE cannot be handled simply as an exponential term and must be obtained by numerical calculation, which complicates all further calculations.

At the beginning of the 1990s decade, the issue was revisited by myself while developing the first version of the DATTPIXE package.¹² After the first approach based on the work of Ahlberg,⁵ a variant was developed taking the model of Richter and Wätjen⁶ as a working base to define a function of depth term for the secondary fluorescence correction, which can be added to the primary X-ray yield prior to integration along the particle penetration path. This model, then named the “penetration function model”, as presented in 1996,¹³ was applicable for thick and half-thick targets and was implemented as such in the DATTPIXE package 1996 version.¹³

As mentioned above, PIXE samples are considered thin if it is possible to assume that the energy loss of incoming particles after crossing the target is negligible. In practice, in many cases, this energy loss is not negligible and the samples must be considered either half-thick, if the beam particles emerge from the target, or thick if they are completely stopped inside it.

If the samples are not homogeneous in depth the simplest case that can be considered is that of layered targets. These are targets that can be modeled as a set of physically distinct layers, each of them being a thin or half-thick target that is crossed by the particles of the beam, which may in the end be stopped in a thick substrate on top of which the layers are successively present. In this case, a more complex situation is faced, both for yield calculation and even more for cases where the secondary fluorescence effect must be accounted for.

In the case of XRF, the handling of secondary fluorescence effects in layered targets has been described in detail by De Boer.¹⁴ In this case, since the primary and secondary excitation processes are identical, major correction terms may be expected in several cases since the ionization cross-section of the radiation inducing secondary fluorescence is higher than the corresponding ionization cross-section of the incident X-ray beam.

This is not the case in PIXE, since the particle collision ionization cross-sections of matrix atoms are, in most (if not all) of the cases, higher or even much higher than the ionization cross-sections of matrix atoms by the primary X-rays produced after the particle collisions.

In many cases, in PIXE experiments, secondary fluorescence enhancement effects in layered targets can, therefore, be neglected since it is reasonably possible to assume that any possible correction is very small. Nevertheless, since the PIXE technique is becoming increasingly used to study layered targets, frequently using a Total-IBA¹⁵ approach, complex problems are starting to emerge and secondary fluorescence calculations in layered targets can no longer be disregarded, even if just to ensure that they are small.

Although, also for PIXE, the problem of secondary fluorescence effects in non-homogeneous samples has been addressed since the beginning of the 1990s,^9,10 still, a systematic and detailed description of the general PIXE case of layered targets, similar to De Boer's work for XRF, could not be found in the standard accessible literature, even though it is mentioned in Ryan et al.’s 1990s papers as “in preparation”.

Besides this difficulty in finding calculation details on the 1990s work on the subject, the present paper focuses on PIXE spectral simulation, while previous work has so far focused on calculating changes that must be taken into account to fit spectral details. In fact, although the two goals share a significant fraction of problems, not all of them are exactly the same and the best solutions for one and other issues are also not fully coincident.

In this work, we revisit the secondary fluorescence correction penetration function model published in 1996 (ref. 13) for homogeneous thick and half-thick targets and extend it to include layered targets.

No limitation is set on the presence of elements in layers, meaning that elements may be repeated in different physical layers and/or generate secondary X-rays due to primary radiation originating in layers where they are not physically present; in such cases the “illusion” of an element being present where the primary radiation originates may emerge.

Finally, to ensure that details on changes in relative intensities of various transitions to the same sub-shell are properly dealt with, calculations and integration over the multilayer structure are carried out for each transition individually.

Taking into account the complexity of the problem, in this work the presentation is limited to the description of the model in the case where secondary X-rays are induced by primary characteristic X-rays, and to its implementation in the DT2 package.^16–19 In related studies, to be published in the near future (parts II and III), applications and the problem of secondary X-rays induced by electrons provenant from the non-radiative transitions following the initial collision of beam particles, will be addressed.

2 PIXE target X-ray yield

2.1 Thin targets

When considering thin targets under particle irradiation, the number of X-rays, N_j,_Zi, detected from rearrangement transitions j (K_α, L_α,…) of element Z_i can be written as:beingwhere

is the detector solid angle fraction, ε_det,j and T_sis,j are the energy dependent detector efficiency and the transmission coefficient of the absorbers placed between the sample and the detector, respectively, for the X-rays emitted by transitions j of element Z_i. N_p is the number of particles used in the irradiation, C_pp is the charge per particle in μC, b_cs is the particle beam cross-section and ψ_inc is the incidence angle defined between the beam direction and the normal to the target surface.

is the target total X-ray yield, for transition j of element Z_i, per μC for a target irradiated by E_p energy particles, and includes the mass fraction of element Z_i in the target, f_Zi. Finally, N_Av is Avogadro's number, M_at,Zi the molar mass of element Z_i, the X-ray production cross-section in barns for particles of energy E_p and ξ is the sample areal mass in μg cm⁻², frequently referred to as thickness even though it does not have the dimensions of a distance. The value of has been calculated from the revised SI standard based on the 2017 CODATA revision.²⁰

It is important to emphasize here that the mass fraction, f_Zi, of element Z_i is not included in ξ, but is kept separate on purpose both to allow it to be treated as an unknown in analytical processes, or to serve as a parameter in system calibration operations.

2.2 The equivalent thickness concept

When dealing with thick or half-thick targets, the calculation of the total X-ray target yield is not so straightforward. In these cases, as the ion beam particles penetrate the target, they lose energy, which changes their X-ray production cross-section, since E_p is reduced, and the induced X-rays are absorbed before exiting the sample. The target total X-ray yield for any transition j originating from any element Z_i must now be determined by integrating the differential effective yield density. Still, introducing the concept of equivalent thickness,¹³ ξ_eq,j,Zi(E_p),

eqn (1) still allows the calculation of the target total X-ray yield, for thick and half-thick targets as:

In eqn (3) x(E) is the penetration depth variable defined as the distance of a given point along the particle penetration path and the sample surface, measured along the ion beam path.

T_j,Zi(x(E)) ≡ T_j,Zi(x) is the absorption of X-rays j of element Z_i while travelling from penetration depth x ≡ x(E) to the surface of the sample, and is the ion beam particle energy loss derivative.

Normalizing to the incident energy X-ray production cross-section allows the total thick target yield to be formally written in the same way as for thin targets by replacing the target thickness by the equivalent thickness. The main difference is that, while the thin target surface area is independent of the X-ray being measured, the equivalent thickness is different for each X-ray.

2.3 Non-homogeneous targets

The use of the equivalent thickness concept allows the expression for the target total X-ray yield to be extended, even for the general multilayer case. Still, it is important to realize that now the elemental mass fraction in each layer must be also included in the definition because it changes from layer to layer.

Making the transmission of element Z_i j transition X-rays from layer n to the surface, the result is:

Last but not necessarily least, even if the sample is not laterally homogeneous and/or if the detector size or detector sample distance leads to transmission terms or layer structure description that depends on the y, z positioning of the beam on the sample, the concept although becoming a bit abstract, can still be used to establish the following general expression for the PIXE target yield of general targets irradiated by particles of E_p energy, if a set of homogeneous (y_a, z_b) regions can be established to describe the sample:

beingand

3 Secondary fluorescence penetration function model

3.1 Primary point emission and cylindrical symmetry

Consider Fig. 1 representing a homogenous target. As the ion beam particles penetrate the target, at any given penetration depth x₁, X-rays are induced and emitted in all directions. A fraction of these, say B^prim_β, travels in the direction of the detector, while others, say A_α, are emitted in another direction and may be absorbed in the target material, say in a volume dV at a distance r from the x₁ position, and also induce the emission of B_β X-rays, which may be emitted in the direction of the detector and contribute, with say B^sec_β X-rays, to the B_β peak in the measured spectrum. In this case, the A_α and B_β X-rays produced at x₁ position are named primary X-rays. The B_β X-rays produced in volume dV at position are named secondary fluorescence X-rays, and some of these may add to the primary B_β X-rays reaching the X-ray detector, enhancing the target total X-ray yield for B_β X-rays.

Fig. 1 Primary X-rays A_α produced at a penetration depth x₁ lead to the emission of secondary X-rays B_β in volume dV that add to primary B_β X-rays, enhancing their target yield. The X-ray emission process is assumed to have cylindrical symmetry and therefore be possible to describe using a simple 2D image (b). This is so, even if the ion beam direction is outside the detection plane defined by the sample normal (in yellow in the images) and the direction defined by x₁ and the detector. The angle α between the incidence plane defined by the beam and the sample normal (a) and the detection plane can take any value (check main text for details). Still, for the presented calculations to be valid, it must be possible to assume the samples as infinite and homogeneous in the plane perpendicular to the sample normal (e.g., planes e and f in the left image).

Fig. 1(a) represents the ion beam incident in a direction that may be not contained in the detection plane defined by the normal to the sample surface (shown in yellow in both Fig. 1(a) and (b)) and the line connecting point x₁ and the detector. Assuming that any relevant distance r is small relative to the distance between x₁ and the detector, so that ψ_det can be assumed as constant and independent of r, the circular symmetry around the sample normal can be assumed for all the detection processes, even if the irradiation beam is not in the detection plane. This is so because point x₁ is the single common point for both irradiation and detection processes. Furthermore, if the target can be considered laterally homogeneous (meaning that layers are infinite and homogeneous in the planes parallel to the sample surface), all points x_n (along the beam path) outside of the detection plane may be assumed, for all calculation purposes, to be equivalent to their projection on the sample normal.

In the case of complex wide angle detector geometries the whole approach still applies, although numerical integration over the various ψ_det values will now be required.

The need for numerical integration in these cases is not a restriction of the secondary fluorescence correction process, but is also required to properly determine matrix corrections processes affecting the primary X-ray yield, as mentioned in the previous section.

3.2 Secondary X-ray fluorescence cross-section

In order to determine the total amount of secondary B^sec_β X-rays, it is necessary to start by writing the differential density function, dX_BβAα(x₁), describing the conversion of X-rays A_α produced at a penetration depth x₁ into secondary X-rays B_β (the “sec” superscript will be omitted for simplicity of writing) that reach the target surface after being induced in the volume element dV. The following expression may be used as a starting point:where:

• is the primary A_α X-ray production density function at penetration depth x₁;

• T_Bβ(x₁, r, θ) is the transmission factor of B_β X-rays from the volume dV up to the target surface, calculated for the detector direction;

• is the conversion probability that A_α X-rays absorbed in element B in sub-shell η are converted into B_β secondary X-rays; and

• is the cross-section for an A_α X-ray to be absorbed at a distance r away from the emission point x₁.

The primary X-ray production term corresponds to the differential terms in the expressions presented in the previous sections, which were also partially addressed in the previous subsection.

It therefore remains to discuss the other terms, whose product may be referred to as the secondary fluorescence cross-section for the conversion of A_α primary X-rays in B_β secondary X-rays.

Still, before any other discussion, it is necessary to address the lack of an explicit term for the element B mass fraction, f_B, in eqn (10), which is needed to add the term resulting from this exercise to eqn (3), to obtain an appropriate expression for an equivalent thickness secondary fluorescence correction, since ξ_eq,j,Zi(E_p) has no mass term.

3.2.1 and specific conversion probabilities. Obtaining this explicit mass term can be done by factoring out the conversion probability component. Only the cases where the absorbing and emitting shell of B are the same will be considered, because the number of secondary B_β X-rays emitted from transitions to a shell different from the shell absorbing the primary A_α X-rays is, in most cases, not relevant compared to the primary B_β X-rays produced in that shell. The cases where this is not valid are just the situations where the ion beam particles either do not reach the fluorescence layer containing the emitter of the B_β X-rays, or reach it having already lost most of its energy, while the B_β X-rays can still significantly emerge from the sample towards the detector. This is a condition which is probably extremely rare in practice.

This being set, the factoring out of the term in the simplest case, an absorbing K-shell, can be obtained based on the following result:

where σ^photo_K,B(A_α) is the K-shell photoelectric ionization cross-section of B for A_α X-rays, ω_K,B is the K-shell fluorescence coefficient of B and k_β,B is the branch ratio of transition β out of all radiative transitions to the K-shell of the B element.

The mass fraction term, f_B, can now be factored out in order to establish a specific probability, independent of both the mass fraction and the mass absorption coefficient of the A_α X-rays, namely:

Eqn (10) can now be written with f_Bβ factored out, namely as:

In the case of L and M sub-shells the situation is a bit more complex; nevertheless, the same reasoning as used for the K-shell applies since changes are only present in the photoelectric cross-section term. The generalization of eqn (12), defining can therefore be made.

K, L and M sub-shell photoelectric ionization cross sections are normally approximated^8,11 based on the total absorption cross-section of the X-ray energy, σ_T,B(A_α) ≡ μ_Aα, and the jump ratios, S_η, for the sub-shell. Taking E_i to be the X-ray or the sub-shell ionization energy, as applicable, the following results apply to K, L and M sub-shell fluorescence (omitting the A_α term for simplicity):

If E_L1 < E_Aα < E_K

If E_L2 < E_Aα < E_L1

If E_L3 < E_Aα < E_L2

If E_M1 < E_Aα < E_L3

If E_M2 < E_Aα < E_M1

If E_M3 < E_Aα < E_M2

If E_M4 < E_Aα < E_M3

If E_M5 < E_Aα < E_M4

3.2.2 Secondary fluorescence production and survival. Using the definitions of the primary A_α X-ray production density function and the conversion probability the differential density function describing the conversion of primary A_α X-rays produced at a penetration depth x₁ into secondary fluorescence B_β X-rays at the volume element dV, which reach the target surface, dX_BβAα(x₁), as defined in eqn (10), may be rewritten as:

In order to obtain the density function for secondary B_β X-rays emerging from the target surface towards the detector, due to secondary emission induced by primary A_α X-rays emitted at penetration depth x₁ it is necessary to integrate eqn (38) over the whole target volume, and we can use this step to define the corresponding specific density function, χ_BβAα(x₁) by dividing by f_B; the result obtained is:

The integral in eqn (39) represents the fraction of primary A_α X-rays that may be converted to secondary B_β X-rays, and, if that happens, will survive until reaching the target surface.

3.2.3 The function and SFC equivalent thickness. In order to properly analyse eqn (39), it is important to focus on the differential under the integral:

This is the differential cross-section for a primary A_α X-ray produced at the penetration depth x₁ to be absorbed at a distance r from x₁ and converted into a secondary B_β X-ray that reaches the target surface along a trajectory that leads to the X-ray detector.

Although it looks simple, there are a few details, including theoretical ones, which are worth taking into account carefully.

The most critical term, even if it may not seem so, is the detailed description of the absorption of A_α X-rays in the differential volume. Using spherical coordinates, there are two main components in this process. A geometrical one that is related to the angular description, which leads to a term in the angular variables, namely r²sin(θ)dθdϕ, and a second term related to the ionization process itself.

Since X-rays vanish when interacting with atoms to produce ionization, as opposed to what is observed with ions, which just lose energy but do not vanish, the number of matrix atoms ionized is proportional to the number of absorbed X-rays.

Considering a small slab of thickness Δr → dr this results in the following expression for the number of absorbed A_α X-rays, using a first order Taylor series approximation:

N_X(Aα)(0) is the number of X-rays reaching the slab. The absorption in volume dV therefore contributes with an overall term given by μ_Aαr²sin(θ)drdθdϕ.

Taking into account that this expression makes use of the number of X-rays reaching the slab, a term describing the loss of intensity of A_α X-rays between the emission point x₁ and the absorbing volume dV, must be considered. Therefore, the differential cross-section for an A_α X-ray to be absorbed in volume dV at a distance r away from the emission point x₁ is:

The remaining term to be mentioned is the probability that the B_β X-rays emitted in the elemental volume dV in the direction of the detector reach the target surface. With μ_Bβ being the target mass absorption coefficient and x₁ and r expressed in consistent units, usually areal mass units, the result is:

Finally, writing the whole term in spherical coordinates, for a homogeneous target (see Fig. 1), the result is:

Therefore the final expression is:

The differential may be referred to as the secondary fluorescence differential cross-section for the conversion of A_α X-rays into B_β X-rays that emerge from the target in the direction of the detector.

The function defined as the integral of over the whole target volume is the secondary fluorescence target yield, emitted in the direction of the detector, originating from the conversion of A_α X-rays into B_β X-rays, and corresponds to the integral in eqn (39).

Using the fact that has cylindrical symmetry, the integral can be immediately integrated in ϕ by taking the x axis as being along the normal to the target surface. Note that the x axis for calculating the integral in eqn (39) is independent of the definition of x₁ along the ion beam penetration path and therefore the x axis for this calculation can be set freely. The result after integrating over ϕ is:

Summing over all primary A_α X-rays produced at the penetration depth x₁ and leading to B_β secondary X-rays, the specific secondary fluorescence correction density function, χ_Bβ(x₁), can be written as:

and added to the equivalent thickness definition, leading to a secondary fluorescence corrected equivalent thickness, which now depends not just on the X-ray being detected, but also on the various other X-ray emitters present in the target:

4 The function analytical solution

4.1 First steps for solving the integral analytically

Considering that r is small relative to the distance to the detector, so that it is possible to assume that the detection angle ψ_det is constant relative to r, the integral in eqn (46) can be solved analytically, as long as it can be assumed that the sample is infinite and homogeneous in all planes normal to the surface normal, at least for the fluorescence process. This means that the model may be easily adapted for application to a small inclusion emitting primary X-rays if the particle beam is kept within it, but it is not applicable to the case of a small inclusion emitting secondary X-rays due to primary X-rays originating in its surroundings.

In order to obtain the analytical solution, it is important to start with a change of variable, namely by setting:

The expression can then be simplified to:

being

Calculating the integral in eqn (50) is better done by separating the full integral in six different cases according to the relations between d and t described in Table 1 (see Fig. 1 for variable references).

Table 1 Integration limits for r and θ for the various integrals (I₁ to I₆)

	i	rmin,i	rmax,i	cos θi	cos θf
d ≤ t/2	1	0	d	1	−1
	2	d	t – d	1	−d/r
	3	t – d	∞	(t – d)/r	−d/r
d > t/2	4	0	t – d	1	−1
	5	t – d	d	(t – d)/r	−1
	6	d	∞	(t – d)/r	−d/r

---

Equation eqn (50) is thus better written as follows (i values according to Table 1):

being:

Setting now:

the first three integrals mentioned above become

Expanding these expressions leads to

Now, Gradshteyn²¹ states

Both Gradshteyn²¹ and Abramowicz²² define the exponential integral as:

Gradshteyn further sets for negative values of x: while Abramowicz²² defines the exponential integral of order 1 for positive values of the variable as:

leading to the relation:

These Abramowicz definitions having x ∈ ]0, ∞[ will be used for the remainder of this work.

The exponential integral and the exponential integral of order 1 may also be presented as power series as follows:

where γ = 0.57721156649… is Euler's constant.

Now in eqn (58) the signs of the constants in the exponential and sinh() function are well defined since both the mass absorption coefficients and the distances are positive.

Before applying eqn (61) to (58) it is still important to obtain a few additional expressions.

Assuming a > 0, b > 0 and x > 0, from eqn (61) using Abramowicz nomenclature, it is important to note that:

If a = b Gradshteyn in its equation 2.484.6 (ref. 21) further states:

that converting to Abramowicz nomenclature becomes:

Before applying these expressions to eqn (58) and other integrals, it is important to check the case where x → 0, since in this condition, ln(x), E_i() and E₁() are divergent. The limit differences are:

Setting a = μ_Aα and and applying these to the definitive integral I₁ the result is:

It is important to note that, from the above equations, it results for all these cases:

These results can be written in a more condensed and physically interesting form, namely:

I₁ = E₁(2μ_Aα·d) + ln(2μ_Aαd) + γ, if g_ = 0 (78)

Addressing the calculation of the definitive integral I₂, eqn (72)–(74) are not applicable to eqn (59) and the definitions in eqn (63)–(65) must be used directly. Setting ζ > 0 and η > 0 and y = at the result is:

Applying this to eqn (59) the result is:

4.2 Infinite thickness targets

In the case of thick targets, t = ∞ and only I₁ and I₂ apply. Adding eqn (76)–(78) and (81) provides

The computational implementation of these results must take into account that for very small values of the argument, the exponential integral diverges due to the term in ln(x) in eqn (65). Still, in the case of small values of d (x₁ still close to target surface) no problems arise since the results are:

Before proceeding to deal with half-thick targets, it is still important to check the theoretical possibility that d is not too small but either |g₋| is too small but not enough to make the product |g₋|·d too small, or the mass absorption coefficient of the A_α X-rays is so small that the product μ_Aα·d → 0. In all these cases numerical calculation problems emerge linked to eqn (82)–(84). Besides, the problematic conditions in |g₋| may also combine with those on μ_Aα and therefore all cases must be addressed carefully.

Taking into account the power series expansions in eqn (65) the results for |g₋| → 0 while the product |g₋|·d does not, are:

As could be expected this expression is identical to that of eqn (84) since in the limit g₊ = 2μ_Aα. In the case where μ_Aα·d → 0, two conditions can be found, namely, |g₋| → 0, or not so and g₋ < 0. In the first case, the limit of eqn (84) is ln(2). In the second case, it is necessary to establish an ad hoc cut-off, say corresponding to a 95% intensity decrease of A_α X-rays, which causes eqn (58) and (59) to become:

which, when summed, provides:

If in this case d → 0 this equation is also not valid. Using the power series expansions leads to:

4.3 Homogeneous half-thick targets

4.3.1 Primary X-rays emitted before half-layer depth. In the general case of half-thick targets, all six integrals must be calculated. As can be seen from Table 1, the six integrals are separated into two distinct cases. Integrals I₁ to I₃ provide the results for the situation where the point x₁ exists at a distance to the target surface less than half of the target thickness, and integrals I₄ to I₆ provide results for the situation where this is not so and therefore d > t/2.

In the case of I₃, eqn (63) should be applied directly to eqn (60), the result being:

therefore, setting the result for I₂ + I₃ is:

I₂ + I₃ = e^b·dE₁(μ_Aα·d) − e^{−b·(t−d)}E₁[μ_Aα·(t − d)] − [E₁(g₊·d) − E₁[g₊·(t − d)]] (94)

For the homogeneous half-thick target and d < = t/2 the sum of I₁, I₂ and I₃ results in:

and since g₋ = 0 ⇒ g₊ = 2μ_Aα,

I₁ + I₂ + I₃ = E₁[2μ_Aα·(t − d)] + ln(2μ_Aαd) + γ + e^b·dE₁(μ_Aα·d) − e^{−b·(t−d)}E₁[μ_Aα·(t − d)] if g₋ = 0 (97)

In this case, when d → 0 the result for all three possibilities is the same, namely:

In what concerns other extreme cases, as in the previous subsection, we may find μ_Aα → 0 while d is not too small. Once again we can have two different conditions for this. In the case when |g₋| → 0, the limit of the sums is 0 because if μ_Aα → 0 and |g₋| → 0 then b → 0. If it is instead g₋ < 0, then b is no longer a vanishing value and a cut-off must be used to calculate the I₃ integral and (t – d) must replace the cut-off in eqn (91), leading to the result:

If now both d → 0 and μ_Aα → 0 eqn (99) must be used to calculate the limit and the result is:

4.3.2 Primary X-rays emitted after half-layer depth. Since the sum of I₁ to I₃ is only valid for d ≤ t/2, when the contrary is true, meaning when d > t/2, the sum of integrals I₄ to I₆ applies.

Based on Table 1 these are:

In the case of I₄ given the formal identity to I₁ once d is replaced by (t – d), the result is:

= E₁[2μ_Aα·(t − d)] + ln[2μ_Aα(t − d)] + γ, if g₋ = 0 (104)

In the case of I₅, expanding the expression in eqn (101) provides:

Taking into account eqn (79) and (80) three results are possible for I₅, namely:

In the case of I₆ the result is:

Adding I₄, I₅ and I₆ the results are now:

In this case, when (t – d) → 0 the results for the three possibilities are:

As in the previous case, it is also important to address the potential extreme conditions where μ_Aα → 0 while (t – d) is not too small. As before the two situations that may be addressed are |g₋| → 0 and g₋ < 0. In the first case, the result of the sum of integrals is 0 as it was also for the condition d ≤ t/2. In the case of g₋ < 0 the result is:

If now also (t – d) → 0 applies, the limit of this eqn (116) must be used, the result being:

4.4 The general layered target case

If the target is more complex than a single homogeneous layer and made up of several physically distinct layers as drawn schematically in Fig. 2, calculating secondary fluorescence processes for PIXE experiments becomes a bit more complex and, as far as the author knows, this work is the first time a systematic, general global solution is presented in standard literature.

Fig. 2 Primary X-rays A_α produced at a penetration depth x₁ in layer t_e induce the emission of secondary X-rays B_β in volume dV in layer t_f. These will appear as if produced at depth x₁, and will either enhance B_β X-ray target yield from layer t_e or create a “phantom” presence of element B in layer t_e. (a) Case 8: the secondary fluorescence is produced in a physical layer located deeper into the target than the layer emitting the primary X-rays. (b) Case 9: the secondary fluorescence is produced in a physical layer located less deep in the target than the layer emitting the primary X-rays. In both cases it is once again assumed that the sample layers are homogeneous and infinite in the plane perpendicular to the sample normal (shown in yellow).

In this case, three different situations can be faced in respect to secondary fluorescence: (a) the secondary X-rays are produced in the same layer as the primary X-rays, (b) the layer emitting secondary X-rays is located deeper into the target than the primary X-rays layer or (c) the layer emitting secondary X-rays is closer to the target surface than the primary X-rays layer.

In case a, or 7 since it follows integral I₆, eqn (47) and (50) need just a slight change to cope with the extra layers that may be present between the emitting layer and the target surface, the result being:

Case a (or 7): making

based on eqn (50)

In cases b, or 8, and c, or 9, the situation is different because it is necessary to account for three facts, namely, (i) the primary A_α X-ray absorption between the emission point x₁ and the absorption volume V_fl is not homogeneous, (ii) the path of B_β X-rays from the integration volume up to the surface of the layer where secondary fluorescence effects are taking place has a different expression from the one defined in eqn (43) used for the case of the single homogeneous layer and case (a) of multilayered targets, and (iii) a single more complex integral expression applies.

In cases (b) and (c) eqn (43) and (45) need to be re-written. In order to simplify the expressions both for easy reading and for a clear understanding, some definitions are presented in Table 2.

Table 2 Integration limits for r and θ for the integrals I₈ and I₉. The integral limits are the same due to the fact that angle θ was defined as the smallest angle to the normal, in both cases

	s	t^sre	t^sref	t^srf
a The value of is an ad hoc cut-off taken as the value above which less than 5% of primary X-rays exit the emission layer.
d < t^begf	8(b)	t^ende – d	t^begf – t^ende	r cos(θ) – (t^8re + t^. 8ref)
d > t^endf	9(c)	d – t^bege	t^bege – t^endf	r cos(θ) – (t^9re + t^. 9ref)

---

xsi=cos(θ)	rs,min	rs,max	ζi = cos θi	ζf = cos θf
d < t^begf			1	^a
d > t^endf			1	^a

---

Based on these definitions and on Fig. 2, eqn (45) can be promptly adjusted (note that t^s_re is the fraction of the emitting layer crossed by A_α X-rays, and t^s_rf is the fraction of the layer absorbing the A_α X-rays, crossed by these) leading to the following results:

The following expression replaces eqn (47):

being in this case,

Integrating the above differential expressions having the integral limits defined in Table 2 provides:

Further simplification will result from applying the following change of variables to the integrals in eqn (125) and (126):

ζ = cos(θ); dζ = −sin(θ)dθ; (129)

which using the changes of variable mentioned above results in:




therefore,which leads to the following integrals that are solved numerically using Gaussian methods.All terms in the exponential having been grouped together to avoid numerical integration problems.

The final expression for the number of B_β X-rays emitted by a layered target, whose structure may be simulated as a set of layers parallel to the surface, and infinite in the directions perpendicular to the sample normal, can now be written as:

beingwhere, with provided by eqn (132), is:

It is important to note that solving eqn (132) numerically adds an additional set of sums to the ones already introduced by eqn (136), combined with eqn (135), which must be carefully implemented.

Note that now, because the mass fraction term must be included in the definition of the equivalent thickness, it cannot be just put in evidence, as was done in eqn (39).

This is not a problem for simulations, but is a complex situation to address if the problem in question is the exact fitting of spectra of unknown samples. In the present work, this issue is not addressed beyond this statement, still it is a subject that will be addressed in the applications part of this trilogy.

5 The general case expression

Summing up all previous results, it is possible to write a global expression for the most general case possible, namely for the PIXE yield of a wide spot or wide detector that requires a generalized sum over a set of (y_a, z_b) pairs.

It is nevertheless important to ensure that homogeneous conditions are verified within each partial spot (y_a, z_b), as otherwise the expression cannot be used without detailed adaptations that have not been presented in this paper, even if they may eventually be derived from the results presented here.

Starting from eqn (7)–(9) and adding up the secondary fluorescence terms, the final result is:

beingand

In these equations, refers to the homogeneous cases and case a (or 7) of eqn (118) and refers to inter-layer secondary fluorescence, cases b and c (or 8 and 9), as described by eqn (123).

6 Implementation and analysis

6.1 Homogeneous targets

Once obtained these results, their computational implementation is reasonably straightforward, the single remaining issue needing some attention being the cases where variables take very small values so that limit expressions must be used.

The implementation was made as additional code to the previous DT2 code,^16,19 which was designed from the start to allow the handling of multilayered targets.²³

6.1.1 The infinite target case. In the infinite target case, eqn (82) and (83) are used as long as the following expression is true:

d > 10⁻⁵ ∧ g₋ > 10⁻⁵ ∧ μ_Aα > 10⁻⁵ (140)

If this expression is not true, then each condition must be taken into account individually. Table 3 lists the conditions, equations and limit cases replacement when dealing with infinite (thick) targets.

Table 3 Equation selection table for the case of infinite targets

d	|g−|	|g−|·d	μAα·d	g−	Equation
>10^−5	>10^−5	—	>10^−5	g− > 0	Eqn (82)
				g− < 0	Eqn (83)
			≤10^−5	g− < 0	Eqn (91)
	≤10^−5	>10^−5	>10^−5	—	Eqn (88) ≡ (84)
			≤10^−5	—	ln(2)
		≤10^−5	>10^−5		Eqn (84)
			≤10^−5	—	ln(2)
≤10^−5	>10^−5	>10^−5	>10^−5	—
			≤10^−5	g− < 0	Eqn (92)
		≤10^−5	—	—	ln(2)
	≤10^−5	—	—	—	ln(2)

---

Simulations corresponding to one of the alloy cases presented in the 1992 paper¹² are shown in Fig. 3. In this case the BCS S387 iron–nickel standard was considered. The spectra shown correspond to simulations assuming 1.65 MeV proton beam irradiation, replicating the experimental conditions used in the 1992 study. Simulations were also carried out for proton beams of 1.1 MeV and 2.5 MeV. In Fig. 4 the changes in percentage correction determined as a function of beam energy are presented for the five elements exhibiting the most significant effects. It can be seen that as ion beam energy increases, the necessary correction also increases. The results are different from those presented in the 1992 (ref. 12) paper because the present work uses a penetration function method and Gaussian integration, which accounts for the whole sample, as used in the 1996 paper¹³ and not the Simpson integration over pairs of irradiated numerical layers (similar to the Ahlberg et al. method⁵) used in 1992. The present results for this homogeneous thick target are, therefore, identical to those found in the 1996 paper. By applying the correction factors presented in Fig. 4 for 1.65 MeV, to the experimental data published in Table 3 of ref. 12, relative differences of 1.7%, 0.78%, 5.0% and 1.33% are found now between secondary fluorescence corrected data and reference values for Ti, Cr, Mn and Fe respectively. Taking into account that the reference values have uncertainties of 4%, 0.64%, 5.0% and 0.55% respectively, it can be concluded that the results obtained after secondary fluorescence correction fully agree with the standard reference data.

Fig. 3 Overlap of simulated spectra for 1.65 MeV proton irradiation of the BCS_SS387 reference material taking into account secondary fluorescence corrections (w/SFC) and not considering these (no SFC). It can be seen that differences are observable in the most intense peaks, but not so much in the others. In this case, the most intense SFC effect is observed for Cr at 5.4 keV, with an effect of 11.6%, while Fe at 6.4 keV presents a SFC effect of 7.4%.

Fig. 4 Change in the percentage of secondary fluorescence correction (%SFC) counts on the total counts in the area of the X-ray peaks simulated for five different chemical elements, as function of proton beam energy. It can be seen that for all these cases, the %SFC increases as ion beam energy increases.

Secondary fluorescence correction situations may, nevertheless, be significantly different from this. Testing as examples some potentially complex cases such as MoP, PbCrO₄, Ti82.5–Mo10–Mn2.5 and Co10–Cu90, under 1.65 MeV proton irradiation, different cases can be observed.

In the case of low energy X-rays, namely P-K, Mo-L and Pb-M, no meaningful secondary fluorescence corrections are observed; the most intense case is Mo-L_β1 that shows a 1.86% increase under irradiation of a bulk Ti82.5–Mo10–Mn2.5 sample. The difference in energy between Pb L lines and the Cr–K absorption edge results in a photo-electric absorption cross section that is too low for a significant effect to be observable in PbCrO₄.

In the Co10–Cu90 case, a different situation applies and secondary fluorescence corrections for Co K_α lines from 18% to 30% are found. The effect visible in the Co K_α peak height, for a proton irradiation at 1.65 MeV, is shown in Fig. 5.

Fig. 5 Simulation of a 10 wt% cobalt alloy in copper: overlap of SFC corrected and not corrected spectra (left) and change in the percentage of secondary fluorescence correction (%SFC) counts on the total counts in the area of the Co X-ray peaks simulated as function of proton and He beam energy. It can be seen that the %SFC increases as a function of ion beam energy are stronger for proton beams than for the He beams.

6.1.2 The half-thick target case. In the case of the half-thick targets, as shown in Section 3.3, calculations are a bit more complex, and so is their implementation. The selection of equations for this case is summarized in Tables 4 and 5, since two different sets of integral sums must be dealt with.

Table 4 Equation selection table for the sum I₁ + I₂ + I₃ (primary emission point located before the layer half-thickness) for the case of homogeneous half-thick layer targets

d ≤ t/2	|g−|	|g−|·d	μAα·d	g−	Equation
>10^−5	>10^−5	—	>10^−5	g− > 0	Eqn (95)
				g− < 0	Eqn (96)
			≤10^−5	g− < 0	Eqn (99)
	≤10^−5	—	>10^−5	—	Eqn (97)
			≤10^−5	—	0
≤10^−5	>10^−5	>10^−5	>10^−5	—	Eqn (98)
			≤10^−5	g− < 0	Eqn (100)
		≤10^−5	—	—	0
	≤10^−5	—	—	—	0

---

Table 5 Equation selection table for the sum I₄ + I₅ + I₆ (primary emission point located beyond the layer half-thickness) for the case of homogeneous half-thick layer targets

t – d ≤ t/2	|g−|	|g−|(t – d)	μAα(t – d)	g−	Equation
>10^−5	>10^−5	—	>10^−5	g− > 0	Eqn (110)
				g− < 0	Eqn (111)
			≤10^−5	g− < 0	Eqn (116)
	≤10^−5	—	>10^−5	—	Eqn (112)
			≤10^−5	—	0
≤10^−5	>10^−5	—	>10^−5	g− > 0	Eqn (113)
				g− < 0	Eqn (114)
			≤10^−5	g− < 0	Eqn (117)
	≤10^−5	—	—	—	Eqn (115)

---

Applying these to the simulation of the most intense case shown in the previous section, namely the cobalt copper alloy, it can be seen that the secondary fluorescence correction in thin targets is not zero, but it decreases significantly with thickness as well as with ion beam energy.

In Fig. 6 it can be seen that the secondary fluorescence correction increases as a function of beam energy (as already observed for thick targets) as well as the target thickness.

Fig. 6 Simulation of 10 wt% cobalt alloy in copper: change in the percentage of secondary fluorescence correction (%SFC) counts of the total counts in the area of the Co X-ray peaks as a function of beam energy and target thickness.

Although not shown in the graph, He ions at 2500 keV are fully stopped in 3.2 and 6.4 mg cm⁻² targets, and the same applies to He 5000 keV and proton 1100 keV beams in the 6.4 mg cm⁻² target. Still, out of these four cases, only for the He 2500 keV beam in the 6.4 mg cm⁻² target is the secondary fluorescence correction identical to that of the thick target.

This results from the fact that secondary fluorescence effects that take place beyond the ion beam range, still affect the overall spectra.

6.2 Layered targets

If the target is not thick enough but composed of more than a single homogeneous layer, secondary fluorescence may be induced in the same region or in regions different from the one where the primary X-rays are emitted.

As presented in the previous section, the complexity of the case requires that in the second case, the integrals involved must be solved numerically.

The first of these cases, which involves calculating secondary fluorescence effects taking place in the same physical layer as the primary X-rays emission, is handled using eqn (118) and apart from the absorption term and the shift in the penetration value relative to the layer surface, nothing is changed relative to the homogeneous half-thick layer target case.

The second of these two conditions involves the emission of secondary fluorescence X-rays from layers different from that emitting the primary X-rays.

In this case, two conditions can arise, namely either the layer emitting secondary X-rays is deeper than that emitting the primary X-rays, or vice versa.

In each of these situations, eqn (132) applies and the only numerical extreme issue that must be overcome is the occurrence of vanishing cosine values, which is resolved by setting an ad hoc cut-off as mentioned in Table 2.

Although this extreme value problem is minor in this case, it is still necessary to take into account and overcome a large number of embedded sums, which must be managed to ensure proper implementation of the general case calculation.

In order to illustrate these types of conditions, simulations were run for a combination of layers and substrate materials, specifically MoP and Co10–Cu90 alloy. As shown, in the case of MoP bulk, secondary fluorescence induced in P by Mo-L lines is small relative to the direct primary induction of P X-rays. If a film of Co10–Cu90 alloy is set on top of it, not much difference is observed even though the secondary fluorescence in P increases to roughly 11%. In Fig. 7 the effect of a 1.6 mg cm⁻² film of Co10–Cu90 placed on top of a bulk MoP substrate is shown.

Fig. 7 Simulation of a 1.6 mg cm⁻² film of 10 wt% cobalt alloy in copper placed on top of a bulk MoP substrate.

Still, if the order of the materials is exchanged, a different image can be found. In Fig. 8 the change of the effect observable as a function of top layer thickness is presented for both the MoP layer on top and the other way around; the figure also includes a comparison of the simulated spectra for a multilayer sequence of 0.8 mg cm⁻² MoP and Co10–Cu90 films starting with MoP, using three times less charge.

Fig. 8 Spectra of a 0.8 mg cm⁻² MoP film on top of a 0.8 mg cm⁻² Co10–Cu90 film, a Co10–Cu90 film on top of a MoP film and a multilayer sequence of 3 pairs of MoP/Co10–Cu90 films. The differences can be observed to be very significant, as expected.

It can be seen that important differences are observed. A systematic validation of these results is necessary to ensure that both theoretical work and software implementation are working properly, before the results presented here can be used systematically. Still, the report of this validation will be presented in part II.

7 Conclusions

Simulation of PIXE spectra is a useful tool for various purposes, ranging from the simplest task of teaching PIXE without access to an accelerator to its unavoidable use for analysis of data from Total-IBA¹⁵ experiments.

PIXE spectral reproduction is available through a few computer codes described in the literature, such as GUPIX,²⁴ GeoPIXE⁹ or LibCPIXE²⁵ but to the best of author's knowledge, up until the present paper, no available computer code was able to deal with simulation and secondary fluorescence corrections of multilayer samples where the same chemical element may be present in more than one layer.

As far as the author is aware, a general and global theory presented here to deal with X-ray induced secondary X-ray fluorescence in PIXE experiments under such general conditions was not previously available in standard and easily accessible literature before this work.

The present algorithms are implemented in the new version of the DT2 code (DT2F_0v9_98), therefore corresponding to a major upgrade of its prior versions.^16,19

Data availability

The data used in the present work are entirely generated by the revised version of the DT2 computer code mentioned in the Conclusions section. At the present moment, the generated data used in this work are available just upon request to the author. Still, in the near future a repository will be created for it and the possibility of making the simulation software available as freeware is under evaluation. In the event of a positive response to this evaluation, the freeware version of the executable and the usage conditions will be made available in the same repository.

Conflicts of interest

The author has no conflict of interests related to this work.

Acknowledgements

The present work was partially funded by the FCT UID/Multi/04349/2020 project. The global PIXE secondary X-ray fluorescence corrections theory and its computational implementation, which might have lacked in PIXE standard and easily accessible literature for a long time now, was requested several times by Chris Jeynes as a needed keystone for the PIXE component in Total-IBA.¹⁵ A special acknowledgement is therefore due to him since his pressure for it was one of the important drives for carrying out such a long, detailed and necessary work. Even though in many practical situations faced the effect in complex multilayered targets is small, as shown, situations do emerge where that is not so, and in any case, uncertainty budgets calculations do need the results presented in this work even for conditions where the correction is small.

Notes and references

1. T. B. Johansson, A. Akselsson and S. A. E. Johansson, X-ray Analysis: Elemental Trace Analysis at the 10⁻¹²g Level, Nucl. Instrum. Methods, 1970, 84, 141–143.
2. N. P. Barradas, M. Mayer, M. A. Reis and F. Schiettekatte, IBA Software, in Ion Beam Analysis, Fundamentals and Applications, ed. M. Nastasi, J. W. Mayer and Y. Wang, CRC Press, Taylor and Francis Group, LCC, Boca Raton, FL, USA, 2015.
3. F. Folkmann, C. Gaarde, T. Huus and K. Kemp, Nucl. Instrum. Methods Phys. Res., 1974, 116, 487.
4. W. Reuter, A. Lurio, F. Cardone and J. F. Ziegler, J. Appl. Phys., 1975, 46, 3194.
5. M. S. Ahlberg, Enhancement in PIXE analysis, Nucl. Instrum. Methods, 1977, 142, 61–65.
6. F. N. Richter and U. Wätjen, PIXE calibration and correction of matrix effects in the case of thick samples, Nucl. Instrum. Methods Phys. Res., 1981, 181, 189–194.
7. B. van Oystaeyen and G. Demortier, Matrix effects in PIXE evaluation of the major components in thick homogeneous samples, Nucl. Instrum. Methods, 1983, 215, 299–313.
8. J. L. Campbell, J.-X. Wang, J. A. Maxwell and W. J. Teesdale, An exact treatment of secondary fluorescence and terciary fluorescence enhancement in PIXE, Nucl. Instrum. Methods Phys. Res., Sect. B, 1989, 43, 539–555.
9. S. H. Sie, W. L. Griffin, G. F. Suter, C. G. Ryan and D. R. Cousens, Quantitative PIXE Microanalysis of Geological Material using the CSIRO Proton Microprobe, Nucl. Instrum. Methods Phys. Res., Sect. B, 1990, 47, 55–71.
10. C. A. Heinrich, W. L. Griffin, S. H. Sie, T. P. Mernagh, C. G. Ryan and D. R. Cousens, Quantitative PIXE microanalysis of fluid inclusions based on a layered yield model, Nucl. Instrum. Methods Phys. Res., Sect. B, 1991, 54, 291–296.
11. T. Shiraiwa and N. Fujino, Theoretical Calculation of Fluorescent X-Ray Intensities in Fluorescent X-Ray Spectrochemical Analysis, Jpn. J. Appl. Phys., 1966, 5, 886–899.
12. M. A. Reis and L. C. Alves, DATTPIXE, a Computer Package for TTPIXE Data Analysis, Nucl. Instrum. Methods Phys. Res., Sect. B, 1992, 68, 300–304.
13. M. A. Reis, L. C. Alves and A. P. Jesus, Matrix Effects Corrections for Quantitative TTPIXE Analysis, Nucl. Instrum. Methods Phys. Res., Sect. B, 1996, 109/110, 134–136.
14. D. K. G. de Boer, Calculation of X-Ray Fluorescence Intensities from Bulk and Multilayer Samples, X-Ray Spectrom., 1990, 19, 145–154.
15. C. Jeynes, V. V. Palitsin, G. W. Grime, C. Pascual-Izarra, A. Taborda, M. A. Reis and N. P. Barradas, External beam Total-IBA using DataFurnace, Nucl. Instrum. Methods Phys. Res., Sect. B, 2020, 481, 47–61.
16. M. A. Reis, P. C. Chaves, A. Taborda, J. P. Marques and N. P. Barradas, Fixed and free line ratio DT2 PIXE fitting and simulation package, Nucl. Instrum. Methods Phys. Res., Sect. B, 2014, 3189, 65–69.
17. P. C. Chaves and M. A. Reis, H⁺, O²⁺, O³⁺ and high resolution PIXE spectra of Yb₂O₃, Nucl. Instrum. Methods Phys. Res., Sect. B, 2017, 410, 193–199.
18. M. A. Reis and P. C. Chaves, DT2 fit of a high-resolution EDS PIXE Yb₂O₃ spectrum, X-Ray Spectrom., 2017, 46, 88–92.
19. M. A. Reis, G. R. Fonseca, P. C. Chaves, M. Herbst, M. Bühler, A. Fleischmann, A. G. Karydas and K. Phelan, High resolution energy dispersive spectrometry (HiREDS), a new tool for X-ray emission work, J. Anal. At. Spectrom., 2023, 38, 66–79.
20. D. B. Newell, F. Cabiati, J. Fischer, K. Fujii, S. G. Karshenboim, H. S. Margolis, E. de Mirandés, P. J. Mohr, F. Nez, K. Pachucki, T. J. Quinn, B. N. Taylor, M. Wang, B. M. Wood and Z. Zhang, The CODATA 2017 values of h, e, k, and NA for the revision of the SI, Metrologia, 2018, 55, L13–L16.
21. Table of Integrals, Series and Products, ed. I. S. Gradshteyn and I. M. Ryzhik, 44th edn prepared by Y. V. Geronimous, M. Y. Tseylin, translated from the Russian by Scripta Technica Inc., Alan Jeffrey Ed. edition, 1965.
22. Handbook of Mathematical Functions, ed. M. Abramowicz and I. Stegun, Dover Inc. Press, New York, 1st edn, 1965.
23. M. A. Reis, N. P. Barradas, P. C. Chaves and A. Taborda, PIXE analysis of multilayer targets, X-Ray Spectrom., 2011, 40, 153–156.
24. J. L. Campbell, N. I. Boyd, N. Grassi, P. Bonnick and J. A. Maxwell, The Guelph PIXE software package IV, Nucl. Instrum. Methods Phys. Res., Sect. B, 2010, 268, 3356–3363.
25. C. Pascual-Izarra, N. P. Barradas and M. A. Reis, LibCPIXE: a PIXE simulation open-source library for multilayered samples, Nucl. Instrum. Methods Phys. Res., Sect. B, 2006, 249, 820–822.

---