Emissivity: A Program for Atomic Emissivity Calculations

Emissivity: A Program for Atomic Emissivity Calculations Taha Sochi∗ April 16, 2012 University College London - Department of Physics & Astronomy - Gower Street - London. Email: t.sochi@ucl.ac.uk. ∗ 1 Contents Contents 2 Nomenclature 3 1 Abstract 3 2 Program Summary 4 3 Theoretical Background 5 4 Long Write-up 7 5 Input and Output 5.1 Input . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Output . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 12 14 6 Acknowledgement 16 7 References 17 8 Nomenclature and Notation 18 2 1 ABSTRACT 1 3 Abstract In this article we report the release of a new program for calculating the emissivity of atomic transitions. The program, which can be obtained with its documentation from our website www.scienceware.net, passed various rigorous tests and was used by the author to generate theoretical data and analyze observational data. It is particularly useful for investigating atomic transition lines in astronomical context as the program is capable of generating a huge amount of theoretical data and comparing it to observational list of lines. A number of atomic transition algorithms and analytical techniques are implemented within the program and can be very useful in various situations. The program can be described as fast and efficient. Moreover, it requires modest computational resources. 2 PROGRAM SUMMARY 2 Program Summary Title of the program: Emissivity Type of program: command line Programming language: C++ Number of lines: 3110 Computer: PC running Linux or Windows Installation: desktop Memory required to execute: case dependent Speed of execution: case dependent (∼ few seconds) Has code been vectorized?: Yes Has code been parallelized?: No Compilers tested: g++, Dev-C++ Number of warnings: 0 Programming methodology: procedural with object oriented 4 5 3 THEORETICAL BACKGROUND 3 Theoretical Background In a thermodynamic equilibrium situation an excited atomic state is populated by recombination and radiative cascade from higher states, and depopulated by autoionization and radiative decay to lower states. Many recombination lines arise from radiative decay and subsequent cascade of strongly autoionizing resonance states near the ionization limit. These lines are dominated by low temperature dielectronic recombination. The populations of such resonance states are determined by the balance between autoionization and radiative decay. When autoionization dominates, the populations are then given by the Saha equation for thermodynamic equilibrium NX(n−1)+ = Ne NXn+ gX(n−1)+  2gXn+ h2 2πme kTe 3/2 −∆Et /kTe e (1) where ∆Et is the energy of the recombined electron in the X(n−1)+ state relative to the ionization threshold, and the other symbols have their usual meaning as given in Nomenclature § 8. The Saha equation, which describes the ratio of different stages of ionization, is based on the assumption of Local Thermodynamic Equilibrium (LTE) in a gas where collision dominates other physical processes. Consequently, the local velocity and energy distributions of particles are given by the Maxwell and Boltzmann distributions respectively and a temperature can be defined locally. The Saha equation is therefore strictly applicable only if elastic collisions are responsible for establishing the energetic distribution of particles. In many practical cases, however, atomic processes such as radiative transitions or dynamic effects are more important than elastic collisions and the assumption of LTE is not justified within the whole energy range. In these cases explicit detailed equilibrium calculations are required to determine the velocity and energy distributions of particles over the 3 THEORETICAL BACKGROUND 6 various energy levels [1]. To measure the departure of the state from thermodynamic equilibrium a departure coefficient, bu , is defined as the ratio of autoionization probability to the sum of radiative and autoionization probabilities. The value of this coefficient is between 0 for radiative domination and 1 for thermodynamic equilibrium. A large value of bu (≃ 1) is then required to justify the assumption of a thermodynamic equilibrium and apply the relevant physics. In thermodynamic equilibrium (TE) the rate of radiationless capture equals the rate of autoionization, giving [5] Ne Ni ̺c = NuS Γau (2) where Ne and Ni are the number density of electrons and ions respectively, ̺c is the recombination coefficient for the capture process, NuS is the Saha population of the doubly-excited state and Γau is the autoionization transition probability of that state. In non-TE situation, the balance is given by Ne Ni ̺c = Nu (Γau + Γru ) (3) where Nu is the non-TE population of the doubly-excited state and Γru is the radiative transition probability of that state. The departure coefficient bu is a measure of the departure from TE, and hence is the ratio of the non-TE population to the Saha population. Comparing Equations (2) and (3) gives bu = Nu Γau = NuS Γau + Γru (4) 4 LONG WRITE-UP 4 7 Long Write-up ‘Emissivity’ code is a command line program that was developed to implement the atomic emissivity model. Its main functionality is to calculate the emissivity of the transition lines from electronic recombination and cascade decay all the way down to the ground or a metastable state, that is all free-free, bound-free and boundbound transitions. ‘Emissivity’ is written in C++ computer language and mixes procedural with object oriented programming. It consists of about 3000 lines of code. The program was compiled successfully with no errors or warnings using Dev-C++ compiler on Windows, and g++ compiler on Cygwin and several Linux distributions such as Red Hat Enterprise. Sample results produced by these three versions were compared and found to be identical. Elaborate internal checks were carried out in all stages of writing and debugging the program and the output was verified. Thorough independent checks on sample emissivity data produced by ‘Emissivity’ were performed and found to be consistent. The program reads from plain text input files and writes the results to a main plain text output file. Other secondary output files can also be produced for particular purposes when required. The necessary input files are • The main file to pass the parameters and inform the program of the required calculations. The main parameters in this file will be outlined in § 5.1. • A file for passing the free-states data and oscillator strengths (f -values) for free-free and bound-free transitions with photon energies. The free-states data include an index identifying the resonance, its energy position, width, configuration, term, 2J, parity, and a flag marking the energy position data as experimental or theoretical. • A file containing the bound-states data. These data include an index iden- 8 4 LONG WRITE-UP tifying the state, its energy, effective quantum number, configuration, term, and a flag marking the energy data as experimental or theoretical. • A file containing the oscillator strengths (f -values) for the bound-bound transitions. The last file is imported from the R-matrix code [2, 4] output while the rest are user made. Two other input data files are also required if comparison and analysis of observational data are needed. One of these is a data file that contains astronomical observations while the other includes mapping information of observational lines to their theoretical counterparts. The astronomical data file contains data such as observed wavelength, observed flux before and after correcting for reddening and dust extinction, ionic identification, laboratory wavelength, multiplet number, lower and upper spectral terms of transition, and statistical weights of lower and upper levels. The mapping file contains the indices of the observational lines and their theoretical counterparts. In this section we summarize the theoretical background for emissivity calculations as implemented in the ‘Emissivity’ code • The program starts by obtaining the radiative transition probability Γrul for all free-free transitions as given by Γrul α3 Ep2 gl ful = 2gu τo (5) where the photon energy Ep is in Ryd, and the other symbols are defined in Nomenclature § 8. This is followed by obtaining the radiative transition probability Γrul for all free-bound transitions, as in the case of free-free transitions. 9 4 LONG WRITE-UP • The total radiative transition probability Γru for all resonances is then found. This is the probability of radiative decay from an upper resonance state to all accessible lower resonances and bound states. This probability is found by summing up the individual probabilities Γrul over all lower free and bound states l for which a transition is possible according to the electric dipole rules; that is Γru = X Γrul (6) l • The departure coefficient, bu , for all resonances is then obtained bu = Γau Γau + Γru (7) where Γau and Γru are the autoionization and radiative transition probabilities of state u, and Γa is given by Γa = ∆r ~ (8) • The next step is to calculate the population of resonances by summing up two components: the Saha capture, and the radiative decay from all upper levels. In thermodynamic equilibrium (TE) the rate of radiationless capture equals the rate of autoionization, giving Ne Ni ̺c = NlS Γal (9) where Ne and Ni are the number density of electrons and ions respectively, ̺c is the recombination coefficient for the capture process, NlS is the Saha population of the doubly-excited state and Γal is the autoionization transition probability of that state. In non-TE situation, the population and depopulation of the autoionizing state due to radiative decay from upper states and to 10 4 LONG WRITE-UP lower states respectively should be included, and hence the balance is given by Ne N i ̺ c + X Nu Γrul = Nl (Γal + Γrl ) (10) u where Nl is the non-TE population of the doubly-excited state, Γrl is the radiative transition probability of that state, Γrul is the radiative transition probability from an upper state u to the autoionzing state l, and the sum is over all upper states that can decay to the autoionzing state. Combining (9) and (10) yields NlS Γal + X Nu Γrul = Nl (Γal + Γrl ) (11) u On manipulating (11) the following relation can be obtained Nl  X Nu Γrul Γal + = Γal + Γrl Γrl + Γal u X Nu Γr ul S = Nl b l + r Γl + Γal u NlS  (12) where bl is the departure coefficient of the autoionizing state. This last relation is used in calculating the population. • The next step is to calculate Γrul for the bound-bound transitions. In these calculations the f -values can be in length form or velocity form, though the length values are usually more reliable, and hence ‘Emissivity’ reads these values from the f -values file produced by the R-matrix code. This is followed by finding Γru for the bound states by summing up Γrul over all lower bound states l, as given earlier by (6) for the case of resonances. • The population of the bound states is then obtained Nl = X Nu Γr ul u Γrl (13) 11 4 LONG WRITE-UP where u includes all upper free and bound states. • Finally, all possible free-free, free-bound and bound-bound transitions are found. The emissivity of all recombination lines that arise from a transition from an upper state u to a lower state l is then computed using the relation εul = Nu Γrul hν (14) where ν is the frequency of the transition line. • Apart from the normal debugging and testing of the program components to verify that they do what they are supposed to do, two physical tests are incorporated within the program to validate its functionality and confirm that no serious errors have occurred in processing and producing the data. These tests are the population-balance test and the metastable test. The first test relies on the fact that the population of each state should equal the depopulation. This balance is given by the relation X Nj Γrji = Ni j>i X Γrik (15) k<i The second test is based on the fact that the total number of the electrons leaving the resonances in radiative decay should equal the total number arriving at the metastable states. This balance is given by the relation X ∀j Nj Γrj = X Nk Γrki (16) ∀i,k>i where i is an index for metastable states and j is an index for resonances. 5 INPUT AND OUTPUT 5 12 Input and Output The detailed technical description of the program input and output files is given in the program documentation which can be obtained from www.scienceware.net website. However, in this section we give a general description of the nature of the input parameters and the expected output data so that the user can decide if the program is relevant for the required purpose. 5.1 Input Various physical and computational parameters are required for successful run. Some of these are optional and are only required for obtaining particular data. The main input parameters are • The temperature in K and the number densities of electrons and ions in m−3 . The program uses a temperature vector to generate data for various temperatures. These values can be regularly or irregularly spaced. • The residual charge of the ion, which is required for scaling some of the data obtained from the R-matrix code output files. • Boolean flag for reading and processing observational astronomical data obtained from an external input data file. A mapping file is also required for this process. If the user chooses to read astronomical data, the data lines will be inserted between the theoretical lines in the main output file according to their observed wavelength. • Boolean flag to run the aforementioned physical tests. • Parameters for generating and writing normalization emissivity data for theoretical and observational lines to the main output file alongside the original 13 5.1 Input emissivity data. The normalization can be internal using the emissivity of one of the transition lines produced by the program. It can also be with respect to an outside set of emissivity values corresponding to the input temperatures. Another possibility is to normalize with respect to an outside set of emissivity data in the form of effective recombination coefficients corresponding to the specified input temperatures. In this case the wavelength of the line to be normalized to is required. The effective recombination coefficient ̺f (λ) is defined such that the emissivity ε(λ) of a transition line with wavelength λ is ε(λ) = Ne Ni ̺f (λ) hc λ (J.m−3 .s−1 ) (17) Another choice for normalization is to be with respect to an outside set of emissivity data in the form of effective recombination coefficients corresponding to a set of temperatures that may be different from those used in the input. In this case, the recombination coefficients corresponding to the input temperature values are obtained by interpolation or extrapolation using a polynomial interpolation routine. • Parameters to control the production and writing of the effective recombination coefficients, ̺f , which are equivalent to the given emissivities. Vacuum wavelength of the transition lines will be used due to the restriction on the air wavelength formula (λ ≥ 2000Å). • Parameters to control the algorithm to minimize the sum of least square deviations of emissivity between the observational lines and their theoretical counterparts over the input temperature range. This sum is given by S= X ∀ lines N εN o − εt
2 (18) 14 5.2 Output N where εN o and εt are the normalized observational and theoretical emissivi- ties respectively. This algorithm also offers the possibility of computing the temperature confidence interval and its relevant parameters. The individual square differences for each one of the least squares and their percentage difference which is given by  N εN o − εt εN o  × 100 (19) can also be written to the main output file. The confidence interval for the goodness-of-fit index χ2 can also be found if required by providing the degrees of freedom η and the change in the goodness-of-fit index ∆χ2 . • Parameters to run and control an algorithm for finding the decay routes to a particular state, bound or free, from all upper states. As the number of decay routes can be very large (millions and even billions) especially for the low bound states, a parameter is used to control the maximum number of detected routes. The results (number of decay routes and the routes themselves grouped in complete and non-complete) are written to an output file other than the main one. The states of each decay route are identified by their configuration, term and 2J. For resonances, the Saha capture term and the radiative decay term are also given when relevant. These temperaturedependent data are given for a single temperature, that is the first value in the input temperature vector. 5.2 Output The program produces a main output file which contains the essential emissivity data. The file starts with a number of text lines summarizing the input data used and the output results, followed by a few commentary lines explaining the symbols 5.2 Output 15 and units. This is followed by a number of data lines matching the number of transitions. The data for each transition includes an index identifying the transition, status (FF, FB or BB transition), the experimental state of the energy data of the upper and lower levels, the attributes of the upper and lower levels (configuration, term, 2J and parity), wavelength in vacuum, wavelength in air for λ ≥ 2000Å, the emissivities corresponding to the given temperatures, the normalized emissivities and the effective recombination coefficients corresponding to the given emissivities. Writing the normalized emissivities and the effective recombination coefficients is optional and can be turned off. As mentioned earlier, if the user chooses to read the astronomical data, the observational lines will be inserted in between the theoretical lines according to their lab wavelength. There are other subsidiary output files. One of these is a file containing the results of the minimization algorithm. This file contains information on η, ∆χ2 , a, temperature at minimum χ2 , and the temperature limits for the confidence interval. This is followed by the temperature array with the corresponding least squares residuals, the aχ2 and the χ2 arrays as functions of temperature. Another output file is one containing mapping data used for testing purposes. A third output file is a decay routes file which was outlined previously. 6 ACKNOWLEDGEMENT 6 16 Acknowledgement The author would like to acknowledge the essential contribution of Prof. Peter Storey in all stages of writing, testing and debugging the program and providing the theoretical framework. Without his help and advice the development of the program would have been an impossible mission. 7 REFERENCES 7 17 References [1] Benoy D.A., Mullen J.A.M. and Schram D.C. (1993) Radiative energy loss in a nonequilibrium argon plasma. Journal of Physics D 26(9): 1408-1413. 6 [2] Berrington K.A., Eissner W.B. and Norrington P.H. (1995) RMATRX1: Belfast atomic R-matrix codes. Computer Physics Communications 92(2): 290-420. 8 [3] National Institute of Standards and Technology (NIST). URL: http://www.nist.gov. 20 [4] The UK Atomic Processes for Astrophysical Plasmas. URL: http://amdpp.phys.strath.ac.uk/UK_RmaX/codes.html. 8 [5] Seaton M.J. and Storey P.J. (1976) Di-electronic recombination. In: Atomic Processes and Applications. North-Holland Publishing Company. 6 8 NOMENCLATURE AND NOTATION 8 18 Nomenclature and Notation α Γa Γr Γrul ∆χ2 ∆r ∆E ǫo ε εul εN o εN t η λ ν ̺c ̺f τo χ2 fine structure constant (= e 2 /(~c 4πǫo ) ≃ 7.2973525376 × 10−3 ) autoionization transition probability (s−1 ) radiative transition probability (s−1 ) radiative transition probability from upper state u to lower state l (s−1 ) change in goodness-of-fit index width of resonance (J) energy difference (J) permittivity of vacuum (≃ 8.854187817 × 10−12 F.m−1 ) emissivity of transition line (J.s−1 .m−3 ) emissivity of transition line from state u to state l (J.s−1 .m−3 ) normalized observational emissivity normalized theoretical emissivity number of degrees of freedom wavelength (m) frequency (s−1 ) recombination coefficient for capture process (m3 .s−1 ) effective recombination coefficient (m3 .s−1 ) atomic time unit (= ~/Eh ≃ 2.418884326505 × 10−17 s) goodness-of-fit index a Å b bu c e e− E Eh Ep f ful scaling parameter in χ2 confidence interval procedure angstrom departure coefficient departure coefficient of upper state speed of light in vacuum (299792458 m.s−1 ) elementary charge (≃ 1.602176487 × 10−19 C) electron energy (J) Hartree energy (≃ 4.35974394 × 10−18 J) photon energy (J) oscillator strength oscillator strength of transition between upper state u and lower state l 8 NOMENCLATURE AND NOTATION g h ~ J k L m me N Ne Ni NS R r S t T Te Xn+ X(n−1)+ statistical weight in coupling schemes (= 2J + 1 for IC) Planck’s constant (≃ 6.6260693 × 10−34 J.s) reduced Planck’s constant (= h/2π ≃ 1.0545717 × 10−34 J.s) total angular momentum quantum number Boltzmann’s constant (≃ 1.3806505 × 10−23 J.K−1 ) total orbital angular momentum quantum number mass (kg) mass of electron (≃ 9.10938215 × 10−31 kg) number density (m−3 ) number density of electrons (m−3 ) number density of ions (m−3 ) Saha population (m−3 ) resonance matrix in R-matrix theory radius (m) total spin angular momentum quantum number time (s) temperature (K) electron temperature (K) ion of effective positive charge n bound state of recombined ion with effective positive charge n − 1 Abbreviations BB C F FB FF IC J K kg LTE (.)l m Bound-Bound Coulomb Faraday Free-Bound Free-Free Intermediate Coupling Joule Kelvin kilogram Local Thermodynamic Equilibrium lower meter 19 8 NOMENCLATURE AND NOTATION Ryd s TE (.)u (.)ul ∀ 20 Rydberg second Thermodynamic Equilibrium upper upper to lower for all Note: units, when relevant, are given in the SI system. Vectors are marked with boldface. Some symbols may rely on the context for unambiguous identification. The physical constants are obtained from the National Institute of Standards and Technology (NIST) website [3]. It should be remarked that most subscripts (e.g. l and u) are dummy variables, that is they are subject in their interpretation to the context and are not specific to the particular cases they are referring to.