C Ankjaergaard

M Jain

P C Hansen

H B Nielsen

Towards Multi-Exponential Analysis in Optically Stimulated Luminescence

Keywords: Quartz, fitting, nonlinear least squares, Levenberg-Marquardt, Fredholm integral equation, first-order kinetics, time-resolved OSL (TR-OSL), Linearly Modulated

Parameters derived using both models are acceptable (statistically consistent and on an average within ~1% of the expected value) and no obvious preference is observed for any particular model, although there may be a suggestion that peak form data shows a smaller mean bias. This conclusion seems to be independent of the type of peak transformations investigated here. Furthermore, it is found that transformation of OSL decay data to a peak form gives better results than direct measurement of peak form data by, for example, varying the stimulation light intensity. The comparison of the two numerical methods suggests that the NLS method performs somewhat better than the FIE method; however, the latter has the advantage that it does not require the

INTRODUCTION

Luminescence is extensively used to determine absorbed dose (J/kg) in wide band gap insulators due to exposure to ionizing radiation. The dosimeters used are usually natural crystals, e.g., quartz and feldspars, or impurity doped artificially grown crystals such as lithium fluoride, aluminium oxide, etc. The irradiation of these crystals results in creation of free electrons and holes which are subsequently trapped in localised states (lattice defects) known as traps and recombination centres within a crystal, respectively. These trapping states thus store information about the absorbed dose, and the information can be read out in the form of luminescence by exposing the crystal to visible or near IR photons having sufficient energy to cause photoionisation of the occupied defects. Electron-hole recombination following photoionisation is the critical light generating step, and the resulting signal carrying dosimetric information is known as optically stimulated luminescence (OSL). Traditionally, OSL measurements use a constant, continuous flux of incident photons (CW-OSL); this results in a signal that usually shows monotonic decay with time. The OSL technique is widely used for estimation of absorbed dose from ionising radiation in wide ranging applications related to nuclear accidents, cosmic radiation in space, radiation facilities in health and power sectors, and in geochronology, i.e., dating sediments during the last 0.5 million years or so (see [START_REF] Bøtter-Jensen | Optically Stimulated Luminescence Dosimetry[END_REF] for an overview).

In addition to being of use for dosimetry, the OSL signal contains information about the distribution of traps and recombination centres in the crystal, and together with other techniques OSL signals can provide insights into the charge excitation, movement, and recombination in crystals. Such insights are vital to our understanding of the luminescence mechanism.

There are at least two instances where the OSL signal has a multi-exponential decay form:

1) CW-OSL of quartz: In case of natural quartz (the most commonly used material in accident dosimetry and sediment dating), for example, it is generally argued that the monotonic decaying CW-OSL emitted during exposure to blue light of constant intensity consists of several transients having an exponential form (Smith and Rhodes, 1994). The lifetimes of these transients are related to the photo-ionisation cross-sections (proportional to the decay constant) of the electron traps that participate in the luminescence process [START_REF] Bailey | Partial bleaching and the decay form characteristics of quartz OSL[END_REF]. In natural quartz from around the world, up to seven electron trapping states have been identified to participate in the OSL process [START_REF] Jain | Characterisation of blue-light stimulated luminescence components in different quartz samples: implications for dose measurement[END_REF]Singarayer and Bailey, 2003).

2) Pulsed OSL: An alternative method to stimulate a crystal is by pulsing the light intensity to obtain pulsed OSL. The signal measured during and between the light pulses is known as time-resolved OSL (TR-OSL). The processes that govern these signals are generally exponential in nature; in case of quartz the TR-OSL measured between the pulses is a decaying signal and consists of a dominant exponential transient having a lifetime between 30 and 45 µs and two relatively minor transients having relatively shorter and longer lifetimes than the main transient when using standard quartz dating procedure (preheating to 260 o C for 10 s and stimulating at 125 o C) [START_REF] Chithambo | Time-resolved luminescence of low sensitivity quartz from crystalline rocks[END_REF]Ankjaergaard et al., submitted;[START_REF] Pagonis | Experimental and modelling study of pulsed optically stimulated luminescence in quartz, marble and beta irradiated salt[END_REF]. These lifetimes reflect the relaxation of the excited state following electron holerecombination [START_REF] Chithambo | The analysis of time-resolved optically stimulated luminescence. I: Theoretical considerations[END_REF]Pagonis et al., submitted).

Thus, the experimental data acquired during either continuous stimulation (CW-OSL) or between the pulses in pulsed stimulation (TR-OSL) of quartz, and perhaps also other dosimeters, can be described by the model:

) ( ) exp( ) ( 1 decay t b t N t I n i i i i ξ λ λ + + - = ∑ = , 0 ≤ t ≤ T ,
Equation 1

where t is the time, N i is proportional to the initial population undergoing decay, and λ i is the decay constant for each of the n components in the signal. The constant b≥0 accounts for any background light and dark counts of the detector, and ξ(t) represents the noise in the data. The purpose of such multi-component analysis is to determine the constants n, N i and λ i from the measured data, as these quantities give information on the relative number density and a physical characteristic of the trap or centre such as photoionisation cross-section or recombination/relaxation lifetime, respectively.

In case of CW-OSL, the model (Equation 1) applies only in the case of negligible retrapping (first order kinetics); this is applicable, e.g., to the fast, medium and some of the slow components in quartz (see Jain and Lindvold, 2007a), but may not be applicable to all dosimeters. In TR-OSL the two dominant processes, i.e., the electronic transition from the conduction band to the recombination centre, and the relaxation of the excited state of the recombination centre, both follow first order kinetics. However, there are some materials that have non-first order processes, e.g., tunnelling, which will follow a different mathematical form. As discussed later a more general case can be defined using a spectral function for a continuous trap distribution; this can be reduced to a sum of n delta functions for discrete traps [START_REF] Istratov | Exponential analysis in physical phenomena[END_REF].

In addition to better understanding the OSL processes, it is also highly desirable that we can estimate the individual transients in the OSL for higher accuracy in retrospective dosimetry (see, e.g., [START_REF] Jain | A single-aliquot regenerative-dose method based on IR (1.49 eV) bleaching of the fast OSL component in quartz[END_REF]. This objective can be achieved by fitting some form of Equation 1 to the CW-OSL or the TR-OSL data. Unfortunately, fitting multiple exponential functions is an ill-posed problem (e.g., see an excellent review by [START_REF] Istratov | Exponential analysis in physical phenomena[END_REF]. [START_REF] Bulur | An alternative technique for optically stimulated luminescence (OSL) experiment[END_REF] suggested that fitting of multiexponential OSL data could be made more robust if the data had a peak form rather than a monotonic decay form. An OSL peak form could be obtained experimentally by linearly increasing the stimulation light intensity during the OSL measurement;

this technique was named as linearly modulated OSL (LM-OSL) [START_REF] Bulur | An alternative technique for optically stimulated luminescence (OSL) experiment[END_REF]. LM-OSL results in a signal that can be described using the single trap/centre model by the following model:

) ( ) 2 exp( ) ( 1 2 peak t b T t t T N t I n i i i i ξ λ λ + + - ⋅ = ∑ = , 0 ≤ t ≤ T , Equation 2
where T is the total measurement time, and N i , λ i , b, and ξ(t) are as described above. [START_REF] Bulur | A simple transformation for converting CW-OSL curves to LM-OSL curves[END_REF] and Poolton et al. (2003) have further shown that the same peak shaped signal could be achieved by transformation of the monotonically decaying multiexponential data (assuming first order kinetics). However, note that the transformation approach (to be discussed below) produces a linearly increasing background as a function of stimulation time t.

Although the shape of the LM-OSL curve is fundamentally different from the CW-OSL curve, the physical process causing them and the information contained in the two types of data are identical; the apparent differences arise because we view data in the time domain rather than in the event domain (Jain and Lindvold, 2007a). If one considers the probability of luminescence light produced per incident stimulation photon, and if there are two components X and Y (n=2) then

[ ] Φ - = ) ( exp 0 0 X Y Y X Y X a a a a σ σ Equation 3
where a X and a Y are the number density of the trapped electrons (cm -3 ) for the two different traps X and Y, a X0 and a Y0 are the initial concentrations, Φ defines the time integrated number of excitation photons per unit area (fluence in cm -2 ), and σ X and σ Y (cm 2 ) are the photoionisation cross-sections. For a given OSL measurement, , but the actual overlap of the signals will not differ for a given fluence (Jain and Lindvold 2007a). This problem has also been investigated by Wallinga et al. (2008) and Bos and Wallinga (2009a) in which several interesting examples of the data forms obtained by changing incident photon flux were studied. Nonetheless, it was concluded that the separation of the components at any given time was the same for all the different stimulation modes. One method of increasing the physical separation of the transients is differentiation of the OSL signal (Bos and Wallinga, 2009b); however, this approach is not suited to low amplitude, slowly decaying components because of low signal-to-noise ratio. Since the nature of multi-exponential decay and multi exponential peak data (in other than differentiation based methods) are fundamentally identical, it has been discussed in the literature whether peak-form data has in fact any advantage over the decay-form data for exponential analysis [START_REF] Huntley | Thoughts arising from "Choi, Duller and Wintle: Analysis of quartz LM-OSL curves[END_REF][START_REF] Huntley | Response to Jain and Lindvold[END_REF]Jain and Lindvold, 2007a, b;Wallinga et al. 2008;Bos and Wallinga, 2009a;Bos and Wallinga, 2009b). However, this question has not been thoroughly investigated to date.

The aim of this work is to carry out a comparative analysis of the performances of the two above-mentioned data collection systems, decay and peak form, in predicting the 'true' trap parameters. This analysis is based on artificially generated multiexponential TR-OSL data with seven parameters. We use two fundamentally different numerical approaches which have previously been used for fitting luminescence signals by different workers viz. a) multi-exponential analysis using a nonlinear least squares (NLS) formulation solved by a Levenberg-Marquardt approach (see, e.g., [START_REF] Jain | Characterisation of blue-light stimulated luminescence components in different quartz samples: implications for dose measurement[END_REF], and b) the spectroscopic analysis using a first-kind Fredholm integral equation (FIE) [START_REF] Agersnap | Dosimetry based on thermally and optically stimulated luminescence[END_REF]Whitley and McKeever, 2001; also see the review by [START_REF] Istratov | Exponential analysis in physical phenomena[END_REF] with a view to derive a general protocol for multiexponential analysis of three discrete OSL components. The NLS method is a standard method for fitting data and is implemented in most mathematical and statistical software packages. This method requires that the number of exponential components is specified beforehand by the user. The FIE method is a radically different approach compared to NLS as this method automatically determines the number (or spectrum) of components without prior user input. These approaches are briefly described in Section 5. We emphasize that our goal is (primarily via simulation studies) to compare the above-mentioned methods for analysing OSL and TR-OSL data, with a focus on the robustness of the estimation of the parameters N i and λ i in the underlying OSL models. We note that related work in Pereyra and Scherer (2009) focused on the fitting (or prediction) abilities of the two methods based on the properties of the residuals; we note that this is an altogether different mathematical problem and leads to a different conclusion than in the present study focussed on the precision and accuracy of the derived parameters.

INSTRUMENTATION AND METHODS

Much of this paper uses artificial data. However, laboratory measurements on a quartz sample are used to test the methods and to generate realistic artificial data for exponential analysis. These measurements were undertaken using a colluvium quartz sample from Tanzania, sample: 963602 consisting of 150-300 µm grains. The grains were extracted from the bulk sample by sieving, heavy liquid separation, and HF treatment as described in Wintle, 1997. Sample measurements were carried out on a

Risø TL/OSL-20 reader with an integrated pulsing option to control the stimulation LEDs, and with a photon timer attachment with detection resolution of 100 ps to record the TR-OSL data [START_REF] Lapp | Development of pulsed stimulation and photon timer attachments to the Risø TL/OSL reader[END_REF]. The apparent bin-resolution can be redisplayed after data collection, as it is specified as 2 n ×100 ps. Blue light stimulation was performed with a 470 ± 30 nm LED array delivering 50 mW.cm -2 at the sample, and a 7.5 mm thick Hoya U340 filter detecting between 340-350 nm was inserted beneath the photomultiplier tube (PMT).

The TR-OSL data are recorded by time-stamping photons emitted from beta irradiated and preheated (260°C for 10 s) samples between the stimulation light pulses at 125 o C.

The net signal is then derived by summing several hundred thousands of these light pulses emitted from the sample. The duty cycle consisted of a 50 µs on-time and a 500 µs off-time. The data are then summed up to be presented at the resolution of 0.4 µs.

In addition to the OSL counts, the PMT will register counts that constitute a background to the OSL signal. Background counts are generated by thermal excitation of electrons in the photocathode, system electronics, and interaction of the detector with cosmic radiation, which together constitute the 'dark counts' of the PMT, and any scattered or stray light from the stimulation source. As it is important to include these in the artificial data model, the background counts were also measured in the same way as above but using an annealed quartz sample (700°C for 5 min to remove nearly all trapped charge) instead of a beta irradiated sample. Note that in pulsed OSL measurements, the stimulation light is switched off, therefore, background would be equivalent to the 'dark counts' of the detector.

DATA SIMULATION FOR DECAY-FORM DATA

Artificial TR-OSL data was created in order to compare the performance of different methods, and the data comprises two contributions -light from the sample and the background counts.

In order to generate artificial data it is first important to understand the PMT counting statistics for these two processes; this involves measurement of a constant count rate by the detector and determination of the mean and spread in the distribution of the counts. Since the TR-OSL is a decaying signal it cannot as such be used for this type of analysis. Therefore, we used a constant light source; blue LEDs operated with attenuated detection (a neutral density filter inserted underneath the PMT instead of the U340 filter) is used to do the statistical analysis applicable to signal counts. We generated a similar count rate as the sample, was used to mimic distribution of luminescence counts. Data was collected during a 500 µs on-time. The background was measured in the same way as TR-OSL but on an annealed sample. Poisson distributed as the dark noise is also contained in this data, and so the condition to test whether the distribution is a single Poisson distribution is in fact false.

For practical purposes, as can be seen in the figure, the data approximates a Poisson distribution, and so we will assume that the (TR-OSL) signal follows a Poisson distribution.

The values of the decay constants and amplitudes for generating the artificial data analogous to the TR-OSL of quartz in the 380 nm emission is based on experimental values found by [START_REF] Chithambo | Time-resolved luminescence of low sensitivity quartz from crystalline rocks[END_REF]. The artificial data is chosen to contain three exponential transients following Equation 1 with amplitudes

A i = N i •λ i (cts/0.408 µs)
and decay constants λ i (µs -1 ) of A 1 = 7000, λ 1 = 0.500, A 2 = 20000, λ 2 = 0.0300, A 3 = 3000, and λ 3 = 0.0125 (the decay constants correspond to lifetimes, τ of 2, 33 and 80 µs). For the analysis, the lifetimes are converted into their reciprocal values, i.e., the decay constants used in Equation 1. The artificial data sets are generated to contain m = 2000 points with a bin width of 0.408 µs giving a decay length between the light pulses of 816 µs. From Equation 1we can express the exact data through the model:

∑ = - = 3 1 ) exp( ) ( i i i t A t I λ Equation 4
The individual points Y j in the artificial data, with background, can therefore be written as: The relative level of background in an artificial TR-OSL curve is estimated by the initial-signal to background ratio (ISBR) defined as the ratio between the maximum signal amplitude at t = 0 and the constant background b, i.e., I(0)/b. For this study, ten different ISBR in the range 6-6000 were used to cover for the variations in the OSL sensitivities of different quartz samples. For each ISBR ratio we generated five instances of the same decaying TR-OSL signal but with different noise realisations.

Y j = I j + b j j = 1,
Thereby for each ISBR level, more realistic error estimates of the parameters were obtained than those estimated by the fitting algorithm. In total, 50 different instances of artificial decay data were produced. Before fitting the artificial data, an average background value _ b found from the average over the last 10 data points is subtracted in order to remove the constant background component, such that the final artificial data is:

_ b Y y j j - = Equation 6
This is standard procedure for analysis of measured TR-OSL data. As long as the background is truly constant (for times t greater than six times the lifetime of the slowest component), the number of data points used in the average background estimation is not critical. We used 10 points for reasonable statistics.

Figure 1(c) shows an artificial TR-OSL curve from Equation 6with an ISBR of 600 (black line) together with the 'exact' curve from the model in Equation 4(grey line).

The curves were calculated using the parameters defined above. The insets to Fig.

DATA SIMULATION FOR PEAK-FORM DATA

Since TR-OSL is always a decaying signal, irrespective of the manner in which the sample was excited during a pulse, it is not possible to directly obtain a peak form signal such as in the LM-OSL. However the decaying TR-OSL can be converted into a peak form such as LM-OSL either by varying the integration intervals (e.g. Poolton et al., 2003) or by variable rescaling of the intensity and time axis (e.g. [START_REF] Bulur | A simple transformation for converting CW-OSL curves to LM-OSL curves[END_REF].

In this work we transform the decay-form data y j , generated as described above, into two different types of peak forms: (i) a simple variable rescaling of the intensity (simple transformation) and (ii) a variable rescaling of the intensity and time space (pseudo LM-OSL of [START_REF] Bulur | A simple transformation for converting CW-OSL curves to LM-OSL curves[END_REF]). Note, however, that in the context of TR-OSL off -time signal, the term LM-OSL does not have any physical connotation; it merely serves as a different mathematical transformation. The two transformations correspond to the two models:

Simple transformation:

∑ = + ⋅ = 3 1 i 2 / 1 ) - exp( ) 2 ( ) ( i i b t A tT t I λ Equation 7 Pseudo LM-OSL: b T u u T A u I i i i + - ⋅ = ∑ = 3 1 2 ) 2 exp( ) ( λ Equation 8
Here we define T = 2t TR-OSL,max , twice the length of the TR-OSL decay length (1632 µs), and we introduce a new independent variable

2 / 1 ) 2 ( tT u = . Moreover, A i = N i •λ i ,
where N i is proportional to the initial population undergoing decay, λ i is the decay constant for each of the n components in the signal, and b≥0 is a term to account for any background light and dark counts of the detector. The two types of peak data thus correspond to the transformations:

y j → (2 t j T) 1/2 y j , t j unchanged , Equation 9
and y j → (2 t j / T) 1/2 y j , t j → u j (by u j = (2t j T) 1/2 ), Equation 10respectively, of the artificial decay data after correction for background. The artificial decay data and the simple transformed peak data share the same linear time domain t, in contrast to the pseudo LM-OSL peak data which is a function of a 'nonlinear' time domain u.

In Fig. 2, the artificial TR-OSL data from Equation 6with an ISBR ratio of 15 is shown together with the same curve after the simple peak transformation in Equation 9, and after the peak transformation in Equation 10 to obtain pseudo LM-OSL. The pseudo LM-OSL curve is 'stretched' compared to the simple peak transformation due to the time domain transformation and is plotted on the top axis as a function of u.

The TR-OSL and the simple peak transformation is plotted on the bottom axis, t.

NUMERICAL METHODS

Two different numerical methods are used for parameter estimation. Both methods are least squares methods in the sense that they minimize the sum of squared residuals between the data and the model. One method is a standard nonlinear least squares (NLS) method, which is perhaps the most used method for exponential fitting. As an alternative, we also consider a method based on a Fredholm integral equation (FIE)

approach, which has the advantage that it does not require the user to select the number of components in the model. To make a fair comparison between the two methods, we do not use any weighting in the least squares problems.

For the decay form in Equation 1, the nonlinear least squares (NLS) method fits the model

∑ = - = n i i i t A t I 1 ) exp( ) (
λ to the m data points y j after the background has been subtracted:

∑ ∑ = = - - m j n i j i i j A t A y i i 1 2 1 , ) ) exp( ( min λ λ Equation 11
Here we have introduced the amplitudes A i = N i •λ i , and t 1 ,…,t m denote the measurement times. A similar expression can be written for the simple peak transformation in Equation 7. For the peak form in Equation 2, the model

∑ = - ⋅ = n i i i T u u T A u I 1 2 ) 2 exp( ) (
λ is fitted to the data points:

∑ ∑ = = ⋅ - m j n i j i j i j A T u u T A y i i 1 2 1 2 , ) ) 2 exp( ( min λ λ Equation 12
where u j is the transformed time from equation 10. The underlying numerical method is known as the variable projection method, which is a Levenberg-Marquardt algorithm tailored to the specific problem; see [START_REF] Golub | Separable nonlinear least squares: the variable projection method and its applications[END_REF] and Pereyra and Scherer (2009) for details. This algorithm is built into the software package SigmaPlot, which was used in this study. Since the data is fitted to a specific model, the user must specify the number of components n in this model prior to fitting the data. We impose the positivity constraints λ i > 0 on the decay constants, but did not enforce positive amplitudes as negative values never were encountered.

In the second method, a spectroscopic analysis approach is used (see, e.g., [START_REF] Agersnap | Dosimetry based on thermally and optically stimulated luminescence[END_REF][START_REF] Istratov | Exponential analysis in physical phenomena[END_REF]Whitley and McKeever, 2001), in which λ is considered as a continuous variable and a(λ) ≥ 0 is an unknown 'amplitude density function' which identifies the decay constants (and thereby the components) present in the data y(t):

) ( ) ( ) exp( max min t y d a t = - ∫ λ λ λ λ λ Equation 13
This is illustrated in Fig. 3 for a set of artificial data y(t), shown inset, containing three components. The areas under the peaks of a(λ) are related to the amplitudes of the decay constants. The similar equation for the peak form is:

) ( ) ( ) 2 exp( 2 max min u y d a T u T u = - ∫ λ λ λ λ λ Equation 14
Equations 13 and 14 have the form of a Fredholm integral equation (FIE) of the first kind (see [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF]. Unlike in the NLS method, the user does not have to specify the number of components; the amplitude density function a(λ) determines this number (see below).

To solve the FIE in Equation 13numerically, it is discretized by means of a quadrature method (see [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF] which leads to the linear least squares problem: is the case then it is also expected that the amplitudes are approximately equal to the area under each of the spikes, i.e.,

∑ ∑ = = - Δ - m j p k k j k k j a t y 1 1 2 ) ) exp( (
∫ Ω ≈ i d a A i λ λ) (
where Ω i denotes a small interval around λ i . For each spike we therefore compute the corresponding amplitude by means of

∑ ∈ Δ = i k k k i a A spike λ .
Equation 16

More than n spikes were occasionally encountered (although a situation with less than n spikes never occurred), and thus a 'rejection criterion' was needed to remove spurious spikes. It is noted that the contribution to the total signal from the ith component can be measured by the integral: 17Spikes with ratios a i /λ smaller than a fixed threshold value determined by the area under the data and the noise level in the data are therefore rejected. In this paper, a threshold of 0.01 was used, i.e., components whose contribution is less than 1% of the total signal are rejected.

i i i i i T i i a T a dt t a λ λ λ λ ≈ - - = - ∫ )) exp( 1 ( ) exp( 0 Equation

ARTIFICIAL DATA RESULTS

In this section the two fitting methods are used for analysing the artificial data sets for both decay and peak forms.

Nonlinear Least Squares Method: Decay vs. Peak Form

As described in Section 5, the nonlinear least squares algorithm, as implemented in SigmaPlot, is not able to identify the number of components present in the data, and therefore a model with a specific set of components must be chosen by the user prior to fitting the data under examination. To eliminate the advantage that the number of components is known in the artificial data sets, a systematic approach was used for each set of data to identify the number of components as if it was unknown. Each data set was fitted with first two, then three, and finally with four (and so on) components to identify the correct number of components present in the data set. This is illustrated in Fig. 4(a)-(c) for an artificial data set with an ISBR of 600 and fitting using a model with 2, 3, and 4 components, respectively, with corresponding parameter values listed in Table 1. Fitting with too few components caused a systematic behaviour in the residuals (Fig. 4(d)) which is undesirable, while the residuals from the three- 1. The fourth fitted component merely splits the third component into two parts and, as a result, the residuals are identical in the 3 and 4 component fits. This was tested for several curves and it is concluded that the condition for the best fit is the minimum number of components for fits with identical residuals.

In Table 1, the A i and λ i parameters in the first three columns are listed with internal uncertainties estimated by SigmaPlot. To test whether these uncertainties are large enough to contain the variations in parameter values caused by the noise, five datasets with an ISBR of 600 were fitted to calculate the mean and standard error on the parameters. These values are presented in the last column of Table 1. The calculated standard errors based on five different realisations of the noisy curve are systematically bigger than those estimated by SigmaPlot for any individual curve. We therefore fit five datasets for each of the different ISBR values as presented in Section 3 to get more realistic estimates of uncertainties on the parameters from the fits.

Using the built-in nonlinear least squares algorithm in SigmaPlot, we fitted the 50 decay formed curves, the 50 simple peak transformation data sets, and the 50 pseudo LM-OSL peak data sets. The results are shown in Fig. 5. The computed amplitudes A 1 , A 2 , and A 3 are compared in (a)-(c), for the 10 different ISBR. Likewise, the decay constants λ 1 , λ 2 , and λ 3 are compared in Fig. 5 (d)-(f). The dashed lines in each plot represent the 'true' values used in creating the artificial data. It is first observed that all amplitudes and decay constants using the two peak forms are identical; the NLS method does not seem to be able to distinguish the two types of data from each other.

Secondly, the results from decay and peak form data are consistent with each other and within the 'true' value. The uncertainties on the parameters derived for fitting peak data are generally smaller than those for the decay form data. For very high ISBR there is a slight tendency for underestimation of the parameters from the peakform data fits, see Fig. 5 (c), (e), and (f).

Spectroscopic Method: Decay vs. Peak Form

The Fredholm integral equation model was used to fit the same data sets from the previous section. It was observed that the number of components (n = 3) present in a data set is often overestimated by one, two, or three components; but with the rejection criterion from Section 5 almost all of these extra components are rejected by our method. This is shown in Table 2, which lists the number of extra components found in the 50 data sets before and after the rejection criterion (RC) has been used.

The two peak forms had identical numbers of components present before and after use of the rejection criteria, and are therefore only described in one column. Out of the 50 decay-form data sets, 38 had too many components and this number got reduced to 5 after using the rejection criterion. Similarly, from the 50 peak-form data sets, 42 had too many components and this number was reduced to 3.

Figure 6 shows the results from fitting with the Fredholm integral equation model, and again the amplitudes are compared in (a)-(c) while the decay constants are compared in (d)-(f). Similar to Fig. 5, the points in Fig. 6 are averages over the parameters from five data sets with the same ISBR. As described above, the results from five decay-form data sets and three peak-form data sets were rejected due to the overestimation of the number of components (see Table 2). Again we note that the parameters from the two peak forms are identical. The parameter estimates generally have larger uncertainties using this fitting method compared to those using NLS.

Furthermore, there is a tendency to underestimate the parameters in Fig. 6(c), (e), and (f), and to overestimate the parameters in Fig. 6(b) for high ISBR levels for both decay and peak shaped data. For all three data forms, several of the parameters are not consistent with the 'true' parameter values at high ISBR levels.

COMPARISON OF THE METHODS

To make an overall assessment of the performance of the two fitting methods and the two types of data, the parameters from the fitted decay-and peak-form data were analyzed by calculating the mean and the standard deviation of every parameter divided by its expected value. The optimum mean value should therefore be 1. The values are given in Table 3; columns 2 and 3 list the parameters found using the FIE method for decay-and peak-form data respectively, whereas columns 4 and 5 list the parameters found using the NLS method for decay-and peak-form. Note that the values given in rows 3 and 5 are valid for both the simple peak transformation and for the pseudo LM-OSL transformation, as these parameter values were identical (see Fig. 5 and6). The number of data sets included in the calculation of the mean and standard deviations is given in row 2; for the FIE method this number is not 50 due to a few instances where the rejection criteria did not remove all over fitted components (see Table 2 andSection 6.1).

If the decay and peak parameter results from the FIE method (columns 2 and 3) are compared, four peak parameters (A 1 , A 3 , λ 2 , and λ 3 ) are slightly better estimated, i.e., have a smaller mean bias than the decay parameters, but the overall spread in the individual parameter estimates (standard deviation) is larger for peak fitting compared to the decay curve fitting. If the decay and peak parameters from the NLS method (columns 4 and 5) are compared instead, then the three peak parameters (A 2 , A 3 , and λ 3 ) have smaller mean bias compared to the decay parameters, and the value of λ 1 is identical for the two data types. Similarly the standard deviations are smallest for the four peak parameters (A 2 , A 3 , λ 2 , and λ 3 ) compared to the decay parameters. There is an overall greater concordance in the mean values from the two models for the NLS method compared to the FIE method. In general both the models perform very well and the mean values are statistically indistinguishable from the expected values if one considers the measured spread in parameter estimates. However, even though there is no clear preference between the two models, our data may hint in favour of the peak model if one only considers the mean bias.

When comparing the two fitting methods, FIE and NLS, on the basis of the six parameters each estimated for each of the two models (decay and peak), the NLS fitted 10 out of 12 parameters more accurately, and 8 out of 12 standard deviations were the smallest compared to the FIE values. Although the FIE method is more objective, the numbers in Table 3 suggests that the NLS method estimates the parameters more accurately and with smaller error.

'TRUE' LM-OSL VS. PSEUDO LM-OSL

We now investigate the relative performance of the peak model for two cases 1) pseudo LM-OSL data derived from transformation of decay form (CW-OSL or TR-OSL) data, and 2) 'true' LM-OSL. There could be differences in the performance of the two data sets as the signals are measured with different signal to background ratios, and moreover, the transformation step in case 1) also involves transformation of the noise. Thus, it is not apparent whether a true peak measurement such as LM-OSL is preferable compared to a pseudo-LM-OSL or otherwise. Note that a true peak measurement is only possible in the conventional OSL where varying the stimulation light intensity gives rise to a peak shaped signal (see Section 4). This is not an option in TR-OSL measurements where the signal always decays with time between the two pulses. In this section we will compare the results from artificial 'true' LM-OSL with the previous results obtained using pseudo LM-OSL.

We create artificial 'true,' i.e., experimentally obtained LM-OSL containing three components with the same parameter values as in Section 4 and including a background b≥0 arising from any dark counts of the detector (possible minor time dependent background from stray stimulation light is neglected in order to simplify the analysis), and a noise component ξ(u) that reflects the Poisson noise in the data:

∑ = + + - ⋅ = 3 1 2 ) ( ) 2 exp( ) ( i i i u b T u u T A u I ξ λ Equation 18
where T is twice the length of the OSL decay length of 1632 µs [START_REF] Bulur | An alternative technique for optically stimulated luminescence (OSL) experiment[END_REF]. Prior to fitting, an average background _ b estimated from the last 10 data points is removed.

As in the previous sections 50 'true' LM-OSL curves are produced, 5 curves for each of the 10 ISBR levels. Fig. 7(a) shows the true LM-OSL (grey curve) and the pseudo LM-OSL (black curve) for the highest ISBR level, 6000 and in Fig. 7(b) the same curves but for the lowest ISBR of 6. In both figures, the 'true' LM-OSL curves appear much noisier than the pseudo LM-OSL and therefore one might presume that the pseudo LM-OSL would be a better choice for fitting.

Using the NLS method the 50 new curves are fitted, and the parameters are shown in Fig. 8 (black circles) together with the values found from fitting the pseudo LM-OSL in Section 6.1 (grey circles). The computed amplitudes A 1 , A 2 , and A 3 are compared in (a)-(c), for the 10 different ISBR. Likewise, the decay constants λ 1 , λ 2 , and λ 3 are compared in Fig. 8(d)-(f). For the low background levels (ISBR<100), the true LM-OSL generally have very large error bars reflecting bad fits and large scatter in the parameter values within each level, and several parameter values are inconsistent with the known true parameter values. It is clear that the noise in the data for ISBR<100 results in very poor fits and for noise levels of this magnitude, pseudo LM-OSL is clearly preferable. To make an overall assessment of the accuracy of the parameters compared to those of the pseudo LM-OSL data (column 5, Table 3), the mean and standard deviation was calculated for each of the six parameters and listed in column 6 in Table 3. Although the mean values of the 'true' LM-OSL parameters are either very similar or better than the pseudo LM-OSL estimated parameters, the standard deviation is bigger for all the parameters, reflecting the much larger scatter found in the parameter estimates.

It is therefore concluded, that using pseudo LM-OSL data results in better fits irrespectively of the noise level in the data compared to 'true' LM-OSL data. This conclusion assumes that the conditions of measuring the OSL decay curve or the LM-OSL data are the same i.e. the same instrument is used and thereby the same light stimulation source. That pseudo LM-OSL is better than 'true' LM-OSL is probably because the signal-to-noise ratio for the decay curve measurement (CW-OSL) is always much higher than that for the corresponding true LM-OSL measurement.

Therefore, the pseudo LM-OSL derived by transforming CW-OSL data is less noisy than the 'true' peak measurement. Moreover, there may be a problem of time varying background arising form the stimulation light in the true LM-OSL measurement. This problem can further complicate the exponential analysis of the true LM-OSL data.

PERFORMANCE OF THE NLS METHOD ON MEASURED QUARTZ TR-OSL DATA

TR-OSL data was measured from the Tanzanian quartz sample for a regenerative irradiation dose of 30 Gy. Furthermore, this dose was repeated but now measured with a neutral density filter inserted beneath the PMT to attenuate the light from the sample by approximately a factor of 16 and thereby obtaining a smaller ISBR. Both curves were transformed into peak-shape according to the model in Equation 8and are shown in Fig. 9(a). The curves are shown on a log intensity scale in the inset.

The data was initially fitted in SigmaPlot with one component, and the corresponding residuals are shown as the grey curve (30 Gy) and black curve (30 Gy, attenuated) in Fig. 9(b), upper and lower graph respectively; the grey residuals shows a strong trend indicating that the data is not fitted adequately, whereas the black curve show a much less visible trend. As the two curves are from the same sample, they must contain the same number of components, and so two components are sufficient to fit the data. The residuals from the two two-component fits are shown in Fig. 9(c), the upper grey curve being the 30 Gy residuals and the lower curve the attenuated 30 Gy residuals.

The trend in the grey residuals has been reduced significantly and the black residuals have been improved slightly indicating that only two components are present in the measured TR-OSL data. Finally the curves were fitted with three components and the residuals from these two fits are shown in Fig. 9(d); the residuals have not improved to the eye, and the third component is clearly in excess as the decay constants are identical in both cases to an already identified component and the amplitudes have been halved. It is therefore clear that only two components are present in the data. The decay constants and amplitudes found from the two-component fits are given in Table 4. The 30 Gy curve fit resulted in two components with lifetimes of 39 µs and 103 µs for amplitudes of 36,959 cts.(0.4 µs) -1 and 685 cts.(0.4 µs) -1 . The attenuated 30 Gy curve fit gave the lifetimes 38 µs and 107 µs for amplitudes 2,219 cts.(0.4 µs) -1 and 65 cts.(0.4 µs) -1 .

In [START_REF] Ankjaergaard | Optimising the separation of quartz and feldspar optically stimulated luminescence using pulsed excitation[END_REF], it was shown that for a wide selection of quartz samples, the main dominant commonly observed component has a lifetime of 37 ± 5 µs found from the mean and standard deviation from the 30 samples investigated. Furthermore they found that in some samples, a longer or a shorter lifetime component was present with lifetimes of 100 ± 40 µs and 6 ± 1 µs. The lifetimes found here are in good agreement with each other despite different noise levels, and are in very good agreement with published values in [START_REF] Ankjaergaard | Optimising the separation of quartz and feldspar optically stimulated luminescence using pulsed excitation[END_REF], [START_REF] Chithambo | A pulsed light-emitting-diode system for stimulation of luminescence[END_REF], and [START_REF] Chithambo | Time-resolved luminescence of low sensitivity quartz from crystalline rocks[END_REF]. Given the above amplitudes, the ISBR ratios for the individual component for the attenuated curve were calculated to be 44 and 1.3. Thus it is possible to estimate the parameters A 2 , λ 2 for a second component even though the initial intensity of this component has a comparable size to the background.

CONCLUSIONS

A thorough study of two different models for multi-exponential OSL data, the decay form and the peak form, has been undertaken to investigate which data form is better suited for fitting noisy data. Each of these data types were analysed using two fundamentally different numerical methods, one based on solving the nonlinear least squares (NLS) data fitting problem by a Levenberg-Marquardt algorithm, and the other based on solving a linear first-kind Fredholm integral equation (FIE) by means of a non-negativity constrained linear least squares problem. We conclude the following:

1) Parameters derived using both models are statistically consistent with the expected value, given the spread in the data. The mean values themselves are on an average within ~1% of the expected values for both the models. Of the two data types, the decay form and the peak form, neither is clearly superior over the other in estimating parameters when fitting with the FIE or the NLS method, although there seems to be a slight tendency for the peak form to show smaller mean bias in our limited dataset. Furthermore, the two numerical methods are insensitive to whether the peak shape is obtained by only variable scaling of the intensity, or by variable scaling of both intensity and the time space (e.g., pseudo LM-OSL) 2) Although both numerical methods tested here perform well, there is a suggestion that the NLS method is marginally better and relatively less sensitive to the noise in the data. However, as demonstrated here, if an appropriate rejection criterion is incorporated into the routine, then the FIE method has an advantage that it does not require user to define the number of components. The computation time was similar for the two methods.

3) For the continuous-wave OSL (CW-OSL) measurements we conclude that the fitting of peak transformed data (i.e., pseudo LM-OSL) is more robust than the fitting of 'true' LM-OSL data. This is probably because, for the same parameters, the CW-OSL can be measured with a higher signal-to-noise ratio than in the corresponding LM-OSL. It should, however, be noted that CW-OSL data should be collected with high enough frequency so as to be able to resolve the fastest decaying component in the signal to obtain the best fit of multi-exponential OSL data.

Based on the above observations, a possible recommendation for analysing multiexponential data could be fitting peak-transformed decaying OSL signal using the NLS approach. This algorithm is already implemented in many commercial types of software such as SigmaPlot (used here).

Table 2 Number of OSL curves with over-fitting by either one, two, or three exponential components, when using the Fredholm integral equation model (Section 5). The numbers are given before and after the use of the rejection criterion (RC).

Table 3 Summary of the mean (% mean bias) and standard deviations for the parameters presented in Figs. 5,6, and 8 based on data for all the ISBR values for the FIE method using decay (column 2) and peak (column 3) data, and for the NLS method using decay (column 4) and peak (column 5) data. The % mean bias has been defined as (expected value -measured value / expected value)*100. Columns 6 lists the values for the 'true' LM-OSL parameters, and these are to be compared to the values pseudo-LM-OSL in column 5. For each of the six parameters, the mean and standard deviation were calculated using the number of data sets given in row 1, for the FIE method, this number is less than 50 due to imperfect removal of components by the rejection criteria (see Table 2). Note that the values in columns 3 and 5 represent both the simple peak transformation parameters as well as the pseudo LM-OSL parameters as these are identical.

Table 4 Amplitude and decay constant parameters from fitting quartz TR-OSL from a 30 Gy dose (columns 1-2) and TR-OSL from a 30 Gy dose and attenuating the signal by a factor of 16 (columns 3-4). Figure 2 Artificial decay data from Equation 6with an ISBR of 15 (black line), together with the same data transformed into peak form by using the simple transformation in Equation 9 (grey line) and transformed into pseudo LM-OSL using Equation 10(dotted black line). Note that the decay data and the peak data belong to different abscissa axes, and that the decay data and the simple transformed data share the same linear time domain t, while the pseudo LM-OSL data has a nonlinear time domain u.

FIGURE CAPTIONS

Figure 3 A plot of the amplitude density function a(λ) given in Equation 13 vs. the parameter variable, λ, for a set of artificial data y(t) with three components. The artificial data y j and its three components ) exp(

j i i t A λ -
are shown inset, see Equation 6. Figure 9 TR-OSL data measured from a colluvium quartz sample from Tanzania. The sample was given a dose of 30 Gy and was heated to 260 o C for 10 s prior to blue LED stimulation at 125 o C for 22 s. The on-time was 50 µs and the off-time was 500 µs.

Prior to fitting, the data was transformed into a peak shape using Equation 7(grey curve). The experiment was then repeated, but with a neutral density filter inserted

  therefore, for a given photon fluence, the ratio of the light output from any two components will be constant, no matter how that fluence was achieved in time. In the time domain, however, different shapes of luminescence intensity can be achieved by varying the fluence-rate or flux during the measurement

  Figure 1(a) shows the distribution of the background noise for a data resolution set to ~0.4 µs, and a Poisson distribution function with a mean of 5.7 cts.(0.4 µs) -1 (derived from integration of 40,000 pulse stimulations) has been fitted to the data (black line). A Matlab Chi-square goodness-of-fit test could not reject the hypothesis that the data is Poisson distributed with the probability of observing the given result: p = 0.9957. Therefore, in our artificial data model, we can assume that our instrumental background noise is Poisson distributed. The distribution of counts from the LEDs is shown in Fig. 1(b), with the same resolution as in 2(a), where it has been fitted to a Poisson distribution with a mean value of 14.8 cts.(0.4 µs) -1 shown as the black line. The Chi-square goodness-of-fit test rejected the hypothesis that the diode noise data is Poisson distributed with a probability of p = 2.21•10 -9 . It is very likely that this distribution is not perfectly

  ~ P(I (t j )) and b j ~ P(b) are numbers randomly selected from Poisson distributions with mean values I (t j ) and b, respectively (recall that b is the background).

  1(c) zoom in on the early and later parts of these curves. InFig. 1(d) the deviation between the 'exact' curve I(t) and the noisy and background-corrected data y j (both from (c)) are shown. As expected from the Poisson statistics the random deviation increases with the magnitude of the data. More importantly, the errors do not show any structure and our model fitting should at best be able to mimic this.

  λ k are the quadrature points, k λ Δ denotes the width of the interval that contains the kth quadrature point, and a k are the approximations to a(λ) at the quadrature points, p = 500 quadrature points, and due to the large range of λvalues needed in OSL problems, the quadrature points are distributed logarithmically between λ min = 0.001 and λ max = 100. The same technique is used to discretize Equation14. No weighting of the data is used and non-negativity constraints are enforced on the solution which has the side effect of causing the spikes in the solution to be narrower than without constraints. The constrained problem is solved by means of the NNLS algorithm(Lawson and Hanson, 1995, Chapter 23) implemented in Matlab. It is beyond the scope of this paper to go into further details of the computational algorithms, and we refer to the above-mentioned books for details.When the data adheres to the model in either Equation 13 or 14, it is expected that the nonnegative function a(λ) has narrow spikes for i λ λ ≈ , seeFig. 3. Moreover, if this

  component fit (Fig 4(e)) are evenly spread around zero and they do not visibly behave systematically, which according to Fig. 1(d) is the desired result. When fitting with four components, and thereby forcing the algorithm to find more components than what is truly present in the data, either of two situations occur: (i) two components have identical decay constants or (ii) the algorithm creates a fourth component with a decay constant very similar to an already existing decay constant and changes the remaining three decay constants slightly. Both of these possibilities resulted in no further visible improvement of the residuals, see Fig. 4(f). The redundancy in the four-component fit becomes clear by looking at the second (3 component fit parameters) and third column (4 component fit parameters) in Table

Figure 1 (

 1 Figure 1 (a) Distribution of dark counts measured during the off-time in 0.4 µs bins. The data was measured for 22 s with an on-time of 50 µs and an off-time of 500 µs. The noise distribution has been fitted to a Poisson distribution with a mean value of 5.7 cts⋅(0.4 µs) -1 shown as the black line. (b) The distribution of photon counts from the blue LEDs. A Neutral Density filter was inserted underneath the diodes to reduce the light level significantly and data was collected for 22 s during an on-time of 500 µs. The

Figure 4

 4 Figure 4 Artificial decay data with an ISBR of 600 are shown together with (a) two, (b) three, and (c) four component fits to the data. The insets show the same on the log intensity scale. (d)-(f) show the corresponding residuals from the fits.

Figure 5 Figure 6

 56 Figure 5 Results from fitting decay-form data (grey filled circles), simple transformed peak data (open black circles), and pseudo LM-OSL data (black diamonds) using the NLS method. In each plot the dashed line indicates the known parameter value, and each

Figure 7 Pseudo

 7 Figure 7 Pseudo LM-OSL (black curve) and 'true' LM-OSL (grey curve) as calculated from Equations 8 and 18 respectively, for (a) the lowest noise-level with an ISBR of 6000 and (b) and the highest noise-level with an ISBR of 6. The parameters are the same as given in the caption of Fig. 1.

Figure 8

 8 Figure 8Results from fitting pseudo LM-OSL data (grey filled circles) and 'true' LM-OSL data (filled black circles) using the NLS method. In each plot the dashed line indicates the known parameter value, and each point is the average over five parameter values for the same ISBR. The left column shows the results for the amplitudes (a) A 1 = 7 000 cts.(0.4 µs) -1 , (b) A 2 = 20 000 cts.(0.4 µs) -1 , and (c) A 3 = 3 000 cts.(0.4 µs) -1 . The right column shows the results for the decay parameters (d) λ 1 = 0.5 µs -1 , (e) λ 2 = 0.03 µs -1 , and (f) λ 3 = 0.0125 µs -1 .

  Figure 1

Table 1

 1

	Before RC	After RC

Table 4

 4

	Table 3

ACKNOWLEDGEMENTS

We are grateful to the three anonymous referees for making several useful recommendations which helped to improve the article.

TABLE CAPTIONS

Table 1 Amplitude and decay constant parameters obtained from fitting an artificial data set with an ISBR of 600 with 2, 3, and 4 components (columns 1-3) using the nonlinear least squares algorithm in SigmaPlot. The uncertainties in the parameter values quoted in these columns are internal error estimates calculated by SigmaPlot. In column 4, the average parameters are calculated from fittings of five different realisations of TR-OSL curves for the identical parameters, noise distribution, and an ISBR of 600. The quoted uncertainties are one σ (standard deviation).