mudirac: a Dirac equation solver for elemental analysis with muonic X-rays

We present a new open source software for the integration of the radial Dirac equation developed specifically with muonic atoms in mind. The software, called mudirac, is written in C++ and can be used to predict frequencies and probabilities of the transitions between levels of the muonic atom. In this way, it provides an invaluable tool in helping with the interpretation of muonic X-ray spectra for elemental analysis. We introduce the details of the algorithms used by the software, show the interpretation of a few example spectra, and discuss the more complex issues involved with predicting the intensities of the spectral lines. The software is available publicly at https://github.com/muon-spectroscopy-computational-project/mudirac.


Introduction
Muonic atoms, in which one electron is replaced by the much heavier muon, have long been known and studied as an interesting laboratory for many kinds of exotic effects and insights in physics [1]. Thanks to its higher mass, a muon surrounding an atomic nucleus has orbitals much smaller than an electron, and thus interacts far more strongly with the nucleus, probing effects such as vacuum polarisation [2] and the internal structure of the proton [3]. More recently, interest in using muonic atoms for elemental analysis has arisen [4,5], with experiments carried out at international muon facilities such as J-PARC and ISIS. Muons are implanted in a sample with a beam, and then the X-rays emitted during their cascade are monitored and used to determine its composition. In this way, muons provide a very powerful probe that's non-destructive and can be tuned in depth by varying the momentum of the incident beam. This technique is especially useful when analysing archaeological artifacts, which could be damaged by any approach that required extracting a sample from them. While there exist tabulated databases of muonic atom X-ray frequencies [6] and dedicated codes have been developed in the past to compute them [7], we still regarded it as useful to develop a new, modern code that could be used to aid with muon elemental analysis at ISIS. With this in mind, we have developed mudirac, an open source C++ solver for the radial Dirac equation, designed specifically for computation of muonic X-ray spectra. mudirac can compute muonic atom eigenenergies up to precision of less than 1 keV, accounting for the effects of finite nuclear size, vacuum polarizability, and electronic shielding, predicting both energies and transition probabilities for X-ray lines. This paper will provide an overview of the algorithms employed, a short documentation on how to use the software, and finally some tests and comparisons with reference sample spectra.

The radial Dirac equation
The Dirac equation for a radial potential can be reduced to a pair of first order one dimensional differential equations [8,9,10]: in which P and Q represent two components of the equation, with P being the 'large' one when close to the non-relativistic regime. Here m is the mass of the orbiting particle, or can be replaced by the effective mass, µ = (1/m + 1/M) −1 , including the effect of recoil from the nuclear mass M if one does not treat the nucleus as having infinite mass. In this equation and all the following atomic units are used, so that e = = 1 and c = 1/α ∼ 137. k is a quantum number that can take on any positive or negative integer values (but not zero). In terms of the angular momentum and total angular momentum quantum numbers used in the Schrödinger equation, l and j, we have that k = −(l + 1) j = l + 1/2 l j = l − 1/2 (2) . The full three-dimensional wave function can then be reconstructed as: where we make use of the spin spherical harmonics: with r|l, µ − s = Y l,µ−s (r) , Y lm being the usual spherical harmonics, and the Clebsch-Gordan coefficients as seen in [8,11], and remembering eq. 2 for the relationship between k and l. Particular attention should also be paid to nodal theorems for this equation, as they are a useful guide when searching for solutions. The rule is that the large component of the wavefunction, P , has always n − l − 1 nodes, while the Q component has n − l − 1 nodes if k < 0, and n − l nodes if k > 0 [12,13].

Potential
The base potential used in eq. 1 is the Coulomb potential, with an assumption of the nucleus being a uniformly charged sphere if a finite radius is used: . The radiuses of most isotopes of physical interest are extracted from the work by Angeli and Marinova [14]. It should be noted that the values tabulated in the paper are values for the root mean square radius, which means they need to be multiplied by 5/3 to get the equivalent radius of a uniformly charged sphere. For any isotope for which the radius is not known, we fall back on the empirical formula [1]: , where the second version is in atomic units. Two possible corrections to the Coulomb potential have been considered relevant for this software tool, as they've been shown to be largely the biggest contributions to transition energies, especially in heavy muonic atoms [1]. The first one is the so-called Uehling potential, representing the first order contribution of vacuum polarisability to the interaction between the nucleus and the muon -in other words, the contribution of a Feynman diagram in which a photon exchange between muon and nucleus is mediated through an electron-positron pair production loop. This potential has a well known form and can be written as for a point charge [15,16], and for a generic charge distribution ρ(r) [17]. For the specific case of a uniform spherical charge density, the integral becomes . The integral in x can be carried out if we consider the two distinct cases of x < r and x ≥ r. For the first we have while for the second . It then follows that the full term can be expressed as [X < (r, u| min(r, R)) + X > (r, u|R) − X > (r, u| min(r, R))] (16) which is the expression used for the Uehling potential in the software. The integration over u is left to be carried out numerically. The second correction adopted is the contribution of an electronic background charge. In this case the potential is found by solving a differential equation , which comes by simple application of Gauss's law in spherical coordinates.

Transition probabilities
To compute the transition probabilities between two states in presence of an electromagnetic field we need to compute the dipole matrix element between them. For spontaneous emission between two states a and b, in atomic units, we have [18] with E ab transition energy, the wave vector has magnitude k = E ab /c, and α is the vector of Dirac matrices α x , α y , α z . If we consider the wavefunction written as seen in eq. 3, then it follows that each individual matrix element can be written as with the σ i being the Pauli matrices, keeping in mind that . It should be recalled that in eq. 19 the minus sign appears because the part in the bra is the adjoint, and not merely the complex conjugate, of the wavefunction. Now, keeping in mind eq. 4 and 5, and shortening the Clebsch-Gordan symbols seen in eq. 7 as c s kµ = c(l, 1 2 ; µ − s, s), we have and therefore, . If we also set for brevity , and we consider that if to a certain k corresponds an l, then to −k corresponds l − sgn(k), then eq. 19 reduces to a| α x j 0 (kr) |b =i c . From these it becomes then possible to compute the transition probability as seen in Eq. 18. The transition rules are easily seen: only transitions with ∆l = ±1 and ∆µ = 0, ±1 are allowed. Given any two states, the integrals J ab can be computed numerically. While individual transition rates for different values of the magnetic quantum number µ can be computed, these are not in practice resolved, since the Hamiltonian we use leaves the energies degenerate in µ. Therefore, when considering a transition between two states with total angular momentum j b and j a , we actually compute an average transition rate , summing over all possible starting states and averaging over all possible 2j a + 1 degenerate end states.

Grid and boundary values
All calculations in mudirac make use of a logarithmic radial grid, so that r(x) = r 0 e x . In this way the radial Dirac equation becomes . Grid parameters are chosen case-by-case based on the parameters of the atom and the quantum numbers and energy of the state that's being integrated. In particular, by default we have , namely, r 0 is chosen as the maximum between the expected radius of an ideal Schrödinger 1s orbital for a system with nuclear charge Z and reduced mass µ and the finite radius of the nucleus. The boundary conditions are found by solving eq. 1 analytically for the r → 0 and r → ∞ limits. For a value of the energy E, defining K = m 2 c 2 − E 2 /c 2 ), we find at infinity [10]: . At zero the situation is more complex, and differs depending whether we're considering a finite or point-size nucleus. For a point-size nucleus we have [19] and for k > 0 .

Energy error estimation for eigenvalue search
The radial Dirac equation is solved by a shooting method, integrating each step with a fourth order Runge-Kutta (RK4) algorithm [20]. When using a shooting method to integrate a boundary value problem, one applies the following iterative process 1. start with a trial value for the eigenvalue, E; 2. integrate the equation forward from zero and backwards from infinity, using that trial value, meeting halfway at the turning point (namely, the point r t such that V (r t ) = E); 3. compute an error based on how well the left and right solutions match at r t ; 4. choose a new value of E to minimize the error and repeat from step 2 until a certain tolerance is reached.
The important thing to decide is how to evaluate the matching described at step 3, and then how to optimise the energy. One might be tempted to simply look at how the values of the forward-and backward-integrated functions match at the turning point, e.g. δP = P f (r t ) − P b (r t ), but that would be a mistake. The functions at this point can not be normalised, and so they are always known up to a constant term. Since this constant will be likely different between left and right integration (as it depends only on the boundary conditions), any absolute comparison is meaningless. Instead, it makes sense to consider the error which only depends on the ratio of Q and P . When E is an eigenvalue of the Dirac equation, namely, when a solution has been found, we expect δy = 0; more in general we expect we can find this energy by applying a steepest descent search and using the derivative of the error in the energy to minimize it. Now, making use of eq. 1 we find that we can write a differential equation for y where the apices indicate differentiation in r and we defined for convenience . Differentiating in E then we obtain (remember that ∂g + /∂E = −∂g − /∂E = 1/c). This differential equation can be easily integrated once we have a tentative Q and P . The boundary conditions are found by differentiating the ones for Q and P and are . Integrating ζ forwards and backwards from the two extremes of the grid, we can then apply eq. 42 and use it to minimise the energy. The only problem is that y has singularities in the poles of P , which causes practical difficulties in the integration of ζ. This is easily fixed though if we consider an equation complementary to eq. 45, namely which can be integrated in all those parts of the domain in which y grows too big. Performing a Runge-Kutta integration alternating between the two as needed it's possible to estimate ζ f − ζ b and thus converge the energy to an eigenvalue. From eq. 41 and 42 one can find that that Figure 1 shows an example of the results produced by this algorithm. For a muonic hydrogen atom and a range of initial energy guesses E 0 it shows what the algorithm predicts as the converged energy, E 0 − δE. If the linear approximation was correct, which of course tends to be true only when E 0 is close to an eigenvalue of the equation, then the figure would only show plateaus, with sudden jumps when the energy guess causes a new node to appear in the wavefunction, thus entering in the convergence basin of a different eigenvalue. These ideal plateaus, based on converged calculated energies and the number of nodes detected in the wavefunction, are shown as a solid line in figure 1. The curve for E 0 − δE appears instead as an approximation to it, but one can see how in practice a series of updates E ′ 0 = E 0 − δE does generally converge quickly to the appropriate eigenvalue. However, convergence tends to get more difficult when considering higher energy levels, where the gaps between successive shells shrink.

Electronic background potential
While in general the higher mass of the muon means its orbitals are much smaller than the electronic ones, and thus the contribution of the repulsion from electrons is negligible, there can be some overlap between some of the higher muonic orbitals and some of the lower electronic ones, which introduces a small correction to the energy. A full treatment of this correction would require solving the many-body problem for both the muon and all the electrons orbiting the atom, and at the moment has been considered unnecessarily complex for the scope of this software. Instead, the approach described by Tauscher et al. [21] has been used here. The scheme fills iteratively the electronic shells of the atom using the following rules: 1. the shell with n = 1 is treated as seeing a nucleus of charge Z − 1, to account for the shielding caused by the muon; 2. each successive shell is treated as seeing a nucleus of charge Z − 1 − q, with q the total number of electrons placed in lower shells.
The orbitals used are the standard hydrogen-like solutions of the Dirac equation [22]. The total background electronic density is then obtained as the sum of the individual density contributions from the radial parts of each of the filled orbitals. The exact density is calculated by using a user-defined electronic configuration. In order to compute the contribution to the potential of the electronic background, one needs to integrate eq. 17. Let's assume the general case of a background charge ρ(r) defined between the two limits R in and R out . On a logarithmic grid the equation becomes . This equation can not be integrated with the Numerov method that is most often used for second order differential equations due to the presence of a first order term. Instead, a third order scheme was used making use of the following finite differentiation approximations: , setting the first three points of V with an approximation of constant charge for r < R in and shooting outwards for all the following ones. In general, the potential can be constructed piecewise at all values of r as: with V grid the potential integrated on the area of the grid in which ρ is explicitly defined, and the convention V bkg (r) = 0 when r → ∞.

Software details
mudirac has been written to comply with the C++ 11 ISO standard [23]. It makes use of the libraries Catch, for testing, and Aixlog, for logging, which are included with the distribution. It is compiled into an executable that can be ran as follows: $> mudirac input.in with the input file input.in (the name is arbitrary) containing all necessary parameters to determine the calculation. These include the atom type, the terms to include in the potential, and which spectral lines should be calculated. A list of possible keywords with their meaning and default values is presented in Appendix A. The software produces the following output files: • a .log and an .err file, containing regular logging of the program's operation and any errors/exceptions thrown during the run. The amount of information in the log can be changed by setting a verbosity parameter in the input file; • a .xr.out file containing a full report of the energy (in eV) and transition probabilities (in s −1 ) for each requested spectral line, and excluding any lines for which calculations did not succeed or which were found to be forbidden transitions; • optionally, a .spec.dat file containing a final simulated spectrum including all lines convoluted with Gaussian functions and with intensities proportional to their transition probability; • optionally, a series of .<name>.out state files containing potential, P and Q components of the wavefunction for each computed state. The state names use the standard IUPAC X-ray spectroscopy convention [24]; • optionally, a series of .<name1>.<name2>.tmat.out files containing the details of the transition matrix between two given levels, namely, the transition rates between each of their respective degenerate states of different m j .
Line  Atomic units are used inside the code and for most input parameters (either energy or masses). The only exception are parameters for the writing of .spec.dat files, in which those that have units of an energy must be set in electronvolts. mudirac is open source software; it has been released publicly under the MIT License, and is available at the URL https://github.com/muon-spectroscopy-computational-p as part of the tools produced by the Muon Spectroscopy Computational Project.

Hydrogen lines
While mudirac was designed with muons in mind, it is possible to set the mass of the particle to arbitrary values, so solving for electronic atoms is not a problem. This is particularly important to test the quality of the transition probabilities, as it is not easy to find reference values for muonic atoms. As a first test, here we present the computed relativistic transition energies and probabilities for the hydrogen atom, compared to the well known and tabulated experimental values in the NIST Atomic Spectra Database. Calculations with mudirac were run using a finite size nuclear model and the Uehling potential correction on. The lines are labelled The results are shown in Table 1   Energies computed with mudirac make use of a finite size nucleus, Uehling correction to the potential, and an electronic background charge using the configuration of the atom with Z − 1 (so, for example, the configuration for 12 Mg is that of neutral Na). Please note that the values from [6] are manually extracted from images of spectra and thus may be less accurate than the others.

Muonic atom lines
Here we present a few examples of comparisons of spectral lines for muonic atoms to values reported in the literature, in particular in [1], [29], and the spectra from the database presented in [6]. For the values found in [1], a reference computed value is available, besides the experimental one. These computed values are extremely accurate, including also field theory effects beyond the Uehling potential, thus it's not surprising that the mudirac computed values differ slightly. Values from the other sources are all experimental. The comparisons can be seen in Table 2.

Isotope effect
It has been shown [30] that muonic X-ray analysis is also a good way to discern between isotopes of the same element, thanks to the effect of finite nuclear size and mass. Here we show the prediction made by mudirac for 204 Pb, 206 Pb, 207 Pb and 208 Pb, using as reference experimental data from the literature [29]. using a finite nuclear size, Uehling correction, but no electron screening. As it can be seen, the trend is reproduced very well, though there is a systematic error of around -50 keV for the K-L transitions, and of around +3 keV for the L-M one. This is likely due to one of the effects that have not been accounted for in the Hamiltonian [1]; it's interesting to note that it is in the same order of magnitude as the Uehling contribution to the line (∼30 keV), which suggests it is not likely to be a higher order vacuum polarisation term. It is also clear that this effect dramatically affects the inner shells, as transitions between higher ones, even in Pb, are far more accurate (as seen in Table 2. These will be considered for inclusion as the software develops, to guarantee higher accuracy. Nevertheless, the isotope effect is clearly well captured as the line shift from one isotope to another is much more accurate than the absolute value of the energies.

Spectra comparisons
Finally, we consider a few example experimental spectra and compare them with mudirac computed results. The experimental spectra were acquired at ISIS with the RIKEN Port 4 instrument, and come in multiple resolutions, ranging from 2 to 0.5 keV between each point. The samples are  pure metal standards. All calculations include a Uehling correction and the standard electronic configuration for the given element. Simulated spectra have been generated simply by using a Gaussian broadening of 1 keV. Energies are in keV, and intensities are arbitrary units. Figure 5 shows a comparison for the spectrum of aluminium, up to 450 keV. The experimental spectrum has resolution of 0.5 keV. The computed one includes all lines involving states up to the O shell (n=5). One thing that can be noticed is how the intensity ratios between far apart group of lines are remarkably different between the experimental and computed spectrum; for example the L-M lines appear much higher than the K-L ones in the experiment, whereas that relationship is inverted in the computed spectrum. The major factor contributing to this is the fact that the sensitivity of the instrument used to acquire the spectrum isn't constant, but goes down as the energy increases, and the spectrum itself hasn't been corrected for this. If one focuses on individual groups of lines, such as the L-M, L-N, and L-O ones, one sees that experimental and computed spectra exhibit similar patterns, as the lines are concentrated in a small enough interval of energies that the sensitivity of the instrument remains nearly constant. Another thing to consider however is that the intensities in the computed spectrum are based only on the probability of a transition to take place. In reality, the intensity will be controlled also by the population of the upper state from which the decay takes place. In other words, the peak intensities of the computed spectrum would be accurate only if all states in the muonic atom were equally likely to be occupied. This is in practice not the case. In fact, the probability distribution in each of the states in the lower shells heavily depends on the cascade process and the starting angular momentum of the muon upon capture [31,32]. The cascade can be computed given the transition rates, but one still needs to make assumptions about the initial angular momentum distribution. This is part of future plans for this project, but is currently not implemented in the software. Figure 6 shows a comparison between experiment and computation for gold, up to 8000 keV. The gold spectrum was computed up to the P shell (n=6). The comparison is split in different panels to highlight the details of the spectra at different scales; because the detector loses efficiency at high energy, the intensity has been rescaled between panels to better show the lines. The range up to 1000 keV uses data with a resolution of 0.5 keV; the 1000-8000 keV range uses data from a detector with a 2 keV resolution, which was the only one able to reach such high energies. The software allows us to identify various lines, from the K-L group at high energies (>5000 keV) to the very low energy transitions in between states in the same shell (<150 keV). Interpretation of the spectrum, however, is in this case somewhat complicated by the fact that this vast range includes high energy phenomena that aren't only muon level transitions. Due to the high overlap between muon and nucleus, multiple nuclear reactions take place, resulting for example in the transmutation of gold to platinum by the reaction µ − + p + → n + ν, with the excess energy ending up carried away by either the neutrino or the neutron, and the resulting nucleus being in an excited state that decays with emission of gamma rays. In particular, using the transition frequencies from the IAEA Nuclear Data database [33], the figure shows strong evidence that the reaction 197 Au + µ − → 196 Pt + n must be taking place, as several lines connected to the decay of excited states of 196 Pt are identifiable (dashed lines in Figure  6). In addition to these, the 511 keV line corresponding to electron-positron annihilation can be seen. We can also see how the K1-L2 and K1-L3 lines, which have energies above 5 MeV, display single and double escape peaks at -511 and -1022 keV, a phenomenon related to the creation of electronpositron pairs and relative loss of energy in the measured photons. This is a well-known problem with germanium detectors in this energy range [34].   Finally, we consider the case of a two-element system; a 70-30 Au-Ag alloy. In this case our objective is to distinguish peaks belonging to each element, allowing for some basic elemental analysis. The Au computed spectrum was the same used in the previous section; the Ag one was computed with the same parameters, and again, up to the P shell (n=6). Figure 7 shows the results. Here we focus on the low energies, as this is where the spectra are more rich in features, and thus interpretation is usually more difficult.
Some peaks are highlighted with their respective assignments. One interesting observation is that, again, intensities are dominated by capture and cascade mechanisms that obviously are not trivial; for example, despite the relative abundance of gold in the alloy, the M5-N7/M4-N6 doublet from silver is much higher than the N-O group from gold, running against the naive prediction based on transition probabilities. One must also pay attention to potential sources of ambiguity: for example, around 216 keV, one peak is assigned as the O9-P11/O8-P10 doublet from gold, but at the same frequency we can also see a peak in the silver spectrum. The structure of the peak, as well as its intensity, here lead us to believe it originates from gold; however it is hardly a foregone conclusion. Finally, another possible source of confusion must be considered in this low energy range: the lack of information about states above a certain energy. Computations were stopped at a given shell, which means that all transitions involving states above that are missing from our interpretation. As states get closer and closer in energy, going up in shells, more transition energies open up, leading to a nigh continuum of possible transitions in energy. This is further complicated by the fact that converging calculations for such higher states is also more difficult; in this example, the P10 and P11 states for Ag did not converge and the peaks they contribute to are missing from the spectrum. Going above this problem becomes only more severe, and affects our ability to properly interpret spectra in the low energy range. Improving these convergence properties is one of the future objectives for the development of the software.

Conclusions
A software to solve the radial Dirac equation specifically designed for muonic atoms has been produced and released for public use under an open source license. The software, called mudirac, is able to compute energies and transition rates for all the main X-ray spectral lines for muonic atoms across the periodic table, accounting for the most important effects such as the finite size of atomic nuclei and the Uehling correction. This provides a useful tool to help with the interpretation of muonic X-ray spectra for elemental analysis. The software has also been designed to easily accommodate future expansions such as the addition of more correction terms for its potential in order to achieve higher accuracy. In addition, future work will focus on using the software's predictions for transition rates in order to help interpret line intensities based on a model of the muonic cascade, as well as improving its ability to compute states in higher shells.

Acknowledgements
Thanks to Adrian Hillier, Dominik Jochym and Albert Bartok-Party for the useful discussions and contributions. The development of this software was funded by the Ada Lovelace Centre and was carried out in collaboration with the ISIS Muon Group.

Appendix A. Input keywords
As explained in section 4, mudirac takes a single input file, containing multiple lines with the format <keyword>: <value> . Here we list all the currently available keywords, divided by type, their purpose, and default values.
Appendix A.1. String keywords These keywords take a string as value; invalid strings (e.g. a chemical symbol that doesn't correspond to a known element) will give rise to errors.
• element: symbol of the element for the calculation. Determines the nuclear charge. Can be any symbol in the periodic table up to Z=111, Roentgenium (Rg). Default is H.
• nuclear model: model used to describe the nucleus. Can be POINT (point charge) or SPHERE (finite size, uniformly charged spherical nucleus). Default is POINT.
• electronic config: electronic configuration to use in order to describe the negative charge background. Can be a full string describing the configuration (e.g. 1s2 2s2 2p2), an element symbol to represent the default configuration of that atom when neutral (e.g. C) or a mix of the two (e.g. [He] 2s2 2p2). Default is the empty string (no electrons).
• ideal atom minshell: for this shell, and all above it, treat the atom as a simple hydrogen-like point charge Dirac atom, using the known analytical solution and discarding all corrections. Mostly useful for debugging, or when very high shell states have difficulty to converge. The shell must use IUPAC notation (K =⇒ n = 1, L =⇒ n = 2, etc.). Default is the empty string (no ideal solutions used).
• xr lines: the transition or transitions for which energy and rates are desired. Each line must be expressed using the conventional IUPAC notation [24]. Multiple lines can be separated by commas. For example xr_lines: K1-L2,K1-L3 . In addition, colons can be used to indicate ranges of lines. The notation K1:L3-M1 would compute the lines K1-M1, L1-M1, L2-M1 and L3-M1. Note that if some of these lines are forbidden by selection rules, they will simply be skipped. A double colon, like K1:L3-K1:L3 would loop on both sides, and not count all repeated lines.
• max dE ratio: maximum ratio between energy step, δE, and current energy E in convergence search. If the suggested step exceeds this ratio times the guessed energy, it will be rescaled. This also serves as a measure to avoid overshooting and can be tweaked to get around cases of bad convergence. Default is 0.1.
• node tol: tolerance parameter used to identify and count nodes in wavefunctions. Very unlikely to need any tweaking. Default is 1E-6.
• loggrid step: step of the logarithmic grid. Default is 0.005.
• logggrid center: center of the logarithmic grid in units of a 0 = 1/(Zm). Default is 1.
• uehling lowcut: low cutoff for Uehling potential, under which the radius will be considered 0. Default is 0.
• uehling highcut: high cutoff for Uehling potential, over which the radius will be considered r >> c. Default is INFINITY.
• econf rhoeps: charge density threshold under which the electronic charge background will be truncated and treated as zero. Default is 1E-4.
• econf rin max: upper limit for the innermost radius of the electronic charge background grid. Default is -1 (no limit).
• econf rout min: lower limit for the outermost radius of the electronic charge background grid. Default is -1 (no limit).
• spec step: energy step for the simulated spectrum, in eV. Only has effect if write spec = TRUE. Default is 1E2 eV.
• spec linewidth: Gaussian broadening width for the simulated spectrum, in eV. Only has effect if write spec = TRUE. Default is 1E3 eV.
• spec expdec: exponential decay parameter E dec for a sensitivity function for the simulated spectrum, in eV. Multiplies the entire spectrum by a function exp(−E/E dec ). Only has effect if write spec = TRUE.
1. will print out only the transition energies and rates in the .xr.out file; 2. will print out also each of the states in a separate ASCII file as well as the transition matrices for each line; 3. is reserved for future uses and currently has the same effect as 2.