The luminescence phenomena of insulators and semiconductors (e.g., natural minerals such as quartz) have various application domains. For instance, Earth Sciences and archaeology exploit luminescence as a dating method. Herein, we present the R package RLumCarlo implementing sets of luminescence models to be simulated with Monte Carlo (MC) methods. MC methods make a powerful ally to all kinds of simulation attempts involving stochastic processes. Luminescence production is such a stochastic process in the form of charge (electron-hole pairs) interaction within insulators and semiconductors. To simulate luminescence-signal curves, we distribute single and independent MC processes to virtual MC clusters. RLumCarlo comes with a modularized design and consistent user interface: (1) C++ functions represent the modeling core and implement models for specific stimulations modes. (2) R functions give access to combinations of models and stimulation modes, start the simulation and render terminal and graphical feedback. The combination of MC clusters supports the simulation of complex luminescence phenomena.
Light is perhaps the most basic everyday experience. Light emission that is not caused by the heating of a substance is called luminescence or ‘cold light’. Various fields exploit this phenomenon. For instance, Earth Sciences and archaeology determine the timing of past events (e.g., last sunlight exposure or heating) with a technique called luminescence dating. Since 2012, the luminescence-dating (or more general trapped-charge dating) community has gradually adapted R as a universal tool to analyze, model, and visualize their data. Relevant related CRAN packages are: BayLum (Christophe et al. 2018; Bayesian modeling: Philippe et al. 2019), Luminescence (luminescence- data analysis, Kreutzer et al. 2012, 2020), numOSL (luminescence-data analysis, Peng et al. 2013; Peng and Li 2018), RLumModel (luminescence-data modeling, Friedrich et al. 2016, 2020), RLumShiny (graphical interface to functions for plotting and calculation in the framework of luminescence-data analysis, Burow et al. 2016, 2019), and tgcd (curve deconvolution, Peng et al. 2016; Peng 2020).
The luminescence production process is a stochastic process involving
discrete random state transitions of subatomic particles. In the case of
luminescence, this translates to electrons (and holes) moving to
different energy levels, e.g., in the crystal lattice of the natural
mineral quartz. Such processes are ideal for Monte Carlo (MC)
simulations, and their application has a long and propelling history in
physics (cf. Landau and Binder 2015). Figure
Our contribution, and so RLumCarlo, sits on precedent work by Pagonis and Chen (2015), Pagonis et al. (2014), Pagonis et al. (2019), and Pagonis et al. (2020). The included collection of energy-band models for different stimulation modes adapted to MC methods are valuable for, e.g., studying the impact of model parameters on the signal-related stochastic uncertainties or statistic effects in tiny, dosimetric systems. Technically, our approach is closely related to the simulation of birth-and-death processes (for a review on birth-and-death process cf. Novozhilov et al. 2006). Each simulation run describes a Markov process. However, in our case, we allow only a reduction of an initial number of particles (i.e., only death processes).
Herein, we will not derive the full theoretical background of the models, but we will focus on the technical aspects of the package design and the integration of the MC methods. Such a presentation was beyond the scope of previous articles (e.g., Pagonis et al. 2019, 2020), but it is likely of interest to a broader community.
We structured our contribution as follows. The introduction continues with a brief paragraph on luminescence and the term ‘cold light’. After that, we detail the rationales for our contribution by recalling conventional modeling approaches in the field. Readers familiar with these topics may safely skip this part. The subsequent section outlines the concept and the implementation of RLumCarlo, including code examples. The remainder addresses (A) the implementation of a virtual dosimetric system to simulate weak spatial correlation of dosimetric cluster groups. (B) We outline how RLumCarlo can simulate more complex models compared to other solutions, with respect to its strengths and limitations. An outlook outlining the potential to implement more interactions between models will close our manuscript.
Light emissions of semiconductors or insulators after exposure to
ionizing radiation is a luminescence phenomenon now and then paraphrased
as ‘cold light’. The term luminescence relates to light production not
purely caused by the heating of a substance, a condition called
‘incandescence’ or black body radiation, but a phenomenon expressing the
inherent capacity of a material (dosimeter) to emit light (energy) in
the ultraviolet to infrared wavelength range (e.g., Newton Harvey 1957; Mahesh et al. 1989). Heat-related luminescence phenomena of
solids have been explored systematically in physics since the 1930s
(Urbach 1930) to characterize materials and understand charge
transfers in dosimeters (e.g., McKeever 1983; Mahesh et al. 1989). The
amount of luminescence, in the context of this manuscript, correlates to
the energy absorbed by a dosimeter during ionizing irradiation. The
closest analogy to a dosimeter is a battery that can be charged by,
e.g.,
Luminescence phenomena have versatile use in the fields of personal, medical, and accidental dosimetry (e.g., Yukihara and McKeever 2011). As aforementioned, in Earth Sciences and archaeology, the luminescence of natural minerals gained considerable attractiveness as a dating method (luminescence dating). First attempts exploiting luminescence signals as a chronological tool reach back to the 1950s (Daniels et al. 1953; Houtermans and Stauffer 1957; Grögler et al. 1958). Nevertheless, it needed a few decades more before the method took off and became today one of the most frequently used dating methods on sediments for the last 250,000 years and beyond (e.g., Aitken 1985, 1998; Bateman 2019).
To explain luminescence production, Johnson (1939) and Randall and Wilkins (1945) introduced the first basic energy-band models. Today, most of the commonly accepted luminescence models use series of more or less complicated systems of differential equations (for an overview, see Chen and McKeever 1997; Bøtter-Jensen et al. 2003; Chen et al. 2011) employing energy-band models. Those models provide a proper phenomenological match with measured data for various experimental designs by simulating electronic transitions. ‘Conventional’ energy-band models available to simulate luminescence production are developed as a set of nonlinear differential equations. This brings some limitations:
The models become complex easily and cannot be solved analytically.
If numerical methods are used, some equations are numerically unstable, which may lead to wrong simulation results.
A convenient assumption in many of such models is a great abundance of spatially uniformly distributed traps and recombination centers. However, this is not always the most prudent assumption. A spatial correlation and cluster formation of centers may exist for various reasons (cf. Mandowski and Świaltek 1992; Chen et al. 2011; Horowitz et al. 2017).
Deterministic models do not consider stochastic uncertainties and simulated curves are ‘noise free’. These limits subsequent analyses for materials where such uncertainties would matter due to the low, finite, number of charge carriers, and in these scenarios, simulation results are used as reference data to test statistical models used for luminescence data analysis in general.
Modeling code for simulating luminescence production was often written with the tools at hand, e.g., Mathematica (e.g., Pagonis et al. 2006), which has led to a fragmentation of incompatible solutions. In 2016, Friedrich et al. (2016) introduced RLumModel (Friedrich et al. 2020), pooling available kinetic (non-MC) models available for the luminescence production in quartz. A tantamount suite of R code was presented simultaneously by Peng and Pagonis (2016). We will compare results from RLumModel and RLumCarlo at the end of this manuscript.
MC simulations offer an alternative and are indispensable if the simulation of defect clusters in combination with the analysis of stochastic uncertainties is desired. Usually, the underlying models are very simple, but can be combined to describe complex systems. Important early work simulating TL using MC methods goes back to Mandowski and Świaltek (1992) and Kulkarni (1994). Mandowski and Świaltek (1992) tried to overcome the prerequisite of a large number of sample carriers, and Kulkarni (1994) investigated MC methods to overcome very long calculation times encountered for numerical calculations in particular scenarios. Kulkarni (1994) (p. 103) also reported a “statistical fluctuation” (noise like scatter) caused by the MC simulations but considered this more as a disadvantage. Later, Pagonis et al. (2020) explicitly exploited this as a feature, similar to birth-and-death processes and their related random uncertainties, to investigate specifically the stochastic uncertainties.
Before we start to detail RLumCarlo, a preceding note of caution: Any attempt to answer the question of whether a particular model may better explain the one or the other effect measured in luminescence studies would open Pandora’s box (e.g., Horowitz et al. 2017). Consequently, we will not engage in such a discussion. What we have implemented so far in RLumCarlo can be modified and exchanged. However, the underlying design concept remains applicable.
RLumCarlo implements energy-band models in a modular approach. Each model can simulate only an isolated effect (e.g., a single curve, see below), but the package design allows various combinations, e.g., in the form of clusters. Hence, RLumCarlo can evolve beyond a specific mathematical model through a combination of simple models.
To that end, RLumCarlo differs fundamentally from RLumModel, where the collected models allow the simulation of complex phenomena and even entire measurement sequences (Friedrich et al. 2016) but are self-contained by design. In other words, simulations cannot evolve beyond a specific mathematical model selected by the user. In RLumCarlo, the implemented energy-band models can simulate only isolated effects (e.g., a single curve, see below), but the package design allows a combination in the form of clusters. Throughout the text, we will use the word ‘clusters’ to (1) ascribe virtual units used in the MC simulation to run independent random processes (henceforth MC clusters) and (2) to define groups of defects (defect clusters) e.g., defined by their spatial distance. Only the latter carry a physical meaning.
To date, RLumCarlo
ships three simple energy-band models (Figure
While the parameters differ from model to model and depend on the
stimulation mode (heat or light, continuous or ramped), key entities
remain alike across the models, such as the trap depth (the energy
difference of the electron state from the conduction band)
The conditions of the simulations are defined through these parameters,
with
Each model supports up to four different stimulation modes (Figure
As an example, we will detail the mathematical background and its implementation for delocalized transitions below. For all other models, we may refer to the cited literature as well as the package manual.
The basic implementation of the MC processes as a software algorithm
consists of two nested loops. The outer loop iterates over a time
n <- 1000
t <- 1:100
P <- 0.2
remaining <- numeric(length(t))
for (t in t) {
for (j in 1:n) {
if (runif(1) < P && n > 0)
n <- n - 1
}
if (n > 0)
remaining[t] <- n
}
For example, the algorithm starts with 1,000 particles. In time instant
The implementation for luminescence production in RLumCarlo is very similar. To exemplify the adaptation of the models to be run as an MC simulation, we have selected the one-trap one-recombination center (OTOR) model (based on Halperin and Braner 1960) for TL. Our description below follows Pagonis and Chen (2015).
The OTOR model for TL can be expressed with the following set of differential equations:
Beyond already mentioned symbols, we used in the equations
By assuming quasi-static equilibrium conditions (Chen et al. 2011)
the resulting TL intensity becomes the general one trap equation, GOT:
with
and the differential equation becomes a difference equation:
The factor run_MC_TL_<model>()
) and isothermal TL
(functions named run_MC_ISO_<model>()
) applying the localized or
delocalized model
and for TL, from tunneling transitions it reads:
with
In Figurerun_MC
. (2) These functions interface one C++ function, each via
Rcpp (Eddelbuettel et al. 2020) (for an
overview, see the vignette of
RLumCarlo). The R
functions provide a convenient user interface, and the C++ functions
constitute the workhorse, as shown in the modeling core (Figure
As indicated above in the example implementation algorithm, each
simulation run (Kulkarni 1994 used the term ‘particle tracking’)
starts with
Simulations start with a call of the respective function, e.g., for TL
using the DELOC
model run_MC_TL_DELOC()
or run_MC_TL_LOC()
for the
LOC
model, respectively (see Figures
results <- run_MC_TL_DELOC(
s = 3.5e12,
E = 1.45,
clusters = 10,
n_filled = 200,
R = 1,
times = 100:450)
The function parameter names follow the terminology used for the
mathematical expression of the models as closely as possible. For
example, E
and s
are R
(n_filled
is times
. The parameter clusters
sets the number
of MC clusters, i.e., the number of Markov chains. High numbers in
clusters
increase the confidence in the simulation output at the cost
of more computation time.
The output can be passed to a dedicated plot function
(plot_RLumCarlo()
). The function supports a couple of standard plot
arguments, such as main
for the title of the plot, which is passed
down to graphics::plot.default()
via ...
(type ?dots
in the R
terminal).
plot_RLumCarlo(
object = results,
legend = TRUE,
main = "(A) Delocalized transition")
The parameter n_filled
can be a vector enabling different starting
conditions for each MC cluster. FigureTL
stimulation using localized (FigureRLumCarlo_Model_Output
, which is a list
comprising a
multi-dimensional array
(one slice per MC cluster) with the resulting
luminescence signal and a numeric
vector for the stimulation time.
Currently we provide S3-generics for summary()
and c()
. The first
one is also used internally by plot_RLumCarlo()
to melt the array into
a data.frame
before plotting. The plot output adapts to the used
stimulation mode provided via an attribute with each output object.
A straightforward application for this kind of simulation is the study of the impact of physical parameters on the luminescence signal output and the estimation of the stochastic uncertainties, which cannot be achieved with the deterministic approach of differential equations.
We provide more, always up-to-date examples with the package vignette, where we also compiled a table with meaningful physical parameter ranges for each model.
The examples so far presented may not appear very sophisticated, and
still, they allow insight that goes beyond a simple educational purpose
of simulating luminescence based on phenomenological models.
Pagonis et al. (2020), who used a preliminary version of
RLumCarlo, addressed
in detail the stochastic uncertainties of TL and OSL models. These
uncertainties come into play in nano-dosimetric materials with a small
number of defect clusters where the “finite-size” (Mandowski and Swiatek 1991) of
the system starts to matter in terms of a presumed spatial correlation
of defect cluster groups. To some extent, this should also be true for
systems exposed to high-energy radiation causing defect clusters (e.g., Mandowski and Swiatek 1991; Mandowski and Świaltek 1992). Previously in this paper, we have
used the term ‘MC clusters’. For a start, in
RLumCarlo, ‘MC
clusters’ entail independent and continuous Monte Carlo Markov chains
employed to simulate luminescence production, starting with a particular
number of electrons in the system. Whether the processes are run in
parallel or sequentially has no impact on the outcome, except for
computation speed. In other words, ‘MC clusters’ carry no meaning
regarding the underlying physics. However, as mentioned above, ‘MC
clusters’ from different models (with the same stimulation mode) can be
concatenated (see Figure
To simulate a three-dimensional (dosimetric) system, we can add meaning to MC clusters by reinterpreting them as dosimetric clusters. From the modeling perspective, nothing changes, but MC clusters gain a connotation of having a physical meaning.
Figure
Such a system can be created in
RLumCarlo via
create_ClusterSystem()
. The function distributes points randomly with
their coordinates:
Then, the Euclidean distance between the points is determined with
stats:dist()
, which is used by stats:hclust()
to group the defect
clusters (stats:cutree()
, with the outcome shown in Figurestats:hclust()
and stats:cutree()
for defining the clusters is somewhat arbitrary and might be refined in
the future. Therefore, more research, supported by measurements, is
needed.
Now, any function from
RLumCarlo can be used,
and the output of create_ClusterSystem()
is taken as input for the
argument clusters
. For example:
run_MC_TL_LOC(s = 3.5e12, E = 1.45, n_filled = 1000,
clusters = create_ClusterSystem(100))
creates a system with 100 randomly distributed defect clusters. If the
simulation is run in such a mode, the meaning of n_filled
changes.
Previously, it defined the number of electrons in each cluster
(
In the remainder, we want to compare simulation results from RLumCarlo, with other types of solutions, such as RLumModel which uses coupled differential equations to simulate luminescence production. RLumModel was selected since it was developed by some of the authors of this contribution. However, in theory, simple scripts using any existing models to simulate luminescence should work as well (as long as the models are comparable).
In contrast to RLumCarlo, RLumModel input values for physical parameters are preset. RLumModel encourages users to write a virtual luminescence signal measurement sequence, which is processed based on a pre-defined model with preset physical parameters.
For the comparison, we have selected a TL curve simulated with the luminescence model for quartz by Bailey (2001).
output <- RLumModel::model_LuminescenceSignals(
sequence = list(IRR = c(20, 10, 1), TL = c(20, 400, 1)),
model = "Bailey2001"
)
The results are shown in Figure
TL110 <- RLumCarlo::run_MC_TL_DELOC(
s = 5e+12, E = 0.97, R = 5e-10, times = seq(20, 400, 2),
N_e = output$`conc. level 1 (TL)`[1,2] / 1e+5)
N_e
was divided by a constant to reduce the computation time. The
dimensionless parameter R
corresponds to TL230
and TL325
) before all three objects were combined via:
object <- c(TL110, TL230, TL325)
and plotted on the top of the curve derived from RLumModel:
RLumCarlo::plot_RLumCarlo(
object = object,
plot_value = "sum",
add = TRUE,
FUN = function(x) {
x * 1/(1 + (1e+7 * exp(-0.61/(8.617e-5 * (object$time + 273)))))}
)
The argument plot_value = "sum"
was used to plot the total count sum
instead of its average. The additional function injected via the
argument FUN
corrects the TL curves for a phenomenon known as thermal
quenching (Wintle 1975). This is a reduction of luminescence
production efficiency at higher temperatures. The chosen quenching
parameters follow roughly data measured for quartz by Friedrich et al. (2018).
In summary, the results in Figure
The modeling of luminescence phenomena (cold light) of semiconductors and insulators after having received ionizing radiation is a challenging task. MC methods allow setting up flexible and simple systems to simulate luminescence with a finite number of charge carriers. This enables users to address effects usually observed for nano-dosimetric systems, and it provides insight into the stochastic uncertainty structure. We presented RLumCarlo, which renders, to our best knowledge, the first open-source and ready-to-use compilation of basic MC luminescence models for different stimulation modes (so far CW-OSL, LM-OSL, ISO-TL, and TL). We showed that the output from different models, which are simulated in separate MC chains in virtual clusters, can be combined to either simulate more complex systems or to mimic simple spatial correlations between cluster groups. The way of the implementation does not limit RLumCarlo to a specific dosimeter (e.g., quartz). In this light, RLumCarlo can be used in education, and research to test the impact of model parameters, such as cluster sizes and related stochastic uncertainties. Furthermore, RLumCarlo can help in in formulating research hypotheses and test them with commonly accepted or new models, still to be developed.
Future work will implement more models to run as MC simulation, e.g., for irradiation processes in crystals (including its luminescence output: radiofluorescence) and for an advanced interaction of clusters.
We thank two anonymous reviewers for their thorough reviews and
constructive suggestions, which helped to improve our manuscript.
Furthermore, we are grateful to the CRAN team in general for their
tireless efforts in keeping such a great resource alive and in
particular for their patience during the initial submission of
RLumCarlo. Alex Roy
Duncan is thanked for his support on the package development during our
stay in Westminster. The development of
RLumCarlo benefited
from the support of various funding bodies. The initial work by JF, SK
and CS was supported by the Deutsche Forschungsgemeinschaft (2015–2018,
DFG SCHM 3051/4-1, “Modelling quartz luminescence signal dynamics
relevant for dating and dosimetry”). Later financial support was
secured through the project “ULTIMO: Unifying Luminescence Models of
quartz and feldspar” granted by the Deutscher Akademischer
Austauschdienst (DAAD PPP USA 2018, ID: 57387041). SK was supported by
the LabEx LaScArBx (ANR - n.
RLumCarlo, BayLum, Luminescence, numOSL, RLumModel, RLumShiny, tgcd, scales, ggplot2, Rcpp, parallel, doParallel, foreach
HighPerformanceComputing, NumericalMathematics, Phylogenetics, Spatial, TeachingStatistics
This article is converted from a Legacy LaTeX article using the texor package. The pdf version is the official version. To report a problem with the html, refer to CONTRIBUTE on the R Journal homepage.
Text and figures are licensed under Creative Commons Attribution CC BY 4.0. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".
For attribution, please cite this work as
Kreutzer, et al., "RLumCarlo: Simulating Cold Light using Monte Carlo Methods", The R Journal, 2021
BibTeX citation
@article{RJ-2021-043, author = {Kreutzer, Sebastian and Friedrich, Johannes and Pagonis, Vasilis and Laag, Christian and Rajovic, Ena and Schmidt, Christoph}, title = {RLumCarlo: Simulating Cold Light using Monte Carlo Methods}, journal = {The R Journal}, year = {2021}, note = {https://rjournal.github.io/}, volume = {13}, issue = {1}, issn = {2073-4859}, pages = {351-365} }