ElecSus: Extension to arbitrary geometry magneto-optics

We present a major update to ElecSus, a computer program and underlying model to calculate the electric susceptibility of an alkali-metal atomic vapour. Knowledge of the electric susceptibility of a medium is essential to predict its absorptive and dispersive properties. In this version we implement several changes which significantly extend the range of applications of ElecSus, the most important of which is support for non-axial magnetic fields (i.e. fields which are not aligned with the light propagation axis). Suporting this change requires a much more general approach to light propagation in the system, which we have now implemented. We exemplify many of these new applications by comparing ElecSus to experimental data. In addition, we have developed a graphical user interface front-end which makes the program much more accessible, and have improved on several other minor areas of the program structure.

Development of computational tools such as ARC [13], The Software Atom [14] and ElecSus [15] plays an important role in the development of these applications, where system parameters can be optimised in theory then tested and verified experimentally.
For example, understanding the absorption and dispersion of an atomic vapour has led to a deeper understanding of atomic filters based on the Faraday effect, and modelling this has led to the optimisation of the linewidth of such filters [16], or optimisation of the filter in the presence of homogeneous broadening [17]. One can then use these optimised filters in other applications, e.g. using the filter to make an intrinsically frequencystable laser system [9], creating a dichroic beam splitter for Raman light [18] or filtering frequency-degenerate photon pairs from an optical parametric oscillator [19].
The previous version of ElecSus has been used in a wide range of experiments, including magnetic field imaging [20,21], Faraday filtering [16,17,9,22,23], characterisation of hybrid atom-cavity systems [24], determination of spin polarization of optically pumped atoms [25] and absolute absorption measurements [26]. Since the first publication of ElecSus in 2015 [15], we have added significant functionality that adds to both the scientific scope and the accessibility of the program. In brief, these are: • Adding support for magnetic fields that are not parallel to the light propagation axis (i.e. non-Faraday geometry).
• Directly calculating the propagation of electric fields via Jones calculus, which allows, for example: magnetic field gradients across the atomic medium; simulating imperfect polarisers; and simulating birefringent optics.
• A graphical user interface (GUI) now allows users with no knowledge of computer programming to use the majority of the program features.
• A rewrite of the fitting methods using the lmfit module [27], which allows the user to impose bounds on each of the fit parameters. In addition, the user can now select the differential evolution fitting algorithm which is an efficient global optimisation routine for multiparameter fits.
• Several minor changes and bug fixes to program operation since original publication In the remainder of the paper, we discuss the physics and computational implementation of the above additions.

Recap of important concepts
The majority of the theoretical background is unchanged from the original ElecSus publication [15]. However, we briefly summarise here the general principles and important equations.
The main effort of the program is to calculate the linear electric susceptibility, χ, of an atomic vapour exposed to a near-resonant weak-probe laser field. In the weak-probe limit [28,29] optical pumping can be neglected and the optical properties of nearby transitions can be treated independently. A single transition, labelled i, is treated as an isolated two-level atom, where neglecting the atomic motion, the susceptibility is given by where d 2 is the reduced dipole matrix element of the transition, C 2 i is the transition strength, Γ is the natural linewidth of the transition, ∆ i = ω − ω i is the difference between the laser angular frequency ω and the resonance frequency ω i = (E e − E g )/ , which is the difference in energy between a ground state |g and excited state |e divided by the reduced Planck constant. N g is the atomic number density of a particular state in the ground manifold, which in thermal equilibrium is given by where N is the total atomic number density, F a is the isotopic fraction, k B is Boltzmann's constant, and T is the absolute temperature in Kelvin. The sum is over all of the nuclear substates, labelled according to the nuclear spin quantum number I and its projection m I . The energy difference ∆E j is measured with respect to the lowest energy in the ground manifold. The fractional weighting of each of the 2(2I + 1) ground states is nearly uniform, i.e. 1/(2(2I + 1)), for most cases, except where the temperature is very low or the energy difference between ground states becomes very large, for example in extremely high ( 1 T) magnetic fields. This fractional weighting is a new addition since the original version of ElecSus [15]. Atomic motion causes an inhomogeneous broadening of the bare atomic lines; the atoms experience a Doppler-shifted frequency according to their component of velocity in the direction of the beam v z , which is a 1D Maxwell-Boltzmann distribution given by equation (6) in the original publication [15]. The resulting atomic lineshape is thus a convolution between the stationary atomic response (the Lorentzian f (∆)) and the Gaussian distribution of velocities g(v), which yields the well-known Voigt profile, In a multi-level system, the total susceptibility at a given global frequency detuning ∆ can be found by summing over each transition, We use a matrix representation of the atomic Hamiltonian in the |m L , m S , m I quantum number basis (m L,S ,I are the quantum numbers associated with the projection of the electronic orbital, electronic spin and nuclear spin angular momenta, respectively) to calculate the resonant frequencies and transition strengths, as described in section 2.2 of the original publication [15]. The Hamiltonian includes details of the internal level structure such as fine and hyperfine structure, and interactions with an external magnetic field (Zeeman effect). A separate Hamiltonian is calculated for the ground state nS and excited state nP; the Hamiltonians are diagonalised to find eigenenergies E j and eigenstates | j , which are in general a superposition of basis states. The transition frequencies are the difference in energy between a ground state |g and excited state |e . The transition strength is calculated from the dipole matrix element | g|er q |e | 2 , where the subscript q denotes the component of the dipole operator in the spherical basis, and denotes the type of electronic transition; π, σ + or σ − which are associated with a ∆m L = 0, +1 or −1 transition, respectively. In the previous version of ElecSus [15] the π transitions were not calculated since, in the Faraday geometry where the magnetic field vector is parallel to the wavevector of the light (B ·k = ±1, where the hat denotes the unit vector,B = B/|B|), π transitions are forbidden [30]. In this version we relax the constraint on the magnetic field geometry, and therefore we additionally calculate the π transitions.
3. Electric field propagation in an atomic medium with an applied magnetic field The major improvement over the original version of Elec-Sus (versions 1 and 2) is relaxing the constraint that the magnetic field axis must be parallel to the light propagation axis (i.e. the Faraday geometry). In making this change, we must also consider how the polarisation of the input light couples to the atomic transitions, and the resulting effect this has on light propagation through the atomic medium. In this section we will present the general approach to this problem, and then the special cases of the Faraday and Voigt geometries. This section follows the work of Palik and Furdyna [31], and Rotondaro, Zhdanov and Knize [32].

The wave equation and the dielectric tensor
We start with Maxwell's wave equation for a non-magnetic dielectric medium, which is given by where k is the wavevector, E is the electric field, ω is the angular frequency of the plane-wave and is the dielectric tensor. We assume that the light is a plane-wave that propagates in the z-axis; the electric field therefore lies in the x − y plane. The applied magnetic field B can take any angle, but it is practically easiest to rotate the coordinate system around the z-axis such that the magnetic field lies in the x−z plane, which is effectively just a polarisation rotation in the x−y plane by an angle φ B . The magnetic field then makes an angle θ B with the z-axis, as shown in figure 1.
After defining the complex refractive index of the medium n, Figure 1: Representation of the geometry of the system under consideration. θ B is the angle that B makes with the z-axis, φ B is the angle between x and the projection of B into the x − y plane, and the arrow-heads indicate the sign convention. θ B = 0, φ B = 0 yields the Faraday geometry, while θ B = π/2 or φ B = π/2 becomes the Voigt geometry. in the coordinate system discussed above, the wave equation can be written in matrix form as [32] where x , xy , and z are elements of the dielectric tensor, and are related to the complex electric susceptibility by and the χ +,−,0 are the susceptibilities associated with σ + , σ − and π transitions, respectively. The two non-trivial (i.e. |E| 0) solutions of eq. (9) are found by setting the determinant of the matrix to zero, which results in a quadratic equation in n 2 . Each of the two solutions n 1 and n 2 can then be substituted in to find a zero-value eigenvector e 1,2 which together represent the (orthogonal) normal modes for propagation of light in the system, which are dependent on θ B .
To calculate the transmitted field, one must transform the x, y, z basis into the normal mode basis, via the rotation matrix where the * denotes the complex conjugate, and after which propagation in the z-axis through a distance L is computed by evolving each normal mode in space with its associated refractive index n 1 and n 2 . This is done via the diagonal propagation matrix Finally, one can use the inverse matrix M −1 to transform back to the cartesian coordinates. In terms of matrix operations, and including rotation around the z-axis, the full set of operations is where for shorthand we combine all these processes into a single effective Jones matrix J χ [33]. For arbitrary angle θ B , the solutions for n 1 and n 2 and their associated eigenvectors are not easy to write down analytically. In this case, we use the symbolic Python (sympy) package to solve these equations. However, for the case where either θ B = 0, π (the Faraday geometry) or θ B = π/2, 3π/2, there are analytic solutions for n 1,2 and e 1,2 which are much more computationally simple to implement. The Faraday geometry was the only case that the original version of ElecSus accounted for. In the next subsections we will describe these special cases in more detail.

The Faraday geometry
The Faraday geometry, wherek ·B = ±1, is the geometry in which the Faraday effect [refs] is observed. In this geometry, we find the solutions for the refractive index are given by ( for θ B = 0) or alternately, that the two indices are associated with σ + and σ − transitions. The corresponding eigenvectors are It is clear in this case that applying the rotation matrix M to the x, y, z coordinate system simply transforms the coordinates into the circular basis, and we find as expected that σ + transitions couple directly to left circularly polarised light, and σ − transitions couple directly to right circularly polarised light. No component of the light couples to π transitions in this geometry. Switching θ B = 0 to θ B = π simply inverts the coupling, i.e. n 1 and n 2 are swapped, resulting in the circular polarisations and σ ± transitions coupling to each other in the opposite way.
In the general case, there is a different refractive index for each of the two circular polarisations (circular birefringence and dichroism), which on propagation leads to a rotation of the plane of polarisation (Faraday roatation). Light that is initially linearly polarised can be decomposed into the cirular basis (in the x − y plane), and on propagation if the medium is resonant with only one of the σ ± transitions, one component of circular polarisation will be absorbed, leading to an effective circular polarisation filter -the output will be just the circular component that is not absorbed.

The Voigt geometry
The Voigt geometry, wherek ·B = 0, is the geometry in which the magnetic field axis is transverse to the light propagation axis. The solutions for the refractive index in this geometry are [32] therefore n 1 is associated with both σ ± transitions, while n 2 is associated with only π transitions. The corresponding eigenvectors are (for θ B = π/2) For e 2 it is clear that the field component that is parallel to the magnetic field (the x-axis in the case where φ B = 0) drives π transitions. The first eigenvector is harder to immediately visualise; we know that the eigenvectors must be perpendicular to both each other and k, and hence e 1 must point along the cartesian axis perpendicular to both k and B (i.e. it points along y in the case where φ B = 0). The normal mode is elliptically polarised in the plane perpendicular to B, i.e. in the y − z plane. Since the atomic quantisation axis lies along B, a linearly polarised beam along y can be decomposed, in the atomic frame, into equal circular components in the y − z plane, and thus it drives (equally) both σ + and σ − transitions.
In contrast to the Faraday case where the medium exhibits circular birefringence and dichroism, in the Voigt geometry the system is linearly dichroic and birefringent, which leads to a different form of magneto-optic rotation known as the Voigt effect [34][35][36]. Note that, as pointed out by Pershan [37], the Voigt effect is subtly different to the Cotton-Mouton effect [38,39], where birefringence emerges as a result of alignment of diamagnetic molecules in a transverse magnetic field.

Jones Matrices for propagating fields
In the previous section, the use of matrices is a simple and computationally convenient method for calculating the propagation of the electric field in the medium. These matrix methods are generally referred to as Jones calculus [33], and similar matrices can be constructed for a variety of common optical elements, including waveplates and polarisers (see table A.1 in the appendix). Using these matrices, one can calculate the output electric field (and intensity) from any arbitrary combination of optics.

Stokes parameters using Jones Matrices
As in the previous version, the four Stokes parameters S 0,1,2,3 are an easily measurable way of characterising the polarisation state of light. The Jones calculus approach offers an intuitive way of calculating these parameters, and only require calculation of the output field. The four Stokes parameters are S 0 ≡ (I LCP + I RCP )/I 0 = (I x + I y )/I 0 = (I + I )/I 0 , (20a) In the above equations, S 0 is the normalised output intensity and requires no extra matrices -all that is required is to sum the mod-squared field components (I j = ε 0 c|E j | 2 /2, with jx, y, LCP, RCP) in whichever orthogonal basis is most appropriate. Each of the components E x , E y , E LCP , E RCP , E and E can be found by applying the respective Jones matrix to the output field -for example, which are required to compute the S 1 parameter, the difference between the linear basis components in the x and y plane. S 2 is the difference between orthogonal linear polarisations at 45 degrees to the x and y axes, and S 3 is the difference between the two circular polarisation components (note the change in notation from the previous publication -in a non-Faraday geometry,

New applications of ElecSus
In this section we illustrate the new features of ElecSus through a series of example data sets. Note that the examples here are all for Rb, but ElecSus will work with the alkali-metal atoms Rb, Cs, K and Na.

Transmission spectroscopy in the Voigt geometry
Transmission spectroscopy in the Voigt geometry is the simplest method to demonstrate the addition of π transitions to ElecSus. We perform weak-probe [29] transmission spectroscopy (intensity approximately 0.03 mW/cm 2 ) through a 1 mm naturally abundant Rb vapour cell. The cell temperature was approximately 95 • C, and the applied magnetic field strength was approximately 0.42 T, along the x-axis. The laser source is a distributed feedback (DFB) laser (quoted linewidth approx. 2 MHz), which is frequency tuned by temperature tuning of the diode. This allows a mode-hop-free laser scan over an extremely large detuning range (up to ∼ 1 THz). The scan is linearised with a Fabry-Perot etalon, and a room temperature 75 mm reference cell is used as an absolute frequency reference, following the procedure outlined in reference [40]. Tuning the diode temperature slightly changes the laser output power, so we stabilise the optical power that is incident on the atomic vapour using a feedback loop linked to the RF power of an AOM, as described in ref. [41]. Figure 2 shows the resulting transmisison spectra for 3 different angles of incident linear polarisation. For a linearly polarised input electric field, we define θ E as the angle the electric field makes with the x-axis for simplicity. Data are shown as purple points. For light polarised with the electric field along x (θ E = 0, top panel), E B and therefore the only allowed transitions are π transitions. When the electric field oscillates in the y-axis (θ E = π/2, third panel), E ⊥ B and hence σ ± transitions are driven. Finally, when the electric field is at 45 degrees to the x, y axes (second panel), all three types of transition are driven. At this magnetic field strength, for naturally abundant Rb, the spectra in all three cases are very complex. However, we see in all three panels the fit to the data using ElecSus, which are in excellent agreement in all 3 cases (the RMS errors are 0.64%, 0.36% and 0.59% for the θ E = 0, π/4, π/2 data, respectively). We also plot the residuals (difference between theory and experiment) for the θ E = π/2 data, multiplied by a factor of 100, on the bottom panel; the lack of any clear structure in the residuals indicates that the theoretical model incorporates all of the underlying physics.

Stokes polarimetry in the Voigt geometry
In the absence of an applied magnetic field, or in a small applied field, Doppler broadening masks the complex atomic structure in a thermal vapour. Is is therefore more instructive to consider the case of a large applied field. This regime, known as the hyperfine Paschen Back (HPB) regime, has generated much recent interest [42-46, 8, 47-49]. For 87 Rb in the HPB regime, due to the large Zeeman splitting, the atomic transitions are all separated by more than the width of the Doppler broadened lines. For some applications with thermal vapours, this greatly simplifies the physics, since true 2-, 3-, and 4-level systems can be easily isolated, allowing for archetypal demonstrations of selective reflection [50], electromagnetically induced transparency [51], electromagnetically induced absorption [52], and four-wave-mixing processes [53,54].
In figure 3 we plot theoretical predictions of the Stokes parameters for a Rb D2 line spectrum (isotopically enriched, 99% 87 Rb) with an applied magnetic field of 1.54 T, in the voigt geometry (θ B = π/2, φ B = 0). The vapour cell length is 1 mm, cell temperature is 120 • C and the input polarisation is linear, where the T denotes the matrix transpose). The choice of magnetic field, vapour cell length and isotopic composition in figure 3 was cho-sen to match an experimental setup available to us, which will be described later.
From the S 0 spectrum ( fig. 3(a)), we can observe 6 main (strong) sets of features in groups of 4, and two sets of visible smaller (weak) features on the far edges of the spectrum. The group of 4 comes from the projection of the nuclear spin quantum number m I , which for 87 Rb (nuclear spin I = 3/2) can take four values: -3/2, -1/2, 1/2 and 3/2. The two inner-most groups of four originate from π transitions (m J = ±1/2 → m J = ±1/2), whilst the outer groups are from σ − transitions on the side of negative detuning (m J = 1/2 → m J = −1/2, and m J = −1/2 → m J = −3/2) and σ + transitions on the side of positive detuning (m J = 1/2 → m J = 3/2, and m J = −1/2 → m J = 1/2). The weak transitions at ∼ ±60 GHz stem from the incomplete decoupling of the ground state 5S 1/2 into the |m S , m I basis (m L = 0 for the all terms in the ground state manifold) -the ground states are not pure eigenstates in this basis and there is therefore a small admixture of other states which results in the weak transitions, as described in ref. [8] for the Faraday geometry. Figure 4: Electric field propagation example. Each optical element, including the vapour cell, can be described by a Jones matrix which can be cascaded to calculate the output of any series of optical elements. This example shows a typical experimental setup to measure the S 1 Stokes parameter: the initial input polarisation E in is rotated by the half-waveplate (λ/2), passes through the cell including the birefringent windows (BRW), after which the field components E x and E y are analysed by the transmission and reflection (respectively) through a polarising beam splitter cube (PBS).

PBS BRW BRW
The S 1 spectrum ( fig. 3(b)) shows the consequence of the linear dichroism of the medium. Off resonance, there is no interaction and hence I x = I y = I 0 /2 and S 1 = 0. At the atomic resonance frequencies, one component of light in the linear basis is completely absorbed, leading to either I x,y = 0 while the other component far off-resonance and therefore unaffected. The S 1 spectrum therefore swings between values of ±0.5 depending on which component is absorbed.
The S 2 spectrum ( fig. 3(c)) is the difference between the linear polarisation components but in the diagonal (I − I ) basis. This spectrum therefore has an off-resonance value of 1 since we input the polarisation. On resonance the absorption dominates, which reduces both I and I and the S 2 value tends towards zero. Between the resonances, however, there is still optical rotation (since neither nor are eigenmodes of propagation); this is most pronounced between the π and σ ± groups at ±15 GHz detuning since the optical rotation adds.
Finally, the S 3 spectrum ( fig. 3(d)) is the Voigt-geometry equivalent of a Faraday rotation spectrum (which are observed in the S 1,2 Stokes parameters), which is intuitively what one might expect since in the Faraday case the medium is circularly birefringent and the polarisation rotation is observed in the linear basis, whilst the Voigt case is the opposite way around -the medium is linearly birefringent and the rotation can therefore be observed in the circular basis.

Modelling cell window birefringence
An issue in experimental polarimetry measurements comes from birefringence in optical elements, which causes unwanted additional optical rotation. In thermal vapour experiments, a common source of this unwanted birefringence is the vapour cell windows, but the amount of birefringence is not usually known a priori. We can model the effect of an unknown birefringent material through a Jones matrix, J BRW (φ BR , θ BR ) (see Appendix for details) which is characterised by 2 parameters, the phase shift φ BR and the orientation with respect to the optical axis θ BR (the subscript BR is used to differentiate these two parameters from the magnetic field angles used earlier). Figure 4 illustrates an example situation. An input electric field E in travels through a half waveplate (λ/2), then through a vapour cell containing an atomic medium of length L. The vapour cell has birefringent windows (BRW) on both ends. After the vapour cell, the polarisation state is analysed by placing a polarising beamsplitter cube (PBS) which analyses the x and y components of the electric field. The Jones calculus approach is to cascade the matrices for all these optical elements. Combining all elements, the outputs E x,y after the PBS cube would be Taking the difference in intensity between the two output ports of the beam splitter, we measure the S 1 Stokes parameter. As a demonstration of the effects of window birefringence, we now show some experimental polarimetry data. Our optical setup is the same as shown in figure 4. The applied magnetic field is produced by a 'magic sphere' configuration of NdFeB permanent magnets [55] which yields a peak field of 1.54 T at the centre of the hollow cylinder, where we place a microfabricated vapour cell [56] inside a small copper heating block. Right-angled prisms allow the light to propagate through the cell at normal incidence to the field, realising the Voigt geometry.
The S 1 spectroscopic data are shown in figure 5. The experimental conditions are similar to those of figure 3, except the cell temperature which is 98 • C. The input polarisation is set so that far off-resonance, the difference signal is zero. In the absence of window birefringence, this would mean an input polarisation angle θ E = ±π/4. However, when including window birefringence, the input polarisation needs to be rotated by a small amount to satisfy this condition. The data (purple points) show the main optical rotation features displayed in the S 1 spectrum of figure 3, except in the wings of the resonance lines. This is most prominent at around ±15 GHz detuning, where the birefringence of the cell windows causes an additional rotation. On the figure we plot two theoretical curves. The dashed olive line is the theory without birefringence, whilst the blue solid line is the result of a fit to the data, assuming both windows have equal birefringent properties, and allowing the birefringence parameters φ BR , θ BR and the input linear polarisation angle θ E to vary, along with the cell temperature and magnetic field strength. We find excellent agreement with the experimental data, as demonstrated by the small residuals on the bottom panel of figure 5 -the RMS error between theory and experiment is 0.8%. This may therefore be a useful technique to practically determine the birefringent properties of such windows.
In some cases the birefringence can be compensated for with the addition of waveplates to the optical setup. In figure 6 we show an experimental S 3 spectrum, with the same conditions as figure 5 apart from the cell temperature, which is 125 • C. The optical setup is similar to figure 4, with the addition of a halfwaveplate and a quarter-waveplate between the vapour cell and the PBS cube. The quarter waveplate and the PBS constitute a circular polarisation analyser, which is the usual experimental technique for measuring S 3 . The extra half-waveplate is used to compensate for the cell birefringence which would otherwise offset the rotation signal, resulting in a spectrum that is qualitatively very similar to that of figure 3 (the expected spectrum without birefringence is the olive dashed line in figure 6 for direct comparison).
For the fit to this dataset, we constrain the cell birefringence properties to be those from the fit in figure 5, and instead allow the angle of the half-and quarter-waveplates after the cell to vary. In this case, we again find excellent agreement between the data and the fit, with RMS residuals of 1.1%.

Arbitrary angle geometry spectroscopy
For the most general case where θ B 0, π/2, the spectral features are very rich. In a study of magneto-optic filtering, Rotondaro, Zhdanov and Knize [32] showed that the bandwidth of atomic filters can be reduced by using a non-standard magnetic field geometry. Again, here we present a comparison of ElecSus to an experimental data set to illustrate its use in these situations. The experimental setup utilises a 1 mm natural abundance Rb cell, with an applied magnetic field provided by two permanent top-hat shaped magnets set up on a rotation platform, such that a range of angles can be formed between the light propagation axis and the magnetic field axis, limited only by the radial extent of the magnets and their mechanical mounts. The maximum field strength is limited to around 0.4 T with this setup. In figure 7 we show an S 0 spectrum with a magnetic field strength of 0.42 T at an angle θ B ≈ π/3. The experimental data are shown as purple points, and the fit to the data using Elec-Sus is shown as a blue line. Clearly, experiment and theory match very well, as shown by the residuals. The RMS error between theory and experiment is 0.8%. Since the field strength is similar to the data in figure 2, the atomic resonance positions are in nearly the same place. We can then compare for similar conditions (i.e. incident polarisation angle) which transitions are driven. The central features within approximately ±7 GHz originate from π transitions, so there is similarity to the top panel of fig. 2. However, the features at larger detuning in fig. 7 come from σ ± transitions and are therefore not present when interrogating the medium with θ E = 0 polarised light in the Voigt geometry. The relative strength of the atom-light coupling is also completely different to any of the Voigt-geometry cases; for the data in figure 7 the π transitions are clearly more strongly driven than the σ ± transitions.
This data set demonstrates that ElecSus is now suitable for use with arbitrary angle magnetic fields. We also note that we can succesfully reproduce the plots from reference [32]; we provide a Python script to reproduce these figures in the tests/ subdirectory of the GitHub repository.

Magnetic field gradients
Since the electric field is now calculated explicitly, and can be returned directly by the program, it is possible to use ElecSus to now predict spectra from non-uniform systems -this could be, for example, the magnetic field gradient across a thermal vapour cell. We place a room temperature, 75 mm naturally abundant Rb vapour cell between two top-hat magnets, which are separated by a little more than the vapour cell length. This creates an axial magnetic field profile (the Fardaay geometry) shown in the top panel of figure 8, which has calculated minimum/maximum/mean values inside the vapour cell of 41 mT, 311 mT and 100 mT, respectively. The middle and bottom panels show transmission spectroscopy and the Faraday rotation signal S 1 in this experimental configuration. Purple points are experimental data. The dashed grey lines show a calculation which assumes the mean value of magnetic field, whilst the blue lines show a calculation which splits the cell into, in this case, 25 segments, and propagates the electric field through each segment sequentially. The two models are very clearly different in both S 0 and S 1 spectra. We fit using this calculation; the fit parameters are the vapour cell temperature, and the position of the two top-hat magnets (relative separation and position offset with respect to the cell; the magnets' remnant field strength is fixed). There is excellent agreement between this model and experimental data (the RMS error is 0.2%).

Faraday filter with field gradient
Though the previous example showed the case of an extreme field gradient, it is not likely to be a practically relevant case. In this subsection we consider the application of field gradients to Faraday filtering. In reference [9], a Faraday filter with high peak transmission was demonstrated with an approximately uniform axial magnetic field, produced by placing the 5 mm vapour cell between two NdFeB ring magnets, separated by a large distance compared to their extent, and the extent of the vapour cell. For applications development, using a smaller, single magnet system to generate the field is attractive for mechanical simplicity, miniaturisation purposes and cost-saving. However, using a single small magnet necessarily means that the field profile becomes non-uniform. However, as we show in figure 9, we can effectively compensate for this field gradient with a suitable design of magnet, and ElecSus can be used as a design tool to optimise magnet specifications for this kind of application. In figure 9, we simulate the filter profile that could be achieved with a small ring magnet placed close to the vapour cell. The magnet parameters were found by allowing the dimensions of the magnet (inner, outer diameter and thickness), the separation from the cell, and the cell temperature to vary, subject to some upper bounds on the magnet dimensions. For each iteration, the magnetic field profile over the cell and the resulting Faraday filter transmission calculated. The Faraday filter was then optimised for peak transmission at line centre using the methods outlined in refs. [57,17]. The simulation shows that a similar filter profile is generated from the non-uniform field (solid blue line in fig. 9: field maximum/minimum over the cell: 222 G / 178 G), when compared with the filter profile used in ref. [9] (dashed line in fig. 9). Most of the features remain, and the peak transmission of the filter is slightly higher than the uniform field.

Program structure
The significant feature changes in ElecSus necessitated some changes to the overall program structure. Figure 10 shows a diagramatic illustration of information flow with the new program structure. The program can be accessed from either the GUI or an external Python script. In either case, the user supplies the simulation parameters, a set of exeperimental data if a fit is required and, optionally, the boundaries on fit parameters. These are passed to the calculate() or fit data() routines in elecsus methods.py, which calculate spectra by finding the energy levels and state vectors of the system Hamiltonian (EigenSystem), calculating the propagation of the electric field (SolveDielectric) and finally applying the relevant Jones matrices to find the transmitted fields and intensities. The fitting methods have been updated and now use the lmfit module [27] which allows the use of boundaries for fit parameterssee below for more details.

Graphical User Interface
A graphical user interface (GUI) was developed for ElecSus which makes using the program much more accessible, particularly for users without knowledge of programming -it is now possible to use most of the program features without using any of the back-end source code. In figure 11 we show a screenshot from the GUI, which is broadly split into two panels. On the left side, an interactive (i.e. axes can be dynamically rescaled) plot panel (using matplotlib [58]) shows the spectral data that has been already calculated or loaded from user-supplied csv files. If fitting has been performed, a second tab in this panel 7UDQVPLVVLRQS 0 'HWXQLQJ*+] 100R Figure 7: Example of arbitrary angle geometry spectroscopy. The magnetic field (strength B = 0.42 T) is oriented at an angle θ B ≈ π/3 from the z-axis, and we take an S 0 spectrum through a 1 mm naturally abundant Rb cell at a temperature T = 90 • C, with linearly polarised input light (θ E ≈ 0). Some similarity can be noted between this data and that of figure 2 (the magnetic field strengths are approximately equal), but there are clear differences, noably around ±7 GHz, from the Voigt geometry with the same input polarisation. The purple points are experimental data, and the blue line is a fit to the data using ElecSus, allowing θ E , θ B , B and T to vary.
shows the result of that fit, plotting residuals between the experimental and theoretical data. Additional tabs show text-based information about fit parameters, program status and any error messages that may have been generated. On the right side of the window are the program input parameters, with two tabs for purely theoretical calculations or fitting data. At the top of the panel, the user may select which outputs are displayed, from a list that includes all Stokes parameters, relative intensities of linearly polarised and circularly polarised light, and a few others. Underneath this are the input paramters. Figure 11 shows the theory calculation panel, and figure 12 shows the fit panel. For both, parameters are sub-divided into general parameters, and parameters specific to the magnetic field, and polarisation parameters. Any fit parameter can be allowed to vary or be held constant, selected via the "Float?" tick-box. On selection, the fit bounds options become active, allowing the user either to avoid unphysical values, or to constrain some parameters to lie within some experimental uncertainty. Finally, at the bottom of the panel, the user can select the fitting algorithm -see section 6.3 for further details.
When experimental data is imported, it can be locally averaged ("binned") or a moving-average smoothing applied using the Data Processing menu option. Data binning is recommended when the number of experimental data points is large ( 5000), since the computation time scales roughly linearly with the data length. After computing the spectra, the data can be exported, either by saving the plot as an image (in any of matplotlib's supported formats: png, ps, eps, pdf, tiff, svg), or by saving the calculated data as a csv file.

Methods file
For integrating into other Python scripts, the elecsus methods.py module allows ElecSus to be called using a functional approach, for either calculation of spectra or fitting data using the calculate() and fit data() methods, respectively.
Experimental parameters are passed to these methods as key:value pairs in a Python dictionary (a list of keys can be found in the code comments). This change has the advantage that parameters can be passed in any order, and unspecified parameters use default values, reducing the complexity of code needed. is produced across the cell by this configuration, which then yields the filter profile shown by the blue solid line on the bottom panel. The non-uniform field filter compares well to the filter profile used in previous work [9] (dashed black line), with largely similar features and a slightly higher peak transmission value.

Update to the fitting methods
Though conceptually the same as the previous version, the implementation of data fitting has been updated for the new version, to make the code more clear and also to make use of new fitting options which are possible using the lmfit module [27]. This module natively supports the ability to fix or vary fit parameters, which greatly simplifies the coding required for a many-parameter fit where not all parameters are allowed to vary. lmfit also allows the user to specify bounds on parameters, which can prevent unphysical values (e.g. negative cell length) or narrow the parameter range when experimental details are known to a good level of accuracy.
The three algorithms from the previous version, Marquardt-Levenburg (ML), Random-Restart (RR) and Simulated Annealing (SA) are retained, and their functionality is the same. In addition to the above algorithms, Differential Evolution (DE) has been added as an option in the GUI, which is a global fitting routine based on stochastic methods developed by Storn and Price [59] that is reported to converge quicker than the SA (Metropolis algorithm [60]) method. In principle, any of the methods that are supported by the lmfit module can be used. However, it is beyond the scope of this work to detail their individual advantages and disadvantages; more information can be found on the scipy documentation pages [61]. The ML  method [62] should be used for the simplest fit problems with few varying parameters. It is the quickest algorithm, but can only find local minima which is an issue for complex problems with a rich parameter space. In these cases, either DE, RR or SA should be used, which are more likely to find the global minimum of the parameter space.

Additional changes since version 1
Owing to the numerous additional features, the previous runcard.py way of using ElecSus is now unsupported, and hence backwards compatibility is broken with version 1. A full list of minor changes and updates can be found on the GitHub page for ElecSus.

Installation and usage
The program is hosted on GitHub at www.github.com/jameskeaveney/ElecSus and the program can be downloaded directly from there either by using git clone if git is installed, or alternately as a zip archive. Installation as a Python module is optional, but can be done by running python setup.py install in a command-line/terminal from the top directory. The GUI can be run from the command-line via python elecsus gui.py Figure 11: Screenshot of the graphical user interface to ElecSus. The left side of the frame shows the calculated spectrum/spectra, whilst the right side contains the experimental parameters to simulate. from the elecsus sub-directory. For integration into user Python scripts, we provide elecsus methods.py which includes two functions: calculate() and fit data(). These take in parameters as Python dictionaries (see source code docstrings for lists of parameter keys), and output a series of numpy arrays which contain the spectral data, and fit parameters with associated uncertainties in the case of fit data().

Test data
Along with the program, we provide another GitHub repository which comprises a series of test data:

www.github.com/jameskeaveney/ElecSusTestData
This test data includes the two examples from the previous version of ElecSus, and also includes all normalised experimental data from the figures in this paper. In the appendix we provide initial parameters for fitting ElecSus to these example data sets.

Conclusions and outlook
We have presented an updated computer program to calculate the electric susceptibility of an alkali-metal vapour. In addition to the previous features of ElecSus (versions 1 and 2), the program is now able to account for magnetic fields with arbitrary orientation with respect to the light propagation axis, and electric field propagation. Together, these allow calculation of susceptibility through non-uniform samples (e.g. magnetic field or density gradients), and the inclusion of optical elements such as birefringent windows. For each of these major changes, we have demonstrated their applications with comparison to example data sets, and find excellent agreement in all cases. In addition, we have developed a graphical interface and new API, which greatly simplifies the workflow for the majority of applications, which we hope will allow ElecSus to reach a wider audience and be more useful to the wider atomic physics community. Quarter-waveplate with fast axis at angle θ to x-axis e iπ/4 cos 2 (θ) + i sin 2 (θ) (1 − i) sin(θ) cos((θ) (1 − i) sin(θ) cos(θ) sin 2 (θ) + i cos 2 (θ) Half-waveplate with fast axis at angle θ to x-axis cos(2θ) sin(2θ) sin(2θ) − cos(2θ) Birefringent material that imprints a phase shift φ oriented with the fast optical axis at an angle θ to x-axis e −iφ/2 e iφ/2 cos 2 (θ) + e −iφ/2 sin 2 (θ) (e iφ/2 − e −iφ/2 ) cos(θ) sin(θ) (e iφ/2 − e −iφ/2 ) cos(θ) sin(θ) e iφ/2 sin 2 (θ) + e −iφ/2 cos 2 (θ)