ZEST: A Fast Code for Simulating Zeeman-Stark Line-Shape Functions

We present the ZEST code, dedicated to the calculation of line shapes broadened by Zeeman and Stark effects. As concerns the Stark effect, the model is based on the Standard Lineshape Theory in which ions are treated in the quasi-static approximation, whereas the effects of electrons are represented by weak collisions in the framework of a binary collision relaxation theory. A static magnetic field may be taken into account in the radiator Hamiltonian in the dipole approximation, which leads to additional Zeeman splitting patterns. Ion dynamics effects are implemented using the fast Frequency-Fluctuation Model. For fast calculations, the static ion microfield distribution in the plasma is evaluated using analytic fits of Monte-Carlo simulations, which depend only on the ion-ion coupling parameter and the electron-ion screening factor.


Introduction
Line shapes are key ingredients of opacity and emissivity codes, as they often serve as a diagnostics of laboratory or astrophysical plasmas.Indeed, profiles contain information on local electric fields produced by electron and ion perturbers (leading to Stark splitting), on the magnetized environment (yielding Zeeman patterns), on macroscopic variables of the plasma (density and temperature of particles), on collective phenomena, etc.When considering radiative transport in optically-thick plasmas, it well known that the Rosseland mean opacities may be sensitive to the line shapes, in particular due to the tails of these functions [1].More generally, there is a close connection between line profiles, radiative transfer and the atomic populations which describe the microscopic states of the plasma [2].Such a relationship is shown for instance by the escape factors formalism [3], which is often used in non-local-thermodynamical-equilibrium (non-LTE) collisional-radiative models to account for opacity effects on the population of states.Based on the fact that photons in the tail of lines can escape the plasma more easily than photons close to the center of lines, an effective reduction factor is introduced in the calculation of radiative transitions, which are then used in the rate equations to obtain the population of states.These escape factors depend closely on line-shape profiles I ν , with integrals of the type dνI ν e −τ 0 I ν /I 0 where τ 0 is the optical depth at the center of the line.It is often important to include Stark broadening in the calculation of these factors, though it is not an easy task [4].
The calculation of line shapes broadened by Zeeman-Stark effects is an old and complex subject [5,6] that is still debated today [7,8].The ZEST code (short name for ZEeman-STark), which is the object of this paper, was developed recently to improve the calculation of line shapes in our spectroscopic models in LTE [9] or non-LTE [10].The code was presented and used, for the first time, in the 4th Spectral Line Shapes in Plasmas code comparison Workshop (20)(21)(22)(23)(24) March 2017, Baden, Austria).The code is based on the Standard Line-Shape Theory where the ions are quasi-static during the time emission process, whereas electrons are treated in the framework of a binary collision relaxation theory.This simplification allows one to perform very fast calculations, while keeping a reasonable accuracy.A general description of the model can be found in Section 2, where we give details on the different assumptions made on the dipole auto-correlation function, the broadening operator, the implementation of Zeeman effect, the calculation of the static-ion microfield distribution and the treatment of ion dynamics effects using the fast Frequency-Fluctuation Model (FFM) [11].In Section 3, we show various numerical applications of the model.Conclusions and future developments will be discussed in Section 4.

Quasi-Static and Binary Collision Approximations for Stark Effect
In general, microfields generated by ion and electron perturbers evolve on different time scales because of their relative mass.Since the typical fluctuation rate of the electron field is greater than for ions by at least two orders of magnitude, one can often separate the contribution of both species: ions are considered as quasi-static during the emission process, whereas electrons are dynamic and treated as fast binary interactions with ions.The impact approximation for electrons assumes that the collisions are complete during the time of interest and is thus valid only near the center of the profile.This approximation may be softened by using a binary collision relaxation theory which can account for incomplete collisions, resulting in a frequency-dependent electron broadening operator (see Section 2.2) which preserves the static limit in the far wings.Moreover, the quasi-static approximation being valid only in the wings of the profile (large frequency detuning or small time interval), corrections due to ion-dynamics effects must be taken into account near the center of the line since this corresponds to radiation over time long enough for the ions to move (see Section 2.4).
The quasi-static line shape profile is given by: where P(F) is the normalized probability that the emitter undergoes the static-ion field F. The field-dependent profile reads: where L(F) = [H er (F), ] is the Liouville operator associated to the electron-radiator Hamiltonian H er , ρ is the normalized radiator density operator, d is the radiator dipole operator, φ(ω) is the electron broadening operator and C is a constant factor ensuring the normalization of the spectral profile.The trace in Equation ( 2) represents the average of the inverse Fourier Transform of the system's dipole-dipole autocorrelation function.It must be performed on the product of matrices ρ, d.d † and which involve tetradic operators, i.e., acting on four sets of internal states of the emitter corresponding to two initial and two final states.The inversion of matrix A in the complex domain is one the most time-consuming task in a line-shape code.In order to speed up the calculations, we make the following assumptions.First, perturbers are not allowed to induce a transition between the lower and upper states involved in the radiator emission process (no-quenching approximation).We only consider linear Stark effect, so that the electron-radiator Hamiltonian reads H er (F) = H 0 − d.F/h, H 0 being the unperturbed (no field) Hamiltonian.We assume that only states belonging to the same Layzer complex (see Section 3 for definition) may be mixed by Stark and/or Zeeman effects (∆n = 0 interaction channels).The broadening operator φ(ω) is calculated to the second-order in the radiator-electron interaction (dipole approximation).The imaginary part of φ(ω), which corresponds to small plasma-generated static level shifts, is neglected.We also neglect the interference terms between upper and lower states in φ(ω).It is important to mention that we plan to relax some of these approximations in future versions of the code (see discussions in Section 4).
Since the trace is invariant under any unitary transformations, it is more convenient to use the matrix P that diagonalizes the Liouville operator, L = P T LP, where L is a diagonal matrix in the dual space, with matrix elements noted i f |L |i f = ω i f (F) corresponding to the energy of the Stark component i → f .Numerical computations can use the fact that the matrix associated to L is real and symmetric, yielding a real transformation matrix P that can be easily computed and inverted (by taking its transpose P −1 = P T ).The problem is that the broadening operator φ is not diagonal with the eigenfunctions of L in general, i.e., when applying the transformation: φ = P T φP.For the conditions where homogeneous broadening is small compared to the inhomogeneous broadening produced by the static ions [12], it is, however, possible to neglect the off-diagonal matrix elements of φ and to retain only the diagonal matrix elements denoted by φ i f (this approximation will be relaxed in future versions of the code with the consideration of the interference terms).With all the above hypothesis, the trace can now be evaluated by the expression which involves easy inversion of a diagonal matrix.If we denote by {|i } the set of orthogonal wave-functions in the new basis, the line shape (2) can be simplified to the expression: The normalization factor, C = ∑ i f | i|d| f | 2 , is the total line strength of the transitions, which is invariant under any unitary transformation (it may be evaluated with the unperturbed wave-functions).
The calculation of matrix elements of Liouville, dipole and electron broadening operators in the new basis requires knowledge of matrix elements of these operators between unperturbed (no field) wave-functions, denoted by |i 0 .For an arbitrary operator O, we have: where P ii are the eigenvectors of the unitary transformation (each column of matrix P).In practice, we use an external atomic database generated with the multi-configuration Dirac-Fock code MCDF [13] or the cFAC code [14,15].The data stored are: the total angular momentum J, energy, parity and configuration label for a collection of levels related to the unperturbed Dirac Hamiltonian H 0 ; and the signed reduced matrix element for electric-dipole (E1) and magnetic-dipole (M1) [See Section 2.3] transitions between these levels.The use of Wigner-Eckart theorem allows one to compute the matrix elements i 0 |O| f 0 from the corresponding reduced-matrix elements in the database.For instance, the dipole matrix elements of the Stark Hamiltonian H S = − d.F/h between two unperturbed states read: where a 0 is the Bohr radius, C (1) is the spherical harmonic tensor operator of rank 1, and using a level description based on the total angular momentum J, its projection M and some other quantum numbers α to account for Pauli exclusion principle.For hydrogenic emitters, the required atomic data (energies, reduced matrix elements) may be evaluated directly using analytical formulas [16].It is worth mentioning that perturbers (ions or electrons) are assumed to produce a randomly oriented microscopic field, so that the emission spectrum is not polarized in general.The quantization axis may be thus chosen arbitrarily.
As concerns the population of states, ρ i = i|ρ|i , we have different options: 1. Local Thermodynamics Equilibrium (LTE): where k B is the Boltzmann constant and T e is the electron temperature.2. Statistical Weight Approximation (SWAP): where g is a constant factor ensuring that ∑ i ρ i = 1.This corresponds to the high-temperature limit of LTE plasmas.3. Non-Local Thermodynamics Equilibrium (non-LTE): where the populations of states j 0 |ρ|j 0 NLTE are derived from results of an external collisionalradiative model [10] which do not account for the splitting of levels by Stark effect (the CR model is based on degenerate levels).Denoting by p γ j J j the non-LTE probability of the parent level of state j, we use a SWAP assumption within each level: where J j is the total angular momentum of the level.
We have the possibility to include other homogeneous broadening mechanisms: Doppler effect, natural width, autoionization or any other radiative or collisional processes in the plasma limiting the lifetime of states.Lifetime of states is estimated from collisional-radiative calculations by using the diagonal terms of the rate matrix.If we assume that the matrix elements of the electron broadening operator do not depend on the photon frequency φ k (ω) φ k (it is evaluated at the line center, i.e., impact approximation: see Section 2.2), then the convolution of expression (2) with other Lorentzian or Gaussian functions leads to a sum of Voigt functions (see Appendix A).The integration over the ion microfield distribution in (1) may be readily estimated by a sum over a well-chosen set of field values.As a result, the quasi-static line shape can be written in compact form as a sum of Voigt profiles arising from different Stark components and different ion field values, where the normalized weight factors f k (such that ∑ k f k = 1) include the probability of the ion microfield and the strength of the relevant Stark component of energy hω k .The width parameters of the Voigt functions are defined by: where T i is the ion temperature, M i is the ion mass, ∆ω exp is the intrumental FWHM (possibly taking into account the spectral resolution of the spectrometer in an experiment), φ k is the electron impact width and Γ k is an additional homogeneous broadening width (natural, autoionization, etc.).

Electron Broadening Operator
The main contribution in φ(ω) comes from weak electron collisions of duration much shorter than the inverse HWHM of the profile, i.e., much shorter than the decay timescale of the dipole autocorrelation function.Using the perturbation theory up to the second-order in the radiator-electron interaction, we obtain for the Maxwell-averaged operator: where ∆ω is the detuning frequency from line center, m is the electron mass, N e is the electron density.
There are many ways to estimate the G-function.It may be related to charge-density fluctuations in the plasma by using the dielectric constant (∆ω, κ) of an equilibrium plasma as follows [17]: where κ is the wavenumber associated to the momentum transfer in the radiator-electron collision.
The dielectric function may be estimated in the Random-Phase Approximation for a Maxwellian plasma [18], where x 0 e t 2 dt is the Dawson function.A cutoff κ m is introduced in order to limit the momentum transfer due to close (and strong) collisions and to avoid the divergence of the integral.For hydrogenic ions of charge Z, radiating from the upper shell of principal quantum number n, we use: The numerical evaluation of the integral ( 9) is however rather painful.Following the model of Lee [19], it is possible to approximate this cumbersome integral by an analytic expression using the asymptotic limits of the G-function at ∆ω → 0 and ∆ω → ∞: where and t dt is the exponential integral of first order and ω p = e 2 N e m 0 is the plasma frequency.The semiclassical GBK model [20] yields a G-function which is slightly different from Lee's model, in particular near the center of the line: This expression is valid only for classical plasmas (low-density and high-temperature limit).In practice, we use a modified formula for the electron-ion screening length which has a wider range of applicability (see Section 2.5).
where Ryd is the Rydberg energy.Using the same cutoff expression (11), we may rewrite the above formula as: which is equivalent to Lee's formula G ∞ (∆ω) when ∆ω → ∞.
A comparison of the G-functions computed for κ m λ D = 30 with different methods and approximations is shown in Figure 1.The G-function is obtained: (i) numerically from Equations ( 9)-( 11) for Dufty's model; (ii) Equations ( 12) and ( 13) for Lee's model; (iii) by Equation ( 15) for GBK.In Dufty's model, we observe a saturation of the G-function at low detuning: this is due to the fact that the exchange of electromagnetic energy h∆ω between the electron subsystem and the radiator is inefficient below the plasma frequency for a classical plasma.We note that this behaviour is well reproduced by the approximate model of Lee, but not in GBK.The discrepancy with the GBK model amounts to roughly 15% near the center of the line, whereas all models have the same asymptotic behaviour for large detuning.We have plotted in Figure 2 the ratio of the electron impact widths at the line center, obtained with the model of Lee and the model of GBK, as a function of the parameter κ m λ D .The latter quantity represents the ratio of maximal over minimal impact parameters for the electron-radiator collision.As we can see, both models agree only when κ m λ D → ∞ (low-density and/or high-temperature limit).We observe a discrepancy of less than 20% for κ m λ D > 20.For smaller parameters, the differences can go up to a factor of 3; however, the validity of semiclassical approaches is questionable as κ m λ D tends to one.For classical plasmas (low-density and high-temperature limit), the parameter κ m λ D is often very large so that one can further simplify the expression of the G-function.In the limit ∆ω → 0, it is easy to show that both models (GBK and Lee) yield an electron collision width depending on the Coulomb logarithm 2 , as expected in the impact approximation: for the GBK model (γ e is the Euler's constant) and for the model of Lee.The noticible difference is the factor − 1 2 which explains why Lee's model yields smaller impact widths than GBK in general (see dashed curve in Figure 2).
The collisions whose impact parameter is lower than the cutoff ρ m = 1/κ m cannot be treated in the framework of the perturbation theory.To account for these so-called strong (or close) collisions, which are missing in the above models, a constant term G n is added to the G-function as suggested in the reference [20], where n is the principal quantum number of the active electron,  It is defined as the logarithm of the ratio of the upper and lower impact parameter cut-offs for a binary elastic collision.
Neglecting the interference terms, the diagonal matrix elements of the broadening operator read: This frequency-dependent expression is used only when the line shape is evaluated through Equations ( 1) and (3).When the Doppler effect is taken into account, we introduce the approximation G(∆ω) G(0) so that the line shape reduces to a simple sum of Voigt functions, as seen in Equation ( 6).This approximation is also necessary in our implementation of the FFM ion-dynamics approach (see Section 2.4).

Zeeman Effect
The effect of a static magnetic field B on the emission spectrum may be taken into account in the code by adding the term to the radiator Hamiltonian H er , where g s 2.0023193 is the anomalous gyromagnetic ratio for electron spin, µ B is the Bohr magneton, S and J are spin and total angular momenta.In the case of a combined Stark-Zeeman calculation, the Hamiltonian to be diagonalized thus reads H er = H 0 + H S + H Z , where H 0 and H S are respectively the unperturbed Hamiltonian and the Stark Hamiltonian, as introduced in Section 2.1.This requires the calculation of additional matrix elements of the type: The reduced matrix elements in the above expression are identical to those encountered in the calculation of magnetic-dipole transitions (M1).They are thus taken from the external atomic database, in a manner similar to the use of electric-dipole transitions (E1) for Stark effect.In principle, the code can take into account the Zeeman effect in the emission spectrum of any radiator ions (not limited to one-electron system) and for a magnetic field of arbitrary value (not limited to Paschen-Back nor weak-field limits), depending only on the ability to generate an appropriate atomic database.
A particularity of Zeeman effect is that the quantization axis is imposed by the direction of the static magnetic field, so that the emission is polarized in general.Therefore, the emission profile depends on the angle α between the direction of B and the line of sight of the observer: where In this notation, I q (ω) is the intensity calculated for a given polarization q = M − M which can take only 3 possible values: q = 0 (π polarization) and q = ±1 (σ ± polarizations).Moreover, the introduction of a preferential quantization axis breaks the spherical symmetry and the integration over the microfield becomes anisotropic and must be performed in the three directions of space [21,22]: where µ = cos(θ) is the cosine of the angle θ between the electric and magnetic fields.

Ion Dynamics Effects
The quasi-static approximation usually breaks down near the center of the line since this region corresponds to radiation over times long enough for the ions to move.Ion dynamics effects are implemented in the ZEST code using the fast Frequency-Fluctuation Model (FFM) [11].The time-dependent ion field fluctuations are treated as a collision-type process causing an exchange between the quasi-static radiative channels.This mixing is modeled by a stationary Markovian process using a single relaxation rate γ, characteristic of the ion field fluctuations.This rate is estimated by taking the ratio of the most probable ion velocity in the plasma (assuming a Maxwellian distribution) to the average ion-ion distance: where Z is the mean ionization of the plasma.The spectral profile accounting for the FFM ion dynamics method is given by the expression: with: In the above expression, I qs (ω) is the quasi-static line shape given by Equation ( 6).If we assume that the electron broadening operator does not depend on the frequency detuning (G(∆ω) ≈ G(0) i.e., impact approximation), then, from Equations ( 6) and ( 22), we note that the function J(ω) involves the convolution of Voigt functions (i.e., the convolution of Lorentzian and Gaussian functions) with the complex function 1  γ+iω .Taking the Fourier transform 3 of J(ω), and using the convolution theorem, we obtain: where Θ(t) is the Heaviside function.Taking the inverse Fourier transform of the last expression, we obtain: with where W(z) = e −z 2 erfc(-iz) is the complex Faddeeva function [23] which may be computed efficiently by different numerical methods.We note that the real part of W(z) is related to the Voigt function (the latter is often evaluated this way; see Appendix A).The Faddeeva function is frequently 3 We use an unconventional definition of the Fourier Transform in order to favor the domain t > 0.
encountered in the theory of the Brownian motion when considering the narrowing of the Doppler broadening due to collisions (Dicke effects) [24].In order to avoid the manipulation of complex numbers, we use an efficient and fast routine [25] which evaluates separately the real and imaginary parts of the Faddeeva function.Introducing the functions J R (ω) = Re J(ω) and J I (ω) = Im J(ω), the dynamic line-shape profile finally reads: In the case γ = 0, the last expression recovers the quasi-static limit as expected: For a very high fluctuation rate, the radiator system is unable to follow the ultra-fast modulations of the ion field, and the line spectrum is expected to tend approximately to the only electrons profile.However, it is known that the FFM approach does not reach this asymptotic form for γ → ∞ [26].In order to force the impact-approximation limit for both electrons and ions, it may be necessary to introduce a field-dependent fluctuation rate γ(F) or to use a semi-empirical amendment [8,27,28] to expression (20).One can also modify the FFM approach following the method described in reference [29].Another possibility is to use the Frequency Separation Technique (FST) [30], which separates the phase space into fast and slow perturbers.In the FST method, the line shape is approximated as a convolution of two profiles: a fast component based on the impact broadening theory, and a slow component for which we can use the FFM approach.

Microfield Distribution
There are many approaches to estimate the static ion fields distribution in a plasma [31], for instance: the expensive Molecular Dynamics simulations (MD) [32]; more affordable numerical models like the Adjustable Parameter EXponential approximation (APEX) [33] or the Monte-Carlo method (MC) [34]; simple analytic formulas with a restricted domain of validity like Holtsmark (for uncorrelated plasmas), Mayer (strongly correlated plasma) or Nearest Neighbour distributions (strong fields).We need a fast and accurate description of the ion microfield distribution which covers a broad range of plasma conditions.To achieve this goal, we use analytic fits of numerous Monte-Carlo simulations (MC) which reproduce P(F) with a great accuracy [35].The MC simulations were performed at a neutral or charged point using Debye quasi-particles interacting through screened Coulomb potentials.The screening length λ e accounts for the shielding of ion-ion interactions by a constant neutralizing electron background.In the MC simulations, a plasma configuration is represented by a collection of N + 1 interacting particles in a cubic box (one emitter and N perturbers) whose dimension L is related to the ion density N i = (N + 1)/L 3 .All perturbers are assumed to have the same charge Z p e.During the MC run, a sampling procedure is applied in the configuration space, by choosing randomly a particle in the box and its displacement.This displacement is limited to a maximum value which accounts for the plasma conditions (ion coupling parameter and screening length).The effects of distant particles are replicated by using periodic boundary conditions.The state of equilibrium is searched by the Metropolis algorithm, yielding a set of configurations of minimal energy from which it is possible to determine the electric fields distribution at the emitter.The microfield distributions obtained from numerous MC simulations covering a broad range of plasma conditions were then fitted with a reduced number of parameters including the ion coupling parameter Γ ii , the screening factor U ei and the reduced microfield variable β = F/F 0 (where F 0 is the normal field between two perturbers).These fits cover the domain: Γ ii ≤ 100 and 0 ≤ U ei ≤ 2.
Since the MC calculations use the same charge Z p e for all the perturbers (and possibly for the charged emitter), we assume that this charge is equal to the mean ion charge of the plasma: Z p = Z.In the average-atom model, the ion-coupling parameter reads: where Z is the mean ionization of the plasma, T i is the ion temperature and is the mean distance between two ions (Wigner-Seitz radius).The screening parameter is defined by: The electron-ion screening length λ e is evaluated in the Thomas-Fermi-Debye-Hückel approach [2], which makes a bridge between the classical description of Debye-Hückel and the Thomas-Fermi quantum theory: (32) where F n (x) = ∞ 0 y n dy 1+e y−x is the complete Fermi-Dirac integral of order n, η = µ k B T e is the reduced chemical potential.This allows one to tackle situations from hot and diluted plasmas (relevant to laser-plasma conditions) to cold dense plasmas (encountered in astrophysics), for which we recover the Debye-Hückel length and the Thomas-Fermi screening length at zero temperature, respectively (see the derivation of these asymptotic limits in Appendix B).
The microfield distribution P(F) is determined in the following manner.The plasma conditions (atomic number Z, T e , T i and N i ) are used to determine the mean ionization Z of the plasma in LTE or non-LTE, and the electron density N e = ZN i (alternatively, it is possible to fix N e and to determine N i ).The chemical potential µ is then calculated in the Thomas-Fermi approximation by solving the equation: We can now determine the parameters Γ ii , U ei and the normal field which read . The microfield distribution is then obtained from the analytical fitting formula P: There are two different formulas: one for a charged emitter (assuming the same charge as the perturbers, i.e., Z) and another one for a neutral emitter.
A comparison of the distributions obtained with different sets of parameters is shown in Figure 3.The left part addresses screened plasma ions (with U ei = 1), whereas the right part is for unscreened (U ei = 0) one-component plasma (OCP).As expected, the effect of the screening is to narrow the distribution and to shift it toward smaller field values.In the case of OCP, two extreme cases are recovered: the Mayer distribution for strongly-coupled plasmas and the Holtsmark distribution for uncorrelated plasmas.

Steps of the Calculation
The main input required for a line shape calculation with the ZEST code are: 1. the hypothesis for the calculation of state populations (SWAP, LTE, non-LTE) 2. the plasma parameters including: the atomic number Z, the ion temperature T i , the electron temperature T e , the electron (or ion) density N e (or N i ) 3. the value of the static magnetic field (where appropriate) 4. the grid for the sampling of the microfield distribution 5. the initial and final Layzer complexes to define the radiator 6. the spectral range on which the line shape is to be determined 7. the atomic database relevant to the considered case We recall that a Layzer complex is a special type of super-configuration where spin-orbitals n j having the same principal quantum number are gathered into a single shell n.The number of electrons is then fixed on average for each shell of the complex: where 0 ≤ p n ≤ 2n 2 are integers.The use of two Layzer complexes allows one to specify very simply how states will be mixed by Stark/Zeeman effects (indeed, we assume that only states belonging to the same complex may interact with each other) and which radiative transitions will be studied.For example: the Lyman-α line is defined by the complexes (2) 1 → (1) 1 ; the He-β line is defined by (1) 1 (3) 1 → (1) 2 .In general, two complexes may be linked by more than one mono-electronic transition (e.g., (3) 1 → (2) 1 which includes 3s → 2p, 3p → 2s and 3d → 2p transitions), so that we need to specify the spectral range of interest.The use of complexes implies that only ∆n = 0 interaction channels are taken into account in the Stark-Zeeman Hamiltonian.This hypothesis might be removed however by using super-configurations that gather the spin-orbitals of two (or more) shells, for example.
The first step is the generation of the atomic database required by the line-shape calculation.This is done using the MCDF [13] or cFAC codes [14,15].The following data are computed and stored in a file: the index, total angular momentum J, energy, parity and configuration label for the levels of both complexes; the signed reduced matrix elements for E1 and M1 transitions between levels of the same complex (relevant to Stark and Zeeman interactions, respectively); the signed reduced matrix elements for E1 transitions between levels of distinct complexes (relevant to the emission spectrum of interest).In the case of non-LTE plasmas, an external collisional-radiative model is used to determine the population of levels of both complexes.
The second step is the calculation of the microfield distribution, described in Section 2.5.We use the fit formulas, which depend on the ion coupling parameter Γ ii and the screening factor U ei , in order to sample the microfields on the predefined grid.
The third step is an initialization phase.We use the atomic database to retrieve the list of levels in each complex.These levels are split to obtain the list of unperturbed basis states for the initial and final complexes.Accounting for the symmetries of the system, we determine an optimal block-diagonal form for the radiator Hamiltonian, and we initialize the corresponding arrays.
In the fourth and final step, we do a loop over the microfield values.For each ion field, the Liouvillian matrix is built by retrieving all required matrix elements in the database, and diagonalized.The corresponding eigenvalues and eigenvectors are used to compute the strength and the energy of the Stark components, and the values for the homogeneous broadening.The population of the initial states are determined using either: the Boltzmann formula (in LTE), the non-LTE level populations of the CR code, or are assumed to be constant (SWAP assumption).We then compute separately two partial spectra, corresponding to the terms f k Re(F k ) and f k Im(F k ) in Equation (26), where F k depends on the complex Faddeeva function and the coefficients f k contain the populations of initial states, the strength of the Stark components and the weight of the current microfield.These contributions are accumulated to obtain the functions J R (ω) and J I (ω) from which we can derive the dynamic or static line shape profile.

Magnetized Plasmas
In Figure 4, we compare the profiles of the Hydrogen Lyman-10 line (transition n = 10 to n = 1), computed for different values of the magnetic field: 0, 10 T and 20 T. One-component electron plasma at T = 1 eV, N e = 10 13 cm −3 is assumed.These cases address typical tokamak edge conditions.The Doppler effect is omitted.The spectra are computed by assuming an isotropic magnetic field (sum of the three polarizations with equal weight).

Linear-Stark Effects
We first consider a hydrogen plasma at N e = 10 17 cm −3 and T = 100 eV.In this low density and high temperature regime, the plasma is in a classical limit with low values for the ion coupling parameter (Γ ii = 10 −3 ) and screening factor (U ei = 0.056).The microfield distribution is thus close to the Holtsmark distribution.The rate of fluctuation of the ion field is estimated to be γ = 54.8 cm −1 , which is quite high when compared to the half-widths of the quasi-static profiles of Ly-α (0.36 cm −1 ) and Ly-β (56.5 cm −1 ) lines.The comparison between the quasi-static line shape and the dynamic profile implementing the FFM ion dynamics effects is shown in Figure 5 for the Ly-β and Figure 6 for Ly-α.To ease the comparison, the Doppler effect has been omitted in the calculations (electron broadening widths were taken into account, however).As expected, the Stark components are merged together by the strong modulations of the microfield, yielding a single unresolved feature for each line.In the case of Ly-β, the dynamic profile is narrower than the quasi-static profile, which is a consequence of the Dicke effect [11].This narrowing is not observed for the Ly-α line because of the central component which is insensitive to the microfield intensity (it depends only on the direction of the fields).We note however that the FFM approach mixes the central component with the other Stark-shifted components, which is a known flaw of the method.Since the microfield fluctuation rate is highly dependent on the ion temperature, it is interesting to compare the dynamic line shapes calculated for H Ly-α and H Ly-β lines at N e = 10 17 cm −3 and various temperatures: T = T e = T i = 1, 10 and 100 eV.For a better insight into the ion dynamics effects, the Holtsmark distribution is used for all calculations and the homogeneous broadening has been omitted (no electron broadening and Doppler width).These only ions calculations thus differ exclusively by the ion dynamics effects on the inhomogeneous broadening.Results for H Ly-β are H Ly-α lines are shown in Figures 7 and 8, respectively.As expected, the dynamic profiles tend to Lorentzian functions as the ion temperature and the microfield fluctuation rate increase.In the case of H Ly-β line, the gap at the line center is filled gradually and we observe again a narrowing of the profile (Dicke effect).

Interpretation of Experiment
The chlorine Helium-β line and its Lithium-like satellites were measured recently on the ORION laser facility [36].The interpretation of the entire spectrum requires the use of a collisional-radiative model, because of the electron-excitation and dielectronic-recombination channels that are particularly important to reproduce the relative intensity between the Li-like structures and the He-like resonance line.The conditions inferred in the publication are around T e = T i = 500 eV and N e = 10 22 cm −3 .However, the shape of the Helium-β line, alone, is relatively insensitive to the assumption on the populations: we can thus use a LTE assumption to simplify the analysis of this line (this is not the place to discuss non-LTE physics).
The line shapes of the Chlorine Helium-β line were computed with the ZEST code at a temperature 4 T i = T e = 440 eV for electron densities varying from N e = 10 21 cm −3 to N e = 5.10 22 cm −3 .The calculations are compared to the normalized experimental profile in Figure 9.The Stark calculations include the ions dynamics effects and a homogeneous broadening including the electron-collision broadening, Doppler and natural width.The atomic database was generated with the cFAC code for the complexes (1) 1 (3) 1 → (1) 2 (a 0.5 eV shift was added to the energy of the 3p − 1s lines).An experimental Gaussian FWHM of 2 eV was also taken into account in the Voigt profiles.As mentioned in the paper, this experimental width is not consistent with the resolving power of the OHREX spectrometer (10,000 equivalent to FWHM=0.3 eV), but this is assumed to be the consequence the time-integration over the full plasma evolution.As we can see, the width and shape of the He-β line is very sensitive to the electron density.In particular, our analysis took advantage of the strong influence of the Stark effect on the relative intensity between the structures at hν = 3266.8eV and hν = 3271 eV, which correspond roughly to the position of the unperturbed lines 1s + 3p − → 1s 2 and 1s + 3p + → 1s 2 , respectively.As the electron density increases, the states in the complex (1s) 1 (3s3p3d) 1 are mixed together more efficiently, resulting in a progressive redistribution of strengths among the 4 The temperature T = 440 eV was taken from a more detailed analysis of the He-like and Li-like spectra using the non-LTE code DEDALE [10] and the present line-shape code ZEST (see Acknowledgments).lines.For a very large electron density, the off-diagonal Stark matrix elements of the Hamiltonian may become even more important than the effects of the fine-structure on the diagonal matrix elements; when this limit is reached, the corresponding profile becomes symmetric.As we can see, this limit is almost reached at 5.10 22 cm −3 .The best agreement with the experimental data is obtained for densities between N e = 5.10 21 cm −3 and N e = 10 22 cm −3 .

Discussions and Perspectives
A number of assumptions have been made to make the code as fast as possible, so that it can be inlined in an opacity or emissivity model.Some of these approximations need to be relaxed, revisited or investigated.This concerns, for example, the improvement of the calculation of the trace and the inversion of the complex matrix A = ω − iL(F) + φ(ω) [37,38].Taking into account the interference terms in the electronic collision operator may be also important, as shown in [29] for lithium-like satellites of argon.It may also be interesting to study different approximations for the collision operator [39], especially in the limit of weakly coupled plasmas [40] or when accounting for a magnetic field [41].Another interesting development is the implementation of penetrating collisions which can be taken into account analytically in the case of hydrogen [42] and hydrogen-like ions [43] through an ad hoc theory.For the moment, the code is limited to the dipole approximation for the Stark effect.However, recent studies have shown that quadrupole corrections can be important in some cases [44].Accounting for such effects requires the estimation of spatial gradients of the microfield distribution.
Another point that we need to improve in the code is the static limit of the profile.Indeed, the present version of the code assumes the impact approximation (constant electron broadening width for all frequencies) because it allows us to use Voigt functions which ease the implementation of the FFM ion-dynamics approach (leading to simple analytical expressions based on the Fadeeva function).However, this fast procedure does not recover the static limit as expected in the far wings (I(ω) ∼ ω −5/2 ω −2 instead of I(ω) ∼ ω −5/2 ).This issue may be solved easily by using a frequency-dependent electron broadening operator (such as the one in Section 2.2) and by doing the convolutions numerically (e.g., by Fast Fourier Transform).This is however computationally much more expensive than what is done now.
Finally, the FFM ion-dynamics approach provides an excellent route to fast and accurate line-shape calculations, but it is not valid for ultra-fast modulations of the ion field.Different methods have been proposed to extend the validity of the FFM approach and we plan to investigate them.
If we now consider plasmas in the highly-degenerate limit (low temperature and high density), the Fermi-Dirac integrals take the following expressions, Introducing the latter expression in Equation (32), we obtain the Thomas-Fermi screening length at zero temperature:

Figure 1 .
Figure 1.The frequency-dependent part of the electron broadening operator is plotted as a function of the reduced detuning with respect to the center of the line.The dimensionless quantity is evaluated using three different models and approximations (Dufty, GBK, Lee: see text) for κ m λ D = 30.

Figure 2 .
Figure 2. Discrepancies on the electron impact widths at the line center (Lee versus GBK), as a function the parameter κ m λ D which represents the ratio of maximal over minimal impact parameters for the electron-radiator collision.

Figure 3 .
Figure 3. Microfield distributions at a charged emitter are calculated as a function of the reduced field for different values of the ion coupling parameter Γ and screening factor U, using the fit functions of Ref [35].

Figure 4 .
Figure 4. Hydrogen Lyman-10 line shapes at N e = 10 13 cm −3 and T e = 1 eV for different values of the static magnetic field: B = 0, 10 and 20 T. The photon energy is relative to the center-of-gravity of the profile.The Doppler effect is omitted in the calculations.

Figure 5 .
Figure 5. Lyman-β line shapes of hydrogen at T = 100 eV and N e = 10 17 cm −3 .The calculations were performed with the ZEST code either in the quasi-static approximation (QS) or accounting for ion dynamics effects (FFM), assuming a Holtsmark microfield distribution.The Doppler effect has been omitted for the comparison (electron broadening only).

Figure 6 .
Figure 6.Same as Figure 5 for the Lyman-α line shapes.

Figure 7 .
Figure 7. Lyman-β line shapes of Hydrogen at N e = 10 17 cm −3 and various temperatures T = 1, 10 and 100 eV.The calculations were performed with the ZEST code accounting for ion dynamics effects (FFM) and assuming a Holtsmark microfield distribution.Doppler and electron broadenings have been omitted for the comparison.

Figure 8 .
Figure 8. Same as Figure 7 for the dynamic Lyman-α line shapes of Hydrogen.

Figure 9 .
Figure 9. Analysis of the Chlorine He-β profile measured on the ORION laser facility [36].The ZEST calculations are performed for a LTE plasma at T i = T e = 440 eV and various electron density N e .An experimental Gaussian FWHM of 2 eV is taken into account in the simulations.