ARC: An open-source library for calculating properties of alkali Rydberg atoms

We present an object-oriented Python library for computation of properties of highly-excited Rydberg states of alkali atoms. These include single-body effects such as dipole matrix elements, excited-state lifetimes (radiative and black-body limited) and Stark maps of atoms in external electric fields, as well as two-atom interaction potentials accounting for dipole and quadrupole coupling effects valid at both long and short range for arbitrary placement of the atomic dipoles. The package is cross-referenced to precise measurements of atomic energy levels and features extensive documentation to facilitate rapid upgrade or expansion by users. This library has direct application in the field of quantum information and quantum optics which exploit the strong Rydberg dipolar interactions for two-qubit gates, robust atom-light interfaces and simulating quantum many-body physics, as well as the field of metrology using Rydberg atoms as precise microwave electrometers.


Introduction
Highly-excited atoms, in so called Rydberg states, provide strong atom-atom interactions, and large optical nonlinearities.
✩ This paper and its associated computer program are available via the Computer Physics Communication homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655).
They are a flourishing field for quantum information processing [1,2] and quantum optics [3][4][5] in the few to single excitation regime, as well as many-body physics [6][7][8][9][10][11], in the manyexcitations limit.Their exploration requires detailed knowledge of both the single-atom properties, such as lifetimes, energies and transition dipoles elements, as well as atom pair properties, such as their interactions and strongly perturbed energy levels for atoms at distances in the micrometre range.Rydberg states are also highly sensitive to external DC and AC fields, making them excellent candidates for precision electrometry and imaging in the microwave [12,13] and terahertz [14] range, as well as performing state engineering to tune pair-state interaction potentials [15][16][17].Although many of the relevant calculations share the same primitives, such as numeric integration of atomic wavefunctions based on measured energy levels and model core potentials, these basic efforts have been repeated by many groups independently so far.To date no single common resource has emerged for building complex calculations, or for performing quick numerical estimates.To this end, we have developed Alkali Rydberg Calculator (ARC) [18], a library of routines written in Python, using objectoriented programming to make a modular, reusable and extendable collection of routines and data for performing useful calculations of single-atom and two-atom properties, like level diagrams, interactions and transition strengths for alkali-metal atoms.
The hierarchical nature of the package helps organise possible calculations into abstraction levels, allowing one to pick the information at the relevant level for the calculation at hand.The data for individual atomic species is provided as classes that inherit calculation methods defined as abstract classes, allowing one easily to check and update relevant data, should future measurements improve some of the experimentally estimated values.Detailed documentation is provided for all the ARC's modules [19].In addition, the code is commented, cross-referenced in-line and uses self-descriptive names.Whenever possible, the class and function names reflect the hierarchical structure of atomic physics knowledge and the natural decompositions of the calculations, not the low-level implementation details.In addition to the documentation, ARC has example snippets provided as an IPython notebook [20], giving an overview of Rydberg physics and how to perform calculations using the package library.This is a good starting point for new users.
To facilitate the initial adoption of the package and to allow quick calculations useful in the research planning stages, we are also providing a web-interface to basic functions of the package [21].This allows any device with a web browser to access the web-server, that will use the ARC package to perform the calculations and obtain results that are transferred back to the users' web browser.In the process, the web service self-generates and outputs the code, so it can be used as an example-on-demand service, providing a starting point for more complex calculations.
This paper is organised as follows.An overview of the ARC architecture is presented in Section 2, where we introduce the theoretical framework for performing Rydberg state calculations, e.g.calculating wavefunctions and diagonalising interaction Hamiltonians.Here we also provide illustrative examples for building up calculations and visualising results with the provided tools.Initial setup for the ARC library is presented in Section 3, and specific details relating to the implementation are discussed in Section 4. Finally, Section 5 briefly summarises the package and outlines future possible expansions for the library, with a complete list of ARC classes, methods and functions detailed in Appendix A. Detailed documentation of the module is provided in .htmlformat in the Supplemental Material [19] or available from the ARC website [21], along with an IPython notebook [20] that contains numerous additional examples and code snippets that reproduce many of the results from the literature.

Overview
The structure of the ARC library is shown in Fig. 1.At the lowest level, the wigner module implements angular momentum algebra (Wigner 3-j and 6-j coefficients and the WignerD rotation matrix), and arc_c_extensions is a Python extension coded in C to provide fast calculation of the radial part of the atomic wavefunctions.On a higher level, alkali_atom_functions uses these low-level functions to build general methods for calculating single-atom properties, which are contained within the abstract class AlkaliAtom that implements calculation of dipole matrix elements, transition rates, energy levels etc.The module alkali_atom_data defines an explicit class for each alkali element (e.g.Rubidium87(),Caesium()) that encodes all relevant physical parameters and inherits the calculation methods from the parent AlkaliAtom class.These atom classes can be passed as arguments to either of the core calculation modules, calculations_atom_single that implements interactive energy level diagrams (LevelPlot) and calculates Stark maps for atoms in external fields (StarkMap), or calculations_atom_pairstate for dealing with two-atom effects such as long-range dipole-dipole interactions.This pair-state module provides a sophisticated interface to automatically identify Förster resonances for atoms in weak electric fields (StarkMapResonances) and calculate atomic interaction potentials at both long and short range including up to quadrupole-quadrupole couplings (PairStateInteractions).In the following we will outline the basic functionality of the ARC module, provide the theoretical framework for the various modules and give examples of relevant calculations implemented in the library.

The AlkaliAtom class
Almost all calculated quantities in the ARC library can be derived from knowledge of the Rydberg state energy levels and matrix elements.The functions to calculate these values along with other primitive single-atom properties are encapsulated within the methods of the abstract class AlkaliAtom defined by the module al-kali_atom_functions.This class is used in alkali_atom_data where each alkali-metal atom (Lithium6, Lithium7, Sodium, Potassium39, Potassium40, Potassium41, Rubidium85, Rubidium87 and Caesium, as well as Hydrogen) inherits calculation methods from this abstract class, and specifies all the necessary numerical values for a given atom.Calculations with the ARC library begin by initiating an atom class associated with the relevant atom, as shown in the example below for caesium.Physical properties can then be determined using the methods in the form atom.functionName(parameters), where a complete list of available methods is listed in Table B.2 and documentation [19].
The sections below outline the key properties.

Rydberg atom energy levels
Energy of the Rydberg state with principal quantum number n and orbital and total angular momentum ℓ and j respectively, are calculated from the Rydberg formula where E IP is the ionisation energy threshold (ionisationEnergy), Ry is the Rydberg constant corrected for the reduced mass (scaledRydbergConstant), Ry = m e /(m + m e )Ry ∞ and δ n,ℓ,j is the quantum defect (getQuantumDefect) given by where δ 0 , δ 2 , . . .are modified Rydberg-Ritz coefficients obtained by fitting precise measurements of the atomic energy levels [22].
Energies E for a given Rydberg state relative to the ionisation limit can be obtained using getEnergy (in eV), as well as methods getTransitionWavelength and getTransitionFrequency to return the wavelength (m) or frequency (Hz) of a transition between two Rydberg states.Energies and transitions are given relative to the centres of gravity of the hyperfine split ground and excited states.Overview of the Alkali Rydberg Calculator (ARC) package for Python.Object-oriented structure with hierarchy reflecting the structure of atomic physics calculations is used.This allows the user to choose the abstraction level at which one wants to work, from low-level implementations of basic functions of alkali_atom_functions.AlkaliAtom for finding dipole matrix elements and lifetimes, to high-level functions that allow automatic Förster resonance finding and exploration of complex energy-level diagrams of atomic pair-states at small inter-atomic separations in calculations_atom_pairstate.

Rydberg atom wavefunctions
In order to calculate matrix elements for electric dipole and quadrupole couplings of the Rydberg states, it is necessary to calculate the radial wavefunctions R(r) by numerically integrating the Radial Schrödinger equation for the valence electron.Using the substitution ρ(r) = rR(r), we can remove the first order differential to obtain an equation in the form: where V (r) is the spherically symmetric central potential (potential).For hydrogen, and high orbital angular momentum states ℓ > 3, this is simply the Coulomb potential of V (r) = −1/r + V so (r), with added (relativistic) spin-orbit interaction V so (r) = α L • S/(2r 3 ), where α is the fine structure constant.However for alkali atoms with closed shells it is necessary to introduce a model potential that gives a Coulomb potential at long-range and at short range accounts for effects of the core penetration of the valence electron.We adopt the core potential (corePotential) of Marinescu et al. [23] given by the form where α c is the core polarisability, r c is a cutoff radius introduced to avoid divergence of the polarisation potential at short range and the radial charge (effectiveCharge) is given by with ℓ-dependent parameters a 1..4 and r c taken from Table 1 of Ref. [23] which were obtained by fitting the model potential to measured energy levels for each element.Using this model potential and the known Rydberg energy from above, the radial wavefunctions can be calculated by numerically integrating Eq. (3), as shown on Fig. 2.This is achieved by performing a transformation to integrate the function X (r) = R(r)r 3/4  in terms of the scaled co-ordinate x = √ r [24] that gives an approximately constant number of points across each period of oscillation in the wavefunction.The result is a second-order differential equation of the form which is efficiently solved using the Numerov method [25,26]  [23], using state energy from quantum defects or NIST ASD database [28].These calculations are the starting step in calculating all atom-light couplings (dipolar and higher order).
where h is the step size, T (x) = h 2 g(x)/12 and It is necessary to truncate the range of integration as at short range the model becomes unphysical and diverges, whilst at longrange the wavefunction decays to zero.Following Ref. [27], the limits of integration are set to use an inner radius of r i = 3  √ α c , and an outer radius of r o = 2n(n + 15) which is much larger than the classical turning point of the wavefunction.To minimise errors introduced by the approximate model potential at short range, the integration is performed inwards, starting at r o .For high orbital momentum states (ℓ > 3) divergence can occur even before r i , which is automatically detected and integration is stopped at the closest wavefunction node before divergence occurred.The wavefunctions are then normalised, and returned by radialWavefunction.

Matrix elements
Relevant properties for alkali atoms such as lifetimes (Section 2.2.5),Stark shifts (Section 2.3.2) and atomic interaction potentials (Section 2.4.1)require evaluation of the electric dipole and electric quadrupole matrix elements from state |n, ℓ, m ℓ ⟩ to |n ′ , ℓ ′ , m ′ ℓ ⟩.For electric dipole transitions, the interaction is dependent upon matrix elements of the form H = −er • ε where ε is the electric field polarisation vector.Expanding the operator in the spherical basis and using the fact that the calculation can be separated into radial overlap, and angular overlap of the electron wavefunction with operator H, the resulting matrix element can be evaluated in terms of angular momentum terms and radial matrix elements using the Wigner-Ekart theorem [29] ⟨n, ℓ, where q represents the electric field polarisation (q = ±1, 0 driving σ ± , π transitions respectively), and the braces denote a Wigner-3j symbol (Wigner3j).We use Condon-Shortley phase convention [30] for spherical harmonics.The reduced matrix element ⟨ℓ||r||ℓ ′ ⟩ (getReducedMatrixElementL) is given by where the radial matrix element (getRadialMatrixElement) is evaluated from using numerical integration of the calculated wavefunctions.
Transforming these to the fine-structure basis Eq. ( 9) can be expressed in terms of states j, m j (getDipoleMatrixElement) as where curly braces denote a Wigner-6j symbol (Wigner6j).
To account for higher order multipole moments in the interaction between atom pairs (see Section 2.4.1), it is also necessary to calculate quadrupole matrix elements (getQuadrupoleMatrix Element) of the form As with the quantum defect model for the energies, these numerical approaches provide accurate estimates of the dipole and quadrupole terms for highly-excited states but have a large error for low-lying states where the electron density is weighted close to the atomic core where integration is most sensitive to the divergence of the model potential.To overcome this limitation, values for dipole matrix elements available in the literature either from direct measurement or more complex coupled-cluster calculations [31,32] are tabulated (accessible through the function getLiteratureDME).Before calculating a matrix element, the ARC library first checks if a literature value exists and if so utilises the tabulated value with the smallest estimated error.Otherwise a numerical integration is performed using the method outlined above.For each element, these tabulated values are contained within an easily readable .csvfile (file name specified in literatureD-MEfilename) that lists the value along with relevant bibliographical reference information, making it easy for users to add new values at a later date.

Rabi frequency
An important parameter in experiments with Rydberg atoms is the Rabi frequency Ω = d • E/h where d is the dipole matrix element for the transition and E is the electric field of the laser driving the transition.For a transition from state |n 1 ℓ 1 j 1 m j 1 ⟩ → |n 2 ℓ 2 j 2 m j 2 ⟩ using a laser with power P and 1/e 2 beam radius w, the Rabi frequency (in rad s −1 ) can be obtained as atom.getRabiFrequency(n1,l1,j1,mj1,n2,l2,j2,mj2,P,w) The related function atom.getRabiFrequency2returns the Rabi frequency calculated from the electric field amplitude.

Excited state lifetimes
Radiative lifetimes of alkali atoms can be calculated using dipole matrix elements to determine the Einstein A-coefficient [33] for transitions (see getTransitionRate), where ω nn ′ is the frequency of the transition from |nℓ⟩ to |n ′ ℓ⟩.The total radiative decay Γ 0 is then obtained by summing over all dipole-coupled states |n ′ ℓ ′ ⟩ of lower energy, For an environment at finite temperature, it is necessary to account for the interaction of the atom with black-body radiation (BBR) [34].This causes stimulated emission and absorption which depend on the effective number of BBR photons per mode, nω , at temperature T given by the Planck distribution nω = multiplied by the A-coefficient resulting in a BBR transition rate where the summation n ′ , ℓ ′ includes also states higher in energy, since BBR can drive these transitions.Finally, the effective lifetime Fig. 3 shows the relative contribution to the excited-state lifetime of Rb 30 S 1/2 due to radiative decay (dominated by low-lying states due to the ω 3 scaling in Eq. ( 14)) and black-body decay at 300 K calculated using ARC (for full code see Appendix A).Comparison of our calculated lifetimes [20] with previous work on radiative [33] and BBR induced lifetimes for alkali-metals [34] yields excellent agreement.

Atomic vapour properties
For experiments using thermal vapours of alkali atoms, a number of useful functions are provided for returning atomic vapour pressure (getPressure) or number density (getNumber Density), as well as the average interatomic distance (getAverage InteratomicSpacing) and atomic speed (getAverageSpeed) at a given temperature.For full listing of functions see Table B.2 and full documentation in Supplemental [19].

Single-atom calculations
The module calculations_atom_single provides calculations of single-atom properties.It supports plotting of atom energy levels with interactive finding of associated transition wavelengths and frequencies (Section 2.3.1),i.e.Grotrian diagrams, and atomenergy-level shifts in electric fields (Section 2.3.2),i.e.Stark maps.

Level plots
The LevelPlot class facilitates simple plotting of atomic energy levels.It provides an interactive plot for exploring transition wavelengths and frequencies.For example, generation and plotting of the caesium energy level diagram, including ℓ states from S to D for principal quantum numbers from n = 6 to n = 60, can be realised with: atom = Caesium() levels = LevelPlot(atom) # parameters: nmin, nmax, lmin, lmax levels .makeLevels(6,60,0,3)levels .drawLevels()levels .showPlot()This simple example also demonstrates that the more complicated calculations are implemented as a compact classes, whose initialisation and methods closely follow naming and stages one would perform in a manual calculation.Yet, working on high abstraction level, one can obtain information directly relevant for research in just a couple of high-level commands [20].

Stark shifts
Calculation of atomic Stark shifts in external static electric fields provides a tool for both precision electrometry and a mechanism for tuning interatomic interactions to a Förster resonance to exploit strong resonant dipole-dipole interactions (Section 2.4.4).To find the energy of the atom in an electric field E applied along the z-axis it is necessary to find the eigenvalues of the system described by the Stark Hamiltonian where H 0 is the Hamiltonian for the unperturbed atomic energy levels and E is the applied electric field.The electric field term causes a mixing of the bare atomic energy levels due to coupling by the Stark interaction matrix elements ⟨n, ℓ, j, m j |E ẑ|n ′ , ℓ ′ , j ′ , m ′ j ⟩ which can be evaluated from Eq. ( 12) with q = 0.The selection rules of the Stark Hamiltonian give ∆m j = 0, ∆ℓ = ±1 such that only states with projection m j are coupled together, so each Stark map can be constructed by taking basis states with the same m j value.Following the method of Zimmerman et al. [27], Stark shifts are calculated by exact diagonalisation of the Hamiltonian.
Stark shift calculations are handled using the StarkMap class, which is initialised by passing the appropriate atom class.For a target state |n, ℓ, j, m j ⟩, a range of n values from n min to n max with values of ℓ up to ℓ max is required to define the basis states for the Stark Hamiltonian.Convergence is typically achieved using ℓ max of 20 and n max − n min ∼ 10, however for large applied fields or higher values of n it will be necessary to increase the basis size to account for the strong mixing of the energy levels.Finally, the Hamiltonian is diagonalised for each value of the electric field (defined in V/cm).As an example, to calculate the Stark map shown in Fig. 4 for the 28 S 1/2 m j = 1/2 state in Cs we include states n min = 23 to n max = 32 with l max = 20 for 600 equidistant electric field values in the range from 0 to 600 V/cm.The corresponding program call is: # parameters: n, l , j , mj, nmin, nmax, lmax calc .defineBasis (28, 0, 0.5, 0.5, 23, 32, 20) calc .diagonalise (np.linspace (0.,600*1e2,600)) calc .plotLevelDiagram()calc .showPlot() In interactive mode, the plot can be clicked to obtain the dominant contribution of the basis states within each eigenstate and by default the figure is highlighted in proportion to the fractional contribution of the target state (in this case 28 S 1/2 m j = 1/2).

Pair-state calculations
The module calculations_atom_pairstate contains classes and methods for the calculation and visualisation of long-range and short range interactions (PairStateInteractions), as well as an automated tool for identifying Förster resonances (StarkMapResonances) for electric field tuning of the interaction potential.

Interatomic interactions
Pair-wise interactions between two atoms with internuclear separation R, and electron coordinates r 1 and r 2 relative to the respective nuclei, can be expanded in multipole form as [35] V where spectively to dipole-dipole, dipole-quadrupole and quadrupolequadrupole interactions, and where ) is the binomial coefficient and Y L,m (r) spherical harmonics.In Eq. ( 20) the quantisation axis is oriented along internuclear axis R. On the other hand, the atomic states' quantisation axis is typically defined with respect to the driving laser, being directed along the laser propagation direction for circularly polarised laser beams, or in plane of the electric field vector, perpendicular to the laser propagation direction, for linearly polarised driving, or in the direction of applied static electric and magnetic fields.When the atomic quantisation axis is along the internuclear axis, matrix (wignerDmatrix), so that inter-atomic axis defines the quantisation direction [29].The coupling can then be easily calculated between the rotated states |s⟩ = w D (θ , φ)|j, m j ⟩.Note that |j, m j ⟩ in the rotated basis |j ′ , m ′ j ⟩ is now, in general, a superposition of m ′ j states.This multipole expansion is valid as long as the wavefunctions of the atoms do not overlap, which is the case for interatomic separations R larger than the Le Roy radius [36] Evaluation of the Le Roy radius can be achieved using the function getLeRoyRadius.For example the Le Roy radius for two Caesium atoms in the n S state is 0.1 µm for n ∼ 20 and reaches 1 µm for n ∼ 60, marking the interatomic distance below which the results of the pair-state diagonalisation become invalid.
To understand the effect of interactions, consider a pair of atoms in Rydberg pair-state |r, r⟩ coupled via V (R) to Rydberg pair-state |r ′ , r ′′ ⟩ which has an energy defect ∆ r ′ ,r ′′ = 2E r − E r ′ − E r ′′ .In the pair-state basis {|rr⟩, |r ′ r ′′ ⟩} their interaction is described with the Hamiltonian ) . ( We see that at short range where V (R) ≫ ∆ r ′ ,r ′′ the splitting of the energy eigen-states ±V (R) is dominated by the pair-state interaction energy.Assuming that V(R) has non-zero dipole-dipole term of the multipole expansion [L 1 + L 2 = 2 in Eq. ( 20)], this corresponds to the resonant dipole-dipole regime, giving eigenstates' energy distance dependence ∝ C 3 /R 3 .At large R where ∆ r ′ ,r ′′ ≫ V (R), the interaction is second order leading to an energy shift of −V (R) 2 /∆ r ′ ,r ′′ = −C 6 /R 6 , known as the van der Waals regime.The cross-over between these regimes occurs at the van der Waals radius R vdw where For a real system, whilst the pair-state with the smallest ∆ r ′ ,r ′′ dominates the resulting interaction shift, it is necessary to consider the effects of all near-resonant pair-states that are coupled by the V (R) interaction term above.Calculations of these pair-wise interaction potentials are handled by the PairStateInteractions class that is initialised by specifying the element name and target pair-state n 1 , l 1 , j 1 , n 2 , l 2 , j 2 , m j 1 , m j 2 whose behaviour we want to explore.

Long-range limit
At large separation, for off-resonantly (∆ r ′ ,r ′′ ̸ = 0) coupled states, the dipole-dipole interactions dominate to give an interaction of the form −C 6 /R 6 van der Waals potential where the sign depends on the energy defect of the closest dipole-coupled pair-states |r ′ r ′′ ⟩, leading to attractive or repulsive interactions accordingly.In the long-range limit V (R) ≪ ∆ (where ∆ is the energy defect of the closest dipole-coupled pair-state), the C 6 coefficient for pair-state |rr⟩ can be evaluated using second-order perturbation theory as where the sum runs over all pair-states |r ′ r ′′ ⟩ whose energy differs from the pair-states |rr⟩ energy for ∆ r ′ ,r ′′ < ∆ E , where ∆ E is some maximal energy defect that provides a truncation of the basis states.This calculation is performed using the method getC6perturbatively, returning C 6 in units of GHz µm 6 .Users must specify the relative orientation of the atoms, the range of principal quantum numbers for the states used in calculation, and the maximal energy defect.For example for the interaction of two rubidium atoms in 60 S 1/2 , m j 1 = 1/2, m j 2 = −1/2 states, whose inter-atomic axis is set at an angle of θ = π/6 with respect to the quantisation axis [Fig.5(a)], the C 6 interaction term can be calculated from states with principal quantum number differing maximally δn = 5 from the n = 60, and maximal energy difference in pair-state energies of h × 25 GHz as # parameters: atom, n1, l1, j1 , n2, l2 , j2 , mj1, mj2 calc = PairStateInteractions (Rubidium(), 60, 0, 0.5, 60, 0, 0.5, 0.5, -0.5) \ # parameters: theta, phi, deltan, deltaE c6 = calc.getC6perturbatively(pi/6,0, 5, 25.e9) ; print("C_6 = %.0fGHz (mu m)^6" % c6) Using this function, the anisotropy of the V (R) interaction can be easily identified as shown in Fig. 5(b) which plots the magnitude of C 6 for a pair of atoms in the 60 D 5/2 state of Rubidium for θ = 0 . . .2π .

Exact interaction potential
To evaluate the interaction potential for arbitrary separation [valid down to the Le Roy radius, Eq. ( 21)], it is necessary to diagonalise the matrix V (R) containing all the interatomic couplings for each separation R to obtain exact values of the eigenvalues and eigenstates describing the interaction between atomic pair-states.Due to the strong admixing of states, this process requires careful choice of atomic basis states, including higher order orbital angular momentum states up to ℓ max .For each distance R, the matrix is diagonalised using an efficient ARPACK package provided through Numpy, and the n eig eigenstates whose eigenvalues are closest to the target pair-state are returned.
Fig. 6 provides an example of the interaction potential for a pair of atoms initially in the |60 S 1/2 m j = 1/2, 60 S 1/2 m j = −1/2⟩ As with the StarkMap method, the figure is interactive allowing users to determine the dominant composition (expressed in the pair-state basis) of a given eigenstate.The default behaviour is to highlight the eigenvalues proportional to the fractional contribution of the original target pair-state.Alternatively, as shown on Fig. 6], the optional parameter drivingFromState is used in the call to diagnolise to give highlighting in proportion to the relative laser coupling strength from a given state, assuming that one atom is already in one of the two states that make the initially specified pair-state.
By default, PairStateInteractions includes only dipole coupling between the states in calculating level diagrams.Interactions up to quadrupole-quadrupole term [including all terms L 1 + L 2 ≤ 4 in Eq. ( 20)] can be included by setting optional parameter interaction-sUpTo = 2 during the PairStateIntearctions initialisation, which can be important for the short-distance structure of level diagram [37].For evaluation of long-range potential curves this additional term makes only small perturbations to the asymptotic C 6 behaviour.However, for accurate determination of the molecular levels of the short-range potential wells it is important to include the higher order multipole terms.Convergence should also be checked by increasing the basis size and re-evaluating the relevant parameters to ensure all the relevant states are included in the calculation basis.
Following diagonalisation of the pair-state interaction matrix, the long-range (C 6 ) and short range (C 3 ) dispersion coefficients can be evaluated using the methods getC3fromLevelDiagram and getC6fromLevelDiagram respectively, which perform a fit to the eigen-energies associated with the state containing the largest admixture of the target state.A method for finding the crossover distance between van der Waals and resonant dipole-dipole interactions, i.e. van der Waals radius R vdw , is also provided (getVdwFromLevelDiagram).

Stark tuned Förster resonances
As outlined in Section 2.4.1, the finite pair-state energy defects ∆ associated with the closest dipole-coupled pair-states lead to a transition from first order resonant dipole-dipole interactions (∝ R −3 ) to second order van der Waals (∝ R −6 ) at the van der Waals radius.Using external electric fields however, it is possible to Stark shift the pair-states into resonance to obtain longrange ∝ R −3 interactions for all values of R, known as a Förster resonance [38][39][40][41].
To identify suitable Förster resonances the StarkMapResonances class is used, taking in a pair of target pair-states and performing diagonalisation of the Stark Hamiltonian of Eq. ( 18) in the pairstate basis.Due to the angular momentum selection rules of the dipole operator V (R), the target pair-state can be dipole-coupled to pair-states with ∆m j = ±1, 0, it is necessary to calculate Stark maps for a range of m j manifolds.Following diagonalisation, only pair-states with energy close to the target state are considered, with any states not dipole-coupled to the target pair-state being discarded.As the electric field coupling leads to significant mixing of the zero-field pair-states, the algorithm identifies the basis state containing the largest target state fraction at a given electric field to test for the closest dipole-coupled pair-states.Finally, an interactive plot routine enables users to identify states that have a Förster In the above example, similarly as when finding a Stark map (Section 2.3.2),we have to specify a range of acceptable principal quantum numbers [n min , n max ] and maximal orbital angular momentum l max for the basis states for Stark map calculations, as well as the electric field range (in this case 0-1 V/cm at 200 equidistant points).An additional argument energyRange defines the energy window within which we will keep the resonant pair-states.In the above given example, this is in range h × [−0.8, 3.0] GHz.From the interactive plot, selection of the pair-state eigenvalues shows a Förster resonance occurring at 0.18 V/cm with the 42 F , 46 P 3/2 pair-state where, due to the electric coupling, the F state is an admixture of the 42 F 5/2,7/2 m j = 3/2 states.

Installation and usage
It is assumed the pre-requisites (Numpy, SciPy and Matplotlib) are installed and can be located by the Python interpreter (see e.g.[42]).Both Python 2.7 and 3.5 are supported.To achieve good performance, it is recommended to use Numpy packages that connect to optimised backends, like ATLAS [43].Prepackaged Python distributions, like Anaconda [42], provide this out-of-the box.The ARC library can be downloaded online [18] as a .zipfile release.Installation is performed by extracting the downloaded .ziparchive and copying the arc subfolder into the root of your project directory.It is important that Python has write access to the folder where the package is located, so that database files (stored in arc/data/) can be updated and used.

Optimised Numerov integrator
Integration of the atomic wavefunctions to obtain dipole matrix elements is numerically intensive and by default ARC uses an optimised Numerov integration routine implemented with C extensions for Python (arc_c_extensions).In the unlikely case the precompiled executable arc_c_extensions.so is incompatible with the installed system, please install a C compiler and compile arc_c_extensions.clocated in the ARC root folder using by calling from command line python setupc.pybuild_ext --inplace Note that path to arc directory should not contain spaces in order for setupc.pyscript to work.Alternatively, the native Python solver can be used by setting the optional argument cpp_numerov=False when initialising the atom class, however this is not recommended for intensive calculations as it results in substantially slower execution.

Getting started
To initialise the library, use the following code at the start of a Python script or interactive IPython notebook with # locate ARC Directory import sys, os rootDir = '/path/to/root/directory/for/arc' sys .path.append(rootDir) os.chdir(rootDir) # import ARC library from arc import * This firsts sets a path to the directory containing ARC package directory on your computer. 1This is the recommended way of using the package in a research environment, since users can easily access, check and change the underlying code and constants according to their needs.ARC is now ready for use, and can be tested using the example code above.Numerous additional examples are provided in IPython notebook ''An Introduction to Rydberg atoms with ARC'' [20].

Physical constants
As mentioned above, the atomic properties are encapsulated in classes which contain the relevant atomic properties determined from the literature along with model potential coefficients taken from Marinescu et al. [23] which have been optimised against measured energy levels.For each atom, asymptotic expansions of the quantum defects are used to determine the energy levels of states with high principal quantum number to high accuracy, using measured values of the ionisation energy, Rydberg constant and quantum defects taken from Li [44,22], Na [22], K [22], Rb [45][46][47][48], and Cs [49][50][51].For the low-lying states, the quantum defects do not accurately reproduce the measured energies and instead energy levels are determined from data in the NIST ASD database [28].By default, the cut-off between tabulated and calculated energies is determined by the point at which the error in the calculated energy exceeds 0.02%.

Tracking calculation progress
For tracking progress of calculations, most functions (see documentation [19] for details) provide optional bool arguments, that turn on printing of additional information about calculations.For example, progressOutput prints basic status information, and is recommended for everyday use.If debugOutput is set to True, more verbose information about basis states and couplings will be printed in the standard output.

Memoization
In order to achieve good performance, memoization is used throughout.That is, when functions receive request for calculation, memorised previous results are checked, and the value is retrieved from memory if it already exists.Both an in-application SQL database (SQLite) and standard arrays are used for this.Memoization is used for dipole and quadrupole matrix element calculations and angular coupling factors, that are independent of principal quantum number of the considered state, including Wignern J coefficients and WignerD matrices.Single-atom and pair-state calculations will automatically update databases if new values are encountered, speeding up the future calculations.Updating dipole and quadrupole matrix element database can be done manually too, by calling updateDipoleMatrixElementsFile.

Exporting data to file
If one wants to use the obtained results in another program, for further analysis and calculation, both PairStateInteractions and StarkMap provide exportData method that allows saving data in .csvformat, that is human-readable and easy to import in other software tools.As a commented header of the exported files, ARC will record details of the calculation parameters, in humanreadable form.Provided that the initial calculation is variable calc, its data is exported by calling calc .exportData("rootFileName") The program will output list of files (typically three files) storing the calculation data, whose names will be starting with "rootFile-Name".

Advanced interfacing of ARC with other projects
In addition to graphical exploration, and exporting relevant data as .csvfiles, advanced users incorporating ARC in their own projects might want to access directly Stark or pair-state interaction matrices, and corresponding basis states that are generated by defineBasis methods of the corresponding calculation classes.Basis (pair-) states are accessible for both methods in basisStates array.Stark matrix can be assembled as the sum of the diagonal matrix mat1 recording state energies, and the off-diagonal matrix that depends on applied field and can be obtained as product of electric field and mat2.For pair-state interactions, the interaction matrix can be obtained as a sum of interatomic-distance independent matrix matDiagonal, recording energy defects of pair-states, and distance dependent matrices stored in matR, where matR[0], matR [1] and matR [2] store respectively dipole-dipole, dipole-quadrupole and quadrupole-quadrupole interaction coefficients for interaction terms [Eq.( 20)] scaling respectively with distance R as ∝ R −3 , ∝ R −4 and ∝ R −5 .For more details about this, examples of use, and other accessible calculated information, we refer the reader to the detailed documentation [19].

Outlook
This paper describes version 1.2 of ARC that can be downloaded from [21].Supplementary information provides full documentation [19] for functions listed in Appendix B, as well as an IPython notebook [20] that contains more elaborate examples, code snippets and benchmarks against published results.This notebook provides an introduction to most of the capabilities of the ARC, as well as a good starting point for users, providing useful code snippets.Numerous examples also provide quantitative interactive introduction for anyone starting in the field of Rydberg atomic physics.In addition, for quick estimates, when one is in the lab, conference or other meeting, we provide a web-interface to the package [21].For the latest versions of the package and documentation, please download material from the ARC GitHub page [18].
In the future, the package can be extended to include calculations of dressing potentials [52], magic wavelengths [53], atom-wall interactions [54,55], photoionisation, collisional crosssections [56], tensor polarisability, molecular bound states [57], effects of magnetic fields, microwave tuning of interactions [17] and other atomic properties.These can be added as calculation tools, in separate classes, built on top of the existing library.Alkaline earth elements can be included too.While these were beyond the scope of the current project, we hope that this project can provide an initial seed for a much bigger community project.The code is hosted on GitHub [18], allowing easy community involvement and improvements.A proposal of some basic philosophy for the development is provided in the documentation [19].

Fig. 1 .
Fig.1.Overview of the Alkali Rydberg Calculator (ARC) package for Python.Object-oriented structure with hierarchy reflecting the structure of atomic physics calculations is used.This allows the user to choose the abstraction level at which one wants to work, from low-level implementations of basic functions of alkali_atom_functions.AlkaliAtom for finding dipole matrix elements and lifetimes, to high-level functions that allow automatic Förster resonance finding and exploration of complex energy-level diagrams of atomic pair-states at small inter-atomic separations in calculations_atom_pairstate.

Fig. 2 .
Fig. 2. Electron wavefunctions are calculated by Numerov integration, here shown for caesium 18 S 1/2 m j = 1/2 state.Solving the radial part of the Schrödinger equation in the model potential from[23], using state energy from quantum defects or NIST ASD database[28].These calculations are the starting step in calculating all atom-light couplings (dipolar and higher order).

Fig. 3 .
Fig.3.Spontaneous decays and black-body induced transitions from Rb 30 S 1/2 to n P 1/2,3/2 for environment temperature of 300 K.These transitions have to be included in calculation of excited-state lifetime for Rb 30 S 1/2 state.

Fig. 4 .
Fig. 4. Example of Stark map calculation, showing caesium 28 S 1/2 m j = 1/2 state perturbation by DC electric field.Colour highlights contribution of the |28 S 1/2 m j = 1/2⟩ state in the atom eigenstates |µ⟩.(color online) are easily evaluated.When the quantisation axis is not oriented along R, their relative orientation can be described with polar angle θ and azimuthal angle φ that the inter-atomic axis makes with the quantisation axis [Fig.5(a)].Atom states |j, m j ⟩ are then rotated with WignerD matrices w D (θ , φ)

Fig. 6 .
Fig. 6.Example of interactions between the atoms, causing level shifts in the atomic pair-states.Here shown in the vicinity of Rubidium |60 S 1/2 m j = 1/2, 60 S 1/2 m j = −1/2 state.Highlighting is proportional to the driving strength from the 5 P 3/2 m j = 3/2 state with coupling laser driving σ − transitions.Atom interatomic direction is oriented along the quantisation axis (θ , φ = 0).Using the provided methods we can find van der Waals radius R vdW = 2.4 µ m, and short-range C 3 = 16.8GHz µ m 3 and long range C 6 = 135 GHz µ m 6 dispersion coefficients.(color online)

Fig. 7 .
Fig. 7. Example plot produced by StarkMapResonances for finding Förster resonances for Rubidium pair-state 44 D 5/2 m j = 5/2 shown as grey (red) line.Pairstates whose main admixed state is dipole-coupled to the original pair-state are shown as black lines.Clicking on one of these states, marked with square (blue) on the plot, reveals its composition in the plot title.In this case, it is a mixture of 42 F 5/2 and 42 F 7/2 states in one of the atoms, and essentially unperturbed 46 P 3/2 state that is almost resonant with original pair-state at electric field of E ≈ 0.18 V/cm.(color online)

Table B . 1
Class and function listing of alkali_atom_functions module.

Table B . 2
Methods and function listing of alkali_atom_functions.AlkaliAtom class.Typical relative uncertanties are obtained from comparison to measured values.

Table B . 3
Class listing of alkali_atom_data module.All these classes inherit properties of alkali_atom_functions.AlkaliAtom from Table B.2. Method listing of calculations_atom_single.LevelPlot(atomType) class.

Table B . 8
Function and class listing of wigner module providing support for angular element calculations.